梅花鹿與歐洲馬鹿染色體水平基因組的比較研究

唐麗昕，王天驕，劉華淼，張然然，邢秀梅

(中國(guó)農(nóng)業(yè)科學(xué)院特產(chǎn)研究所 特種經(jīng)濟(jì)動(dòng)物分子生物重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)春 130112)

梅花鹿(Cervusnippon)和馬鹿(Cervuselaphus)是哺乳綱、偶蹄目、鹿科、鹿屬的兩個(gè)物種，其種間分化時(shí)間很短，一般很難得到確切的演化關(guān)系[1]。早在20世紀(jì)80年代，王宗仁和杜若甫[2]就從細(xì)胞遺傳學(xué)角度，通過染色體核型來(lái)解釋梅花鹿和馬鹿的起源演化關(guān)系，指出鹿科動(dòng)物染色體的進(jìn)化方式主要是羅伯遜斷裂[3]，主張馬鹿(2n=68)是從梅花鹿(2n=66)演化而來(lái)的。而Frank等[4]就線粒體全基因組的研究指出，天山馬鹿、東北馬鹿和梅花鹿亞種之間的分化時(shí)間小于0.01 MYA(millions of years ago，距今百萬(wàn)年)，推測(cè)它們可能不是兩個(gè)物種，可能屬于梅花鹿的變種。現(xiàn)有的細(xì)胞遺傳學(xué)和線粒體全基因組方面的研究，并沒有很好的解釋梅花鹿和馬鹿演化關(guān)系的問題。隨著生物技術(shù)的飛速發(fā)展，越來(lái)越多的學(xué)者從全基因組角度出發(fā)來(lái)分析物種進(jìn)化關(guān)系[5-7]，這為梅花鹿和歐洲馬鹿的演化關(guān)系研究提供了新的思路。

比較基因組學(xué)(comparative genomics)是在基因組圖譜和測(cè)序基礎(chǔ)上對(duì)已知的基因和基因組結(jié)構(gòu)進(jìn)行比較，來(lái)了解基因的功能、表達(dá)機(jī)理和物種進(jìn)化的學(xué)科，其涉及物種內(nèi)及物種間的基因組同源性、組成和功能的評(píng)估，并且已用于揭示從基因組進(jìn)化到疾病調(diào)控等多方面的相關(guān)過程[8]。比較基因組學(xué)分為種內(nèi)基因組比較和種間基因組比較兩大類。種內(nèi)基因組比較通常用于進(jìn)行與物種適應(yīng)性等特性有關(guān)的功能基因及特異性位點(diǎn)的研究[9-13]；種間基因組比較可以揭示基因潛在的功能、為闡明物種進(jìn)化關(guān)系及基因組的內(nèi)在結(jié)構(gòu)提供理論基礎(chǔ)[14-20]。

本研究對(duì)梅花鹿和歐洲馬鹿染色體水平基因組進(jìn)行比較基因組學(xué)分析，通過兩個(gè)物種基因組染色體共線性及染色體上基因共線性的分析以及對(duì)染色體倒位涉及的基因進(jìn)行功能富集，揭示了兩個(gè)物種基因組的同源性，將染色體倒位與嗅覺表型進(jìn)行相關(guān)分析，以期為梅花鹿與歐洲馬鹿起源進(jìn)化關(guān)系的研究提供新的理論基礎(chǔ)。

1 材料與方法

1.1 梅花鹿與歐洲馬鹿基因組序列數(shù)據(jù)

梅花鹿的基因組序列來(lái)自本實(shí)驗(yàn)室采用Pacbio三代測(cè)序技術(shù)結(jié)合Hi-C(high-throughput/resolution chromosome conformation capture)輔助組裝技術(shù)組裝得到的基因組序列(暫未發(fā)表)。另外，從NCBI數(shù)據(jù)庫(kù)下載得到歐洲馬鹿的基因組序列Cerla1.0(GCA_002197005.1)[21]。梅花鹿(♀)基因組組裝得到33條染色體序列，包括32條常染色體序列以及1條X性染色體序列，基因組大小為2 481 763 803 bp，其中Scaffold N50的長(zhǎng)度為78 786 809 bp；下載的歐洲馬鹿(♂)基因組得到34條染色體序列，包含33條常染色體序列以及X、Y兩條性染色體序列，基因組大小為3 438 623 608 bp，其中Scaffold N50的長(zhǎng)度為107 358 006 bp。兩個(gè)基因組序列均已組裝至染色體水平(表1)。

