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

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

(青海大學，西寧，810016)

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

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

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

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

2 結(jié)果與分析

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

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

通過Illumina HiSeq X10高通量測序平臺，采用2×150 bp雙末端策略，對8個樣品的文庫進行線粒體基因組測序。根據(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)基于單一基因位點構(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個蛋白編碼基因突變位點個數(shù)

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

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

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

分別比較了13個編碼蛋白的核苷酸和氨基酸序列，發(fā)現(xiàn)除了COX3基因沒有突變，COX1、COX2、ND4L、ATP6和ATP8雖然在基因水平上發(fā)生了突變，但并沒有發(fā)生相應的氨基酸改變，即同義突變(same sense mutation)，共有41個基因位點發(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這些基因某些位點的改變導致了氨基酸的改變，即非同義突變(non-synonymous mutation)，這7個基因里共有24個非同義突變(表4)。

表4 蛋白編碼基因非同義突變位點

續(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可見，黃南的非同義突變位點個數(shù)最多，有13個，如HN-2、HN-3的ND4基因第335位C→T，導致ND4氨基酸序列的第112位A→V(丙氨酸→纈氨酸)，ND5基因第22位T→C，導致ND5氨基酸序列的第8位S→P(絲氨酸→脯氨酸)；玉樹的次之，有8個，YS-1、YS-2的ND5基因第50位T→C，導致ND5氨基酸序列的第17位I→T(異亮氨酸→蘇氨酸)；樂都的最少，只有3個，LD-3的Cytb基因第514位T→C，導致Cytb氨基酸序列的第172位Y→H(酪氨酸→組氨酸)，ND2基因第617位A→G，導致ND2氨基酸序列的第206位N→S(天冬酰胺→絲氨酸)，ND6基因第361位T→A，導致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個蛋白編碼基因，計算不同海拔喜馬拉雅旱獺堿基的Ka、Ks、Ka/Ks，統(tǒng)計結(jié)果表明(圖6)，Ka中ND5的值最高，ND3次之，COX1、COX2、COX3、ND4L、ATP6、ATP8的Ka值為0；Ks中ATP8最高，COX2次之，COX3最小為0，ND5的Ka/Ks>1，表明ND5基因存在正選擇效應，可作為適應高海拔低氧環(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的使用率最高，達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)生能量，在生物進化和物種形成過程中起著重要作用，因此成為研究生物適應高原低氧環(huán)境的分子基礎[17]。本研究通過對不同海拔8個樣品的喜馬拉雅旱獺線粒體基因組測序，發(fā)現(xiàn)不同海拔地區(qū)的喜馬拉雅旱獺在線粒體基因組長度和核苷酸組成上并無差異，可能在長期的進化過程中，喜馬拉雅旱獺已適應了高原地區(qū)低溫低氧氣候，獲得了穩(wěn)定的遺傳特性。對其13個蛋白編碼基因進行統(tǒng)計分析，發(fā)現(xiàn)隨著海拔的增高，突變位點個數(shù)也隨之增加，其中ND5基因的突變位點個數(shù)和非同義突變位點個數(shù)都為最多，且ND5基因Ka/Ks>1，認為有正選擇效應，表明ND5基因在適應低氧中發(fā)揮著重要作用。跟玉樹相比，樂都和黃南的喜馬拉雅旱獺親緣關(guān)系相對較近，且每個海拔梯度之間的親緣關(guān)系較近，表明地理隔離是影響喜馬拉雅旱獺種群交流的重要因素，這與馬英等[18]的研究結(jié)果一致。本研究篩選出ND5這一低氧適應相關(guān)基因，為進一步研究喜馬拉雅旱獺高原低氧適應機制提供了新的方向和思路。