不同海拔喜馬拉雅旱獺線粒體基因組比較分析

南新營 李 優(yōu) 李耀東

(青海大學(xué)，西寧，810016)

喜馬拉雅旱獺(Marmotahimalayana)屬哺乳綱(Mammalia)嚙齒目(Rodentia)松鼠科(Sciuridae)旱獺屬，體形肥大，是松鼠科中體型最大的一種，主要棲息于我國青藏高原及其毗鄰的印度、尼泊爾及巴基斯坦等地的各類高寒草甸草原[1]。作為世居高原的一個特有物種，喜馬拉雅旱獺經(jīng)過長期的自然選擇和自身適應(yīng)，在生理生化上獲得了適應(yīng)高原地區(qū)低溫低氧的穩(wěn)定遺傳學(xué)特性，是研究高原低氧適應(yīng)的天然理想動物模型[2]。

早在20世紀(jì)80年代末，陳欽銘等[3]首次用喜馬拉雅旱獺作為動物模型來研究其適應(yīng)高原低氧的機(jī)理。隨后，賈榮莉[4]從心臟、血液和內(nèi)分泌等角度分析，發(fā)現(xiàn)喜馬拉雅旱獺的紅細(xì)胞數(shù)及心肌細(xì)胞的形態(tài)結(jié)構(gòu)、生理功能與生化酶均未受到高原低氧的不良影響，表明對高原低氧有很好的適應(yīng)性。隨著二代測序技術(shù)的迅猛發(fā)展，其高通量、快速、低成本的特點(diǎn)成為越來越多的學(xué)者解決生物學(xué)問題的首選。線粒體基因組中編碼蛋白質(zhì)的所有基因都位于電子轉(zhuǎn)移呼吸鏈上，在氧的利用和能量代謝中起著關(guān)鍵作用。Xu等[5]認(rèn)為線粒體基因的進(jìn)化可能反映了高原適應(yīng)的一些過程；Chao等[6]通過對喜馬拉雅旱獺線粒體基因組測序，已獲得喜馬拉雅旱獺線粒體基因組全序列。線粒體能量代謝在生物進(jìn)化和物種形成過程中起著重要作用，本研究以不同海拔喜馬拉雅旱獺線粒體為供試材料，利用二代測序技術(shù)組裝線粒體基因組，通過線粒體差異基因的比較分析，從mtDNA一級結(jié)構(gòu)上篩選喜馬拉雅旱獺低氧適應(yīng)相關(guān)的基因，為進(jìn)一步在分子水平上闡明高原低氧適應(yīng)機(jī)理提供理論依據(jù)。

1 材料與方法

1.1 樣品采集

在青海省海東市樂都區(qū)引勝鄉(xiāng)(海拔2 279 m，3例)、黃南州同仁縣瓜什則鄉(xiāng)(海拔3 273 m，3例)和玉樹州曲麻萊縣麻多鄉(xiāng)(海拔4 452 m，2例)3個地區(qū)采集喜馬拉雅旱獺肝臟組織共8例，于液氮中保存并帶回，冰箱內(nèi)-80 ℃儲存。

1.2 喜馬拉雅旱獺基因組DNA提取

用痕量樣品基因組DNA提取試劑盒(Qiagen公司)提取基因組DNA。

1.3 喜馬拉雅旱獺線粒體基因組測序