表1 梅花鹿和歐洲馬鹿基因組的組裝概況

1.2 梅花鹿與歐洲馬鹿比較基因組學(xué)分析

1.2.1 染色體共線性分析 為了初步獲得梅花鹿與歐洲馬鹿兩個(gè)基因組序列的共線性，本研究使用MUMmer v3.1(http：//mummer.sourceforge.net/)軟件[22]，以梅花鹿基因組為參考，對(duì)梅花鹿和歐洲馬鹿的基因組序列進(jìn)行比對(duì)，參數(shù)采用默認(rèn)，從而獲得兩個(gè)物種間的正、反向最大精確匹配，然后借用show-coords實(shí)用工具獲得比對(duì)的坐標(biāo)位置、百分比等標(biāo)識(shí)，以獲得各條染色體上基因組序列的比對(duì)情況。

同時(shí)，使用標(biāo)準(zhǔn)參數(shù)運(yùn)行Mummerplot程序，得到gnuplot腳本，用于生成繪圖命令，完成基因組序列染色體水平共線性圖譜的繪制。

1.2.2 基因功能富集分析 將梅花鹿與歐洲馬鹿基因組比對(duì)得到的基因組結(jié)構(gòu)變異(染色體倒位)涉及的基因進(jìn)行功能富集分析。使用KOBAS 3.0數(shù)據(jù)庫(kù)[23]以鹿科動(dòng)物的近源物種牛作為背景數(shù)據(jù)集進(jìn)行GO[24]以及KEGG[25]功能富集分析。以明確與這些基因高度相關(guān)的功能和生物學(xué)通路。

1.3 梅花鹿與歐洲馬鹿基因共線性分析

利用共線性分析工具M(jìn)CScanX，獲得梅花鹿和歐洲馬鹿兩個(gè)基因組間基因序列的共線性同源區(qū)。首先使用BLASTP將梅花鹿基因組的蛋白編碼序列建庫(kù)，與歐洲馬鹿基因組全部蛋白質(zhì)序列進(jìn)行all-vs-all比對(duì)[26]，為了避免由于串聯(lián)陣列而產(chǎn)生大量的局部共線基因?qū)Γ珺LASTP的E值默認(rèn)為10-5[27]。全基因組BLASTP結(jié)果用于計(jì)算所有可能的染色體和scaffolds對(duì)的共線性區(qū)塊。接著使用MCScanX程序，將GFF格式文件和BLASTP表格文件作為輸入，根據(jù)所有可能的染色體和scaffold對(duì)的位置，在兩個(gè)轉(zhuǎn)錄方向上對(duì)基因之間的匹配進(jìn)行分類。相鄰的共線基因連接成共線性區(qū)塊，將成對(duì)的共線塊和串聯(lián)基因?qū)Ψ謩e輸出。

同時(shí)，對(duì)基因共線性數(shù)據(jù)統(tǒng)計(jì)得到兩個(gè)基因組間共線性區(qū)域的數(shù)量、共線區(qū)域內(nèi)包含的基因?qū)?shù)量以及梅花鹿與歐洲馬鹿基因組間共線性基因占各自基因組的百分比，利用TBtools[28]軟件輸入兩個(gè)物種同源基因?qū)Φ腃DS序列，計(jì)算出同源基因的同義突變率(Ks)以及Ka/Ks的值。

2 結(jié) 果

2.1 梅花鹿與歐洲馬鹿染色體共線性分析

2.1.1 基因組序列的同源性 為了獲得梅花鹿和歐洲馬鹿染色體序列的共線性，將兩個(gè)物種進(jìn)行基因組序列比對(duì)，結(jié)果顯示，梅花鹿與歐洲馬鹿兩個(gè)基因組序列表現(xiàn)出同源性(圖1)。梅花鹿基因組序列長(zhǎng)度為2 481 763 803 bp，歐洲馬鹿基因組序列長(zhǎng)度為3 438 623 608 bp,梅花鹿與歐洲馬鹿同源基因序列長(zhǎng)度為1 865 425 381 bp，占梅花鹿基因組序列的75.17%，占?xì)W洲馬鹿基因組序列的55.10%。

染色體共線性圖譜顯示了兩種鹿基因組染色體的同源性。紅色線代表兩條染色體之間是正向比對(duì)，藍(lán)色線代表兩條染色體之間是反向比對(duì)。圖中橫坐標(biāo)代表梅花鹿染色體，縱坐標(biāo)代表歐洲馬鹿染色體

本研究以梅花鹿的基因組作為參考基因組，將梅花鹿與歐洲馬鹿染色體序列的比對(duì)結(jié)果進(jìn)行統(tǒng)計(jì)分析，篩選染色體間同源區(qū)相似性為95%及以上，且對(duì)同源區(qū)長(zhǎng)度大于1 kb的數(shù)據(jù)進(jìn)行統(tǒng)計(jì)。以梅花鹿基因組為標(biāo)準(zhǔn)，統(tǒng)計(jì)歐洲馬鹿對(duì)應(yīng)染色體的比對(duì)情況發(fā)現(xiàn)，梅花鹿的33條染色體與歐洲馬鹿對(duì)應(yīng)的染色體均有較好的同源性，且有27條染色體的同源性在95%以上，其中5、10、11、20、21、25、26、27號(hào)染色體比對(duì)率在99%以上，表明梅花鹿與歐洲馬鹿對(duì)應(yīng)染色體有較高的同源性(表2)。

2.1.2 梅花鹿基因組X染色體的確定 本試驗(yàn)在梅花鹿基因組33條染色體中，一直沒有找到與X染色體對(duì)應(yīng)的染色體。結(jié)合兩個(gè)物種的基因組序列對(duì)比結(jié)果發(fā)現(xiàn)，梅花鹿的23號(hào)染色體與歐洲馬鹿的23號(hào)染色體比對(duì)率極低，僅有0.4%(表2)。為了進(jìn)一步驗(yàn)證梅花鹿的23號(hào)染色體是否對(duì)應(yīng)歐洲馬鹿X染色體，將梅花鹿33條染色體與歐洲馬鹿的兩條性染色體(X和Y染色體)分別進(jìn)行比對(duì)的統(tǒng)計(jì)結(jié)果顯示，梅花鹿的23號(hào)染色體與歐洲馬鹿的X染色體比對(duì)同源性最高，為81 051 312 bp(表3)。基于以上兩方面數(shù)據(jù)支持認(rèn)為，梅花鹿23號(hào)染色體對(duì)應(yīng)歐洲馬鹿基因組中的X染色體。

表3 梅花鹿基因組與歐洲馬鹿X、Y染色體序列比對(duì)統(tǒng)計(jì)

2.1.3 基因組的結(jié)構(gòu)變異 從梅花鹿與歐洲馬鹿的基因組序列進(jìn)行比對(duì)得到的染色體共線性圖中發(fā)現(xiàn)了梅花鹿與歐洲馬鹿各條染色體間序列的同源關(guān)系，除了像1號(hào)、21號(hào)染色體有很好的同源性，在圖中呈現(xiàn)平滑連續(xù)的紅色斜線和藍(lán)色斜線外(圖1紅框)；還有像6號(hào)、12號(hào)染色體同源性并不好，在圖中表現(xiàn)出有一段斷開的反向序列，即連續(xù)紅色斜線和連續(xù)藍(lán)線中出現(xiàn)一段斷開的藍(lán)色斜線或者紅色斜線(圖1綠框)。這種染色體間序列中斷開的反向序列是基因組比對(duì)過程中染色體發(fā)生的一種重排現(xiàn)象，被稱為染色體倒位，屬于基因組結(jié)構(gòu)變異的一種[29]。