采用Illumina HiSeqX10[7]測序技術(shù)對8個樣品的DNA進(jìn)行paired-end測序，用fastQC軟件(http：//www.bioinformatics.babraham.ac.uk/projects/fastqc/)觀測測序質(zhì)量，然后使用NGSQC[8]軟件對測序數(shù)據(jù)質(zhì)控。測序由百邁客生物科技有限公司完成。

1.4 喜馬拉雅旱獺線粒體基因組組裝

使用SPAdes[9]對質(zhì)控過的clean data初步拼接，經(jīng)過同源比對和蛋白序列比對后，利用PRICE(Paired-Read Iterative Contig Extension)[10]進(jìn)行迭代延伸、gap補(bǔ)洞及重拼接，直至獲得完整的線粒體基因組序列。

1.5 數(shù)據(jù)分析

獲得完整的線粒體基因組序列后，使用MEGA X軟件[11]分析不同海拔喜馬拉雅旱獺線粒體基因組大小和堿基組成。采用mafft軟件[12]進(jìn)行序列比對，用Mega X軟件采用最大似然法(Maximum likelihood，ML)[13]構(gòu)建系統(tǒng)發(fā)育樹。用SPSS 19.0對3個海拔喜馬拉雅旱獺線粒體基因組中13個蛋白編碼基因的突變位點(diǎn)個數(shù)進(jìn)行統(tǒng)計(jì)學(xué)分析。使用MEGA X做基于密碼子的比對。用KaKs_Calculator 2.0[14](https：//sourceforge.net/projects/kakscalculator2/)計(jì)算Ka/Ks，計(jì)算方法選擇γ-YN，若Ka/Ks>1，則認(rèn)為有正選擇效應(yīng)；若Ka/Ks=1，則認(rèn)為存在中性選擇作用；若Ka/Ks<1，則認(rèn)為是純化選擇作用。密碼子偏好性分析使用 CodonW 軟件，同義密碼子相對使用度(relative synonymous codon usage，RSCU)等于1，則密碼子使用沒有偏好；若大于1則表明其使用頻率相對較高，反之則該密碼子使用頻率相對較低。

2 結(jié)果與分析

2.1 喜馬拉雅旱獺線粒體基因組測序

2.1.1 測序數(shù)據(jù)統(tǒng)計(jì)

通過Illumina HiSeq X10高通量測序平臺，采用2×150 bp雙末端策略，對8個樣品的文庫進(jìn)行線粒體基因組測序。根據(jù)原始數(shù)據(jù)過濾后結(jié)果，Q20值均在97%以上，Q30值均在93%以上(表1)。表明測序數(shù)據(jù)整體質(zhì)量較好，可用于后續(xù)線粒體基因組的組裝。

表1 喜馬拉雅旱獺樣本過濾后的測序數(shù)據(jù)

2.1.2 線粒體基因組序列

測序獲得了3個海拔共8個樣品的喜馬拉雅旱獺線粒體基因組(表2)，發(fā)現(xiàn)不同海拔喜馬拉雅旱獺線粒體序列長度之間并無差異，均為16 442 bp，其核苷酸組成也無差異，A、T、G、C平均含量為32.1%、31.4%、12.8%、23.7%，有明顯的AT偏倚。每個樣品由37個片段拼接而成，最大片段為1 827 bp，最小片段為59 bp，GC含量為36.5%，8個樣品之間沒有差異。

表2 喜馬拉雅旱獺線粒體DNA序列信息

2.2 系統(tǒng)發(fā)育分析

黃婭琳等[15]研究發(fā)現(xiàn)基于單一基因位點(diǎn)構(gòu)建系統(tǒng)發(fā)育關(guān)系容易造成偏差，使用線粒體基因組全序列構(gòu)建的系統(tǒng)發(fā)育樹更為精確可靠。基于線粒體基因組組裝好的序列，以溫哥華旱獺(Marmotavancouverensis)[16](序列號：MK859 897.1)為外群，用MEGA X軟件采用最大似然法構(gòu)建系統(tǒng)發(fā)育樹(圖1)，從系統(tǒng)樹可看出，黃南的2、3號聚為一支，再與樂都的3號聚為一支，樂都的1、2號聚為一支后又同黃南1號聚為一支，最后與玉樹的1、2號聚為一支，表明與玉樹相比，樂都和黃南的喜馬拉雅旱獺親緣關(guān)系相對較近。

圖1 ML法構(gòu)建的不同海拔喜馬拉雅旱獺發(fā)育樹Fig.1 Phylogenetic relationship of eight Marmota himalayana constructed by ML approach

2.3 喜馬拉雅旱獺線粒體DNA蛋白質(zhì)編碼基因差異

2.3.1 13個蛋白編碼基因突變位點(diǎn)個數(shù)

對不同海拔喜馬拉雅旱獺線粒體基因組中13個蛋白編碼基因的突變位點(diǎn)個數(shù)進(jìn)行統(tǒng)計(jì)，發(fā)現(xiàn)共有65個突變位點(diǎn)，其中COX3基因最為保守，沒有發(fā)生變異，而ND5基因有18個突變位點(diǎn)，比例最高(27.7%)。玉樹的突變位點(diǎn)個數(shù)極顯著高于樂都和黃南(P<0.01)，即隨著海拔的增加，突變位點(diǎn)個數(shù)顯著增多(表3)。

表3 不同海拔喜馬拉雅旱獺13個蛋白編碼基因突變位點(diǎn)個數(shù)

2.3.2 功能基因蛋白編碼序列差異

分別比較了13個編碼蛋白的核苷酸和氨基酸序列，發(fā)現(xiàn)除了COX3基因沒有突變，COX1、COX2、ND4L、ATP6和ATP8雖然在基因水平上發(fā)生了突變，但并沒有發(fā)生相應(yīng)的氨基酸改變，即同義突變(same sense mutation)，共有41個基因位點(diǎn)發(fā)生了同義突變，如玉樹-1、玉樹-2的COX1基因第1 380位G→A，COX2基因第648位T→A，ND4L基因的第513位T→C，ATP6基因第324位A→G(圖2，圖3)，ATP8基因第6位C→T。而ND1、ND2、ND3、ND4、ND5、ND6和Cytb這些基因某些位點(diǎn)的改變導(dǎo)致了氨基酸的改變，即非同義突變(non-synonymous mutation)，這7個基因里共有24個非同義突變(表4)。

表4 蛋白編碼基因非同義突變位點(diǎn)

續(xù)表4

圖2 ATP6的核苷酸序列比較Fig.2 Comparison of nucleotide sequences of ATP6 注：從上到下依次為YS-1、YS-2，LD-1、LD-2、LD-3、HN-1、HN-2、HN-3，下同 Note：From top to bottom are YS-1，YS-2，LD-1，LD-2，LD-3，HN-1，HN-2 and HN-3，the same as below

圖3 ATP6的氨基酸序列比較Fig.3 Comparison of amino acid sequences of ATP6

由圖4和圖5可見，黃南的非同義突變位點(diǎn)個數(shù)最多，有13個，如HN-2、HN-3的ND4基因第335位C→T，導(dǎo)致ND4氨基酸序列的第112位A→V(丙氨酸→纈氨酸)，ND5基因第22位T→C，導(dǎo)致ND5氨基酸序列的第8位S→P(絲氨酸→脯氨酸)；玉樹的次之，有8個，YS-1、YS-2的ND5基因第50位T→C，導(dǎo)致ND5氨基酸序列的第17位I→T(異亮氨酸→蘇氨酸)；樂都的最少，只有3個，LD-3的Cytb基因第514位T→C，導(dǎo)致Cytb氨基酸序列的第172位Y→H(酪氨酸→組氨酸)，ND2基因第617位A→G，導(dǎo)致ND2氨基酸序列的第206位N→S(天冬酰胺→絲氨酸)，ND6基因第361位T→A，導(dǎo)致ND6氨基酸序列的第121位L→M(亮氨酸→蛋氨酸)。ND5基因的非同義突變位數(shù)最多，有7個。

圖4 ND6的核苷酸序列Fig.4 Comparison of nucleotide sequences of ND6

圖5 ND6的氨基酸序列Fig.5 Comparison of amino acid sequences of ND6

2.3.3 13個蛋白編碼基因的Ka/Ks

通過13個蛋白編碼基因，計(jì)算不同海拔喜馬拉雅旱獺堿基的Ka、Ks、Ka/Ks，統(tǒng)計(jì)結(jié)果表明(圖6)，Ka中ND5的值最高，ND3次之，COX1、COX2、COX3、ND4L、ATP6、ATP8的Ka值為0；Ks中ATP8最高，COX2次之，COX3最小為0，ND5的Ka/Ks>1，表明ND5基因存在正選擇效應(yīng)，可作為適應(yīng)高海拔低氧環(huán)境的候選基因，其余基因均小于1，表明這些基因存在純化選擇作用。

圖6 喜馬拉雅旱獺蛋白編碼基因Ka/Ks值Fig.6 Ka/Ks values of protein-coding genes of Marmota himalayana

2.3.4 喜馬拉雅旱獺線粒體密碼子偏好性

喜馬拉雅旱獺13個蛋白質(zhì)編碼基因的核苷酸序列共編碼氨基酸3 812個，去掉起始密碼子和終止密碼子后共編碼氨基酸3 802個，在20種氨基酸中，Leu的使用率最高，達(dá)16.18%，Gys的使用率最低，僅有0.71%。在喜馬拉雅旱獺中RSCU值>1的密碼子有32種，RSCU值<1的密碼子有32種，其中AGG的RSCU為0，即在編碼過程中沒有使用該密碼子(圖7)。

圖7 喜馬拉雅旱獺線粒體密碼子偏好Fig.7 Mitochondrial codon bias analysis of Marmota himalayana

3 討論

線粒體作為動物的“能量工廠”，通過氧化磷酸化產(chǎn)生能量，在生物進(jìn)化和物種形成過程中起著重要作用，因此成為研究生物適應(yīng)高原低氧環(huán)境的分子基礎(chǔ)[17]。本研究通過對不同海拔8個樣品的喜馬拉雅旱獺線粒體基因組測序，發(fā)現(xiàn)不同海拔地區(qū)的喜馬拉雅旱獺在線粒體基因組長度和核苷酸組成上并無差異，可能在長期的進(jìn)化過程中，喜馬拉雅旱獺已適應(yīng)了高原地區(qū)低溫低氧氣候，獲得了穩(wěn)定的遺傳特性。對其13個蛋白編碼基因進(jìn)行統(tǒng)計(jì)分析，發(fā)現(xiàn)隨著海拔的增高，突變位點(diǎn)個數(shù)也隨之增加，其中ND5基因的突變位點(diǎn)個數(shù)和非同義突變位點(diǎn)個數(shù)都為最多，且ND5基因Ka/Ks>1，認(rèn)為有正選擇效應(yīng)，表明ND5基因在適應(yīng)低氧中發(fā)揮著重要作用。跟玉樹相比，樂都和黃南的喜馬拉雅旱獺親緣關(guān)系相對較近，且每個海拔梯度之間的親緣關(guān)系較近，表明地理隔離是影響喜馬拉雅旱獺種群交流的重要因素，這與馬英等[18]的研究結(jié)果一致。本研究篩選出ND5這一低氧適應(yīng)相關(guān)基因，為進(jìn)一步研究喜馬拉雅旱獺高原低氧適應(yīng)機(jī)制提供了新的方向和思路。