將梅花鹿與歐洲馬鹿基因組序列的比對(duì)結(jié)果中發(fā)現(xiàn)的基因組結(jié)構(gòu)變異進(jìn)行統(tǒng)計(jì)發(fā)現(xiàn)，梅花鹿的33條 染色體上均有一定數(shù)量的倒位，共計(jì)倒位數(shù)目為37 847個(gè)。4、23、28號(hào)染色體上的倒位數(shù)目較多，超過了4 000個(gè),其中23號(hào)染色體上倒位最多，為5 238個(gè)；而3、10、16、17、19、21、25號(hào)染色體上的倒位數(shù)目較少，低于100個(gè)(表4)。37 847個(gè)倒位片段中，片段長(zhǎng)度在1～5 kb的倒位數(shù)目最多，共有25 281個(gè)；倒位片段長(zhǎng)度>50 kb的大片段倒位共有15個(gè)(圖2)。

圖2 梅花鹿基因組倒位片段大小統(tǒng)計(jì)

表4 梅花鹿各條染色體倒位數(shù)量

2.2 梅花鹿與歐洲馬鹿染色體倒位涉及基因的功能富集分析

通過將梅花鹿與歐洲馬鹿的全基因組序列進(jìn)行比對(duì)，得到兩物種染色體間序列的共線性及對(duì)應(yīng)染色體上基因組的結(jié)構(gòu)變異，這些基因組結(jié)構(gòu)變異主要以染色體上序列片段倒位為主。將因大片段倒位而受到影響的基因大致分為兩類：被倒位斷點(diǎn)截?cái)嗟幕蚝驮诘刮粌?nèi)部的基因。對(duì)這兩類基因分別進(jìn)行了功能富集分析。被倒位截?cái)嗟幕騁O功能富集分析主要涉及的生物學(xué)過程及分子功能主要有：嗅覺中化學(xué)刺激的檢測(cè)(GO:0050911)；嗅覺感覺(GO:0007608)；嗅覺感受器活動(dòng)(GO:0004984)；信號(hào)傳感器活動(dòng)(GO:0004871)；信號(hào)傳導(dǎo)(GO:0007165)。KEGG功能富集分析主要涉及的信號(hào)通路是嗅覺傳導(dǎo)(hsa04740)(圖3，表5)。未被倒位斷點(diǎn)截?cái)嗟幕騁O功能富集分析主要涉及的生物學(xué)過程及分子功能主要有: G蛋白偶聯(lián)受體信號(hào)通路(GO:0007186)；嗅覺中化學(xué)刺激的檢測(cè)(GO:0050911；嗅覺感覺(GO:0007608)；G蛋白偶聯(lián)受體活性(GO:0004930)；嗅覺感受器活動(dòng)(GO:0004984)。KEGG功能富集涉及的主要信號(hào)通路是嗅覺傳導(dǎo)(hsa04740)(圖4，表6)。富集在通路的基因主要是嗅覺受體家族的成員，其中O10 J1除了與嗅覺信號(hào)通路有關(guān)，還可能在受精過程中參與趨化[30]。

表5 被倒位截?cái)嗟幕騁O和KEGG主要的富集統(tǒng)計(jì)

表6 倒位內(nèi)部未被截?cái)嗟幕騁O和KEGG主要的富集統(tǒng)計(jì)

圖4 在倒位內(nèi)部未被截?cái)嗟幕蚬δ芨患?/p>

值得注意的是，被倒位斷點(diǎn)截?cái)嗟幕蚝臀幢坏刮粩帱c(diǎn)截?cái)嗟膬深惢蛑饕墓δ芨患治龌疽恢拢ǎ盒嵊X中化學(xué)刺激的檢測(cè)(GO:0050911)；嗅覺感覺(GO:0007608)；嗅覺感受器活動(dòng)(GO:0004984)；信號(hào)傳感器活動(dòng)(GO:0004871)以及嗅覺傳導(dǎo)(hsa04740)信號(hào)通路。

2.3 梅花鹿與歐洲馬鹿基因共線性的結(jié)果

2.3.1 梅花鹿與歐洲馬鹿直系同源區(qū)鑒別 為了獲得梅花鹿和歐洲馬鹿兩個(gè)基因組間的直系同源區(qū)，將兩個(gè)物種的全部基因序列進(jìn)行共線性分析，在比對(duì)結(jié)果中共檢測(cè)到梅花鹿和歐洲馬鹿基因組間79個(gè)同源區(qū)段，占梅花鹿基因組序列長(zhǎng)度的95.54%。由于組裝的歐洲馬鹿基因組空位共有1.5 Gb，這些同源區(qū)段只占?xì)W洲馬鹿基因組長(zhǎng)度的41.89%；平均每個(gè)同源區(qū)段含有161個(gè)基因?qū)Γ瑔蝹€(gè)同源區(qū)段最多涵蓋777個(gè)基因?qū)Γ采w基因組長(zhǎng)度104 886 123 bp。單個(gè)同源區(qū)段最少涵蓋14個(gè)基因?qū)Γ采w基因組長(zhǎng)度3 494 220 bp。同源區(qū)段平均長(zhǎng)度30 063 070 bp，最大的同源共線性區(qū)域在5號(hào)染色體上，共線性區(qū)段內(nèi)同源基因的平均Ks值0.096(表7)。

表7 梅花鹿與歐洲馬鹿直系同源區(qū)統(tǒng)計(jì)

2.3.2 梅花鹿與歐洲馬鹿進(jìn)化分析 在基因共線性分析中，共檢測(cè)到12 629個(gè)直系同源基因(表7)，為了估計(jì)梅花鹿和歐洲馬鹿兩個(gè)物種的分化時(shí)間，計(jì)算了12 629個(gè)直系同源基因的同義突變率Ks為0.007，每位點(diǎn)每年中性突變率為1.1×10-8，根據(jù)公式：t=Ks/2 μ(其中，t為物種分化時(shí)間，Ks為同義突變率，μ為中性突變率)計(jì)算得出梅花鹿和歐洲馬鹿分化時(shí)間大概是在0.318 MYA。

基于兩個(gè)物種基因共線性結(jié)果分析，計(jì)算出檢測(cè)的同源基因?qū)Φ腒a/Ks值，可以發(fā)現(xiàn)，同源基因?qū)Φ腒a/Ks值主要分布在0～1之間，且最大分布是在0.5～0.75這個(gè)范圍內(nèi)(圖5)。表明梅花鹿與歐洲馬鹿的大部分基因在進(jìn)化過程中受到純化選擇。

圖5 同源基因?qū)a/Ks值的分布

3 討 論

目前，有關(guān)梅花鹿與馬鹿演化關(guān)系的研究報(bào)道僅限于細(xì)胞遺傳學(xué)和系統(tǒng)發(fā)生學(xué)的研究，并沒有比較基因組學(xué)的相關(guān)研究。Wang等[5]通過對(duì)芥屬的黑芥和芥菜的B亞基因組進(jìn)行基因共線性的分析，找到了兩個(gè)物種間染色體上的結(jié)構(gòu)變異，為芥屬植物的起源進(jìn)化提供了理論基礎(chǔ)。而本研究以梅花鹿和歐洲馬鹿兩個(gè)鹿屬動(dòng)物的基因組為基礎(chǔ)，通過對(duì)兩個(gè)物種基因組的同源性以及基因序列的共線性分析，同樣也發(fā)現(xiàn)了兩個(gè)物種基因組的結(jié)構(gòu)變異以及染色體上序列的共線性關(guān)系，證明從基因組層面研究梅花鹿與歐洲馬鹿起源進(jìn)化關(guān)系是可行的，而且類似的研究在鳥類進(jìn)化起源研究中也有涉及[31]。

本研究選取兩個(gè)不同性別物種基因組進(jìn)行分析，就不同性染色體帶來(lái)的基因差異是否會(huì)對(duì)研究結(jié)果產(chǎn)生影響進(jìn)行分析，首先，在研究前期已經(jīng)將歐洲馬鹿的Y染色體與梅花鹿的各條染色體進(jìn)行了比對(duì)，發(fā)現(xiàn)同源區(qū)很少，且同源片段很短。因此，在基因組比對(duì)分析中，歐洲馬鹿Y染色體的存在基本不會(huì)對(duì)常染色體以及X染色體的比對(duì)分析產(chǎn)生影響。其次，對(duì)基因共線性的分析發(fā)現(xiàn)，由于Y染色體上基因數(shù)相較于本研究發(fā)現(xiàn)的12 629個(gè)同源基因來(lái)說(shuō)很少，因此，缺少Y染色體的基因分析對(duì)研究結(jié)果產(chǎn)生的影響基本上是可以忽略不計(jì)的。綜上所述，本研究中因?yàn)椴煌匀旧w帶來(lái)的差異基本不會(huì)對(duì)研究結(jié)果產(chǎn)生影響。

本研究利用基因共線性結(jié)果分析得到直系同源基因的同義突變率Ks，計(jì)算出梅花鹿與歐洲馬鹿的分化時(shí)間為0.318 MYA。這個(gè)結(jié)果與Kuwayama和Ozawa[32]基于細(xì)胞色素b(Cytb)進(jìn)行系統(tǒng)發(fā)育分析得到的歐洲馬鹿與梅花鹿之間的分化時(shí)間為0.80 MYA稍有差距。這可能是因?yàn)榫€粒體基因組代表的是母系遺傳，在研究物種的起源分化時(shí)，物種的全基因組分析會(huì)比單獨(dú)的母系遺傳研究更加系統(tǒng)，更加全面。也可能是兩個(gè)物種基因組測(cè)序或者組裝質(zhì)量的問題，但是基于兩個(gè)物種染色體共線性和基因共線性的同源關(guān)系來(lái)看，二者的分化時(shí)間應(yīng)該是極為接近的。

染色體倒位(chromosomal inversion)是常見的基因組結(jié)構(gòu)變異，常常導(dǎo)致基因排列順序的改變，能夠抑制基因重組，在物種進(jìn)化過程中扮演了重要的作用[33]。Yang等[7]通過對(duì)陸地棉At和Dt兩個(gè)亞基因組非共線性區(qū)域的基因富集到幾種代謝途徑，包括維生素B6代謝和糖胺聚糖降解，這些合成代謝促進(jìn)胚性愈傷組織發(fā)育。本研究對(duì)梅花鹿與歐洲馬鹿染色體上倒位基因進(jìn)行GO和KEGG功能富集到的主要的生物學(xué)過程及分子功能有：嗅覺中化學(xué)刺激的檢測(cè)(GO:0050911)、嗅覺感覺(GO:0007608)、嗅覺感受器活動(dòng)(GO:0004984)、信號(hào)傳感器活動(dòng)(GO:0004871)以及嗅覺傳導(dǎo)(hsa04740)信號(hào)通路。通過功能富集結(jié)果推測(cè)，梅花鹿與歐洲馬鹿在進(jìn)化過程中嗅覺表型可能受到染色體倒位事件的影響，從而發(fā)生了不同的選擇性進(jìn)化。Fan等[34]對(duì)大熊貓、狗和貓的染色體倒位基因進(jìn)行功能富集，同樣發(fā)現(xiàn)富集到的基因功能主要與嗅覺感覺表型相關(guān)，而對(duì)于梅花鹿與歐洲馬鹿嗅覺表型的差異暫時(shí)沒有相關(guān)內(nèi)容的研究報(bào)道，但是因?yàn)椴煌纳L(zhǎng)環(huán)境和采食偏好性[35]，梅花鹿與歐洲馬鹿在嗅覺表型上應(yīng)該會(huì)發(fā)生不同的選擇性進(jìn)化。就像馴鹿為了可以在寒冷的環(huán)境下生存，嗅覺進(jìn)化的非常發(fā)達(dá)，以至于可以發(fā)現(xiàn)埋在厚厚積雪下的食物[36]。關(guān)于導(dǎo)致梅花鹿與歐洲馬鹿嗅覺表型差異的主效基因的研究將會(huì)是我們后續(xù)的研究?jī)?nèi)容。

4 結(jié) 論

本研究基于比較基因組學(xué)的方法，以梅花鹿與歐洲馬鹿兩個(gè)物種染色體水平基因組為研究對(duì)象進(jìn)行分析，初步獲得了兩個(gè)物種基因組的同源關(guān)系以及在進(jìn)化過程中發(fā)生的染色體倒位現(xiàn)象，估算出梅花鹿與歐洲馬鹿的分化時(shí)間為0.318 MYA，并借助功能富集分析將染色體倒位現(xiàn)象與嗅覺表型的差異聯(lián)系起來(lái)，為鹿屬動(dòng)物染色體進(jìn)化研究提供更多的理論基礎(chǔ)。