鄭海錢，徐志旺，徐康寧，李卉，李樹然，張永普

（溫州大學(xué)生命與環(huán)境科學(xué)學(xué)院,浙江 溫州 325035）

受人類活動的影響，地球環(huán)境正面臨前所未有的巨變，工業(yè)革命后劇烈的氣候變化對全球生態(tài)系統(tǒng)產(chǎn)生了巨大的負(fù)面影響（Hoegh-Guldberg.，2019）。由于獨(dú)特的生活史特征，兩棲動物不僅是水生生態(tài)系統(tǒng)中重要的生物類群，還是聯(lián)系水生和陸生環(huán)境的代表（徐士霞等，2004）。受環(huán)境變化的影響，全球超過2 000種兩棲動物的生存正面臨威脅，其中一半可能極度瀕?；?yàn)l臨滅絕（González-Del-Pliego.，2019）。群體滅絕率的升高與遺傳多樣性的喪失密切相關(guān)（Saccheri.，1998），低遺傳多樣性的群體通常表現(xiàn)出低適合度和高滅絕率（Markert.，2010）。研究兩棲動物的遺傳多樣性有助于評估現(xiàn)有物種或群體的生存狀態(tài)，從而為保護(hù)策略的制定提供科學(xué)的理論指導(dǎo)（崔朝霞等，2011）。

mtDNA以結(jié)構(gòu)簡單、進(jìn)化速率較快、嚴(yán)格遵循母系遺傳，以及幾乎不發(fā)生重組的特點(diǎn)，被廣泛應(yīng)用于系統(tǒng)進(jìn)化研究和生物遺傳多樣性評價(jià)（Avise.，1987；宋潔等，2008；Chen.，2009；楊文嘉，2015；杜民等，2021）。其中，部分mtDNA單個(gè)基因的分析已經(jīng)得到普遍應(yīng)用，但是不同基因進(jìn)化速率往往存在差異，c是唯一由mtDNA編碼的細(xì)胞色素，參與線粒體呼吸鏈電子傳遞，進(jìn)化速率適中；基因是mtDNA中比較保守的部位，進(jìn)化速率略低于基因；而線粒體基因和基因最保守（Hwang.，1999）。這導(dǎo)致同一物種采用不同基因分析可能會得出不同的結(jié)論，因此采用多基因聯(lián)合分析物種系統(tǒng)發(fā)育與遺傳多樣性，以期得到更準(zhǔn)確的結(jié)果（董依萌等，2020）。

青藏高原以其特殊的地理位置和氣候條件（如溫度、氧濃度、紫外線輻射等），孕育了一大批獨(dú)特的生物資源，使青藏高原成為中國重要的生態(tài)安全屏障，以及研究生物對高原環(huán)境適應(yīng)性進(jìn)化的熱點(diǎn)地區(qū)（魯春霞等，2004；高偉等，2019）。高原林蛙是青藏高原常見的一種兩棲動物，分布海拔為2 000～4 400 m，背面皮膚較粗糙，呈灰褐色、棕褐色、棕紅色或灰棕色。對高原林蛙的研究集中于活動特征（齊銀等，2007a，2007b）、種間競爭（張晉東等，2007；陳偉等，2018）、生活史（Zhao.，2014；Yu.，2018）和種間遺傳關(guān)系（Yang.，2017；Wang.，2020）等方面，遺傳多樣性研究相對較少（Zhou.，2013），海拔梯度差異研究更少。本研究基于線粒體、、和基因?qū)Σ煌０危? 000 m、2 600 m、3 200 m和3 800 m）的高原林蛙群體進(jìn)行遺傳多樣性研究，并通過對和2個(gè)蛋白編碼基因的正選擇分析，探討高原林蛙的高原適應(yīng)機(jī)制，以期為其保護(hù)提供科學(xué)依據(jù)。

1 材料與方法

1.1 樣品采集

樣品于2020年7—8月和2021年7月采自青海省民和縣 巴 州 鎮(zhèn)（100°45′E，36°12′N，海 拔2 000 m；=22）、海南藏族自治州共和縣龍羊峽水庫（100°41′E，36°4′N，海拔 2 600 m；=30）、共和縣湖東種羊場（100°46′E，36°42′N，海拔 3 200 m；=24）和果 洛 藏 族 自 治 州 瑪 沁 縣（100°14′E，34°27′N，海拔3 800 m；=24），分別定義為2 000 m群體、2 600 m群體、3 200 m群體和3 800 m群體。剪取高原林蛙左后肢的第三趾（Perry.，2011），95%乙醇固定后-80℃保存。

1.2 DNA提取，PCR擴(kuò)增與測序

采用EZUP柱式動物基因組DNA抽提試劑盒（上海生工生物工程股份有限公司），并參考說明書，提取樣本總DNA，將DNA溶于CE Buffer中，用超微量核酸蛋白測定儀（Nano-600，上海嘉鵬科技有限公司）檢測濃度，合格樣品-80℃保存。

PCR 反應(yīng)體系為 25 μL：2×Taq PCR StarMix（北京康潤誠業(yè)生物科技有限公司）12.5 μL，去離子水9.5 μL，DNA模板1 μL，上、下游引物（表 1）各1 μL。PCR反應(yīng)條件為：95℃預(yù)變性3 min；95℃變性30 s、退火30 s（退火溫度分別為：57.1℃、：50℃、：50℃、：58℃）、72℃延伸30 s，35個(gè)循環(huán)；72℃延伸5 min，通過1%瓊脂糖凝膠電泳（150 V，25 min）檢測PCR產(chǎn)物，送北京擎科生物技術(shù)有限公司雙向測序。

表1 本研究中使用的引物Table 1 The primers used in this study

1.3 數(shù)據(jù)處理

將所得雙向序列經(jīng)ContigExpress拼接后與NCBI數(shù)據(jù)庫進(jìn)行BLAST比對，排除假基因的可能，用DNAMAN（Ji.，2020）進(jìn)行序列對比，剪去多余的片段，保存成fasta格式，使用Sequence Matrix將4個(gè)基因合并成聯(lián)合數(shù)據(jù)集。

利用DNASP 5.0（Librado.，2009）分析多態(tài)位點(diǎn)數(shù)、單倍型多樣性（）、核苷酸多樣性（）、平均核苷酸差異數(shù)、錯(cuò)配分布；利用MEGA7.0（Kumar.，2016）分析堿基組成，計(jì)算變異位點(diǎn)數(shù)、簡約信息位點(diǎn)數(shù)、單突變位點(diǎn)；以倭蛙和桓仁林蛙為外 群（NCBI），用 IQtree構(gòu) 建 最 大 似 然 樹（Nguyen.，2015）。

基于Easycodeml（Gao.，2019）分別對線粒體的2個(gè)蛋白編碼基因（和）進(jìn)行選擇壓力分析，在分支模型上選擇位點(diǎn)模型，位點(diǎn)模型主要選用M0（單一比率）、M1a（近中性）、M2a（正選擇）、M3（離散）、M7（beta）、M8（beta& ω）、M8a（beta& ω=1），利用似然比檢驗(yàn)（LRT）統(tǒng)計(jì)評估成對模型：M0 vs.M3、M1a vs.M2a、M7 vs.M8、M8 vs.M8a，檢測正選擇壓力。使用Arlequin3.5（Excoffier&Lischer，2010）計(jì)算群體遺傳分化系數(shù)（），并利用分子變異分析（analysis of molecular variance，AMOVA）檢測遺傳變異的來源，通過1 000次重復(fù)抽樣檢驗(yàn)不同遺傳結(jié)構(gòu)水平上協(xié)方差的顯著性，進(jìn)行 Tajima’s和 Fu’s中性檢驗(yàn)（Tajima，1989）評估高原林蛙群體動態(tài)，基因流（）=（1-）/4（Slatkin.，1987）。通過 POPART（Leigh.，2015）繪制單倍型TCS網(wǎng)絡(luò)圖，推算各個(gè)單倍型間的進(jìn)化關(guān)系。

2 結(jié)果與分析

2.1 線粒體基因序列特征遺傳多樣性分析

共得到100只高原林蛙的、、和基因序列片段，將每只個(gè)體的4條線粒體基因片段合并為1條長度為4 055 bp的線粒體序列，即4個(gè)線粒體基因聯(lián)合數(shù)據(jù)集。檢測到63個(gè)變異位點(diǎn)，其中單一信息位點(diǎn)和簡約信息位點(diǎn)分別有31個(gè)和32個(gè)，A+T含量為55.4%。

4個(gè)群體線粒體基因聯(lián)合數(shù)據(jù)集定義了34個(gè)單倍型，其中，2 000 m群體、2 600 m群體、3 200 m群體、3 800 m群體分別為10個(gè)、9個(gè)、5個(gè)、11個(gè)，共55個(gè)多態(tài)位點(diǎn)，、和平均核苷酸差異數(shù)分別為0.001 62±0.000 08、0.935±0.012和 6.443。3 200 m群體的遺傳多樣性最低，基于4個(gè)線粒體基因聯(lián)合數(shù)據(jù)集的和分別為0.587±0.102和0.000 19±0.000 05；2 600 m群體的遺傳多樣性最高，基于4個(gè)線粒體基因聯(lián)合數(shù)據(jù)集的和分別為0.860±0.036和0.001 21±0.000 21；2 000 m群體的遺傳多樣性高于3 800 m群體（表2）。

表2 基于線粒體基因分析青海省4個(gè)高原林蛙群體的遺傳多樣性Table 2 Genetic diversity of 4 Rana kukunoris populations in Qinghai Province based on mitochondrial genes

2.2 群體遺傳結(jié)構(gòu)與遺傳分化

TCS單倍型網(wǎng)絡(luò)圖顯示（圖1），高原林蛙2 000 m群體和2 600 m群體組成一個(gè)進(jìn)化單元。4個(gè)基因聯(lián)合建樹結(jié)果表明（圖2），系統(tǒng)發(fā)育樹分為2個(gè)支系，2 000 m群體多數(shù)樣本（15/22）、2 600 m群體多數(shù)樣本（24/30）和3 200 m群體全部樣本為1個(gè)支系，2 000 m群體少數(shù)樣本（7/22）、2 600 m群體少數(shù)樣本（6/30）和3 800 m群體全部樣本為另一支系。

圖1 青海省4個(gè)高原林蛙群體線粒體基因的TCS網(wǎng)絡(luò)Fig.1 TCS network of mitochondrial gene of 4 Rana kukunoris populations in Qinghai Province

圖2 基于線粒體基因構(gòu)建的青海省4個(gè)高原林蛙群體最大似然法系統(tǒng)發(fā)育樹Fig.2 Maximum-likelihood phylogenetic tree of 4 Rana kukunoris populations in Qinghai Province based on mitochondrial genes

群體間的為0.090～0.859，除了2 000 m群體和2 600 m群體外，其他兩兩群體間均發(fā)生了極顯著的遺傳分化（表3）。兩兩群體間的為0.09～2.53。AMOVA結(jié)果顯示（表4），線粒體基因群體間分子變異和群體內(nèi)分子變異分別占63.01%和36.99%，分子變異主要發(fā)生在群體間?？倿?.630，且＜0.001。群體間的是 2.52?；诤突虻倪x擇壓力分析結(jié)果表明（表5），M8模型更靈敏，M8模型中，基因檢測到2個(gè)正選擇位點(diǎn)，包括2 600 m群體檢測到1個(gè)正選擇位點(diǎn)（天冬氨酸D→天冬酰胺N），3 200 m和3 800 m群體檢測到1個(gè)正選擇位點(diǎn)（天冬酰胺N→賴氨酸K），基因沒有檢測到正選擇位點(diǎn)。

表3 基于線粒體基因的青海省4個(gè)高原林蛙群體間遺傳分化系數(shù)（對角線上）和基因流（對角線下）Table 3 Genetic differentiation coefficient（above the diagonal）and gene flow（below the diagonal）among 4 Rana kukunoris populations in Qinghai Province based on mitochondrial genes

表4 基于線粒體基因的青海省高原林蛙4個(gè)群體分子方差分析結(jié)果Table 4 Results of molecular variance analysis of 4 Rana kukunoris populations in Qinghai Province based on mitochondrial genes

表5 青海省高原林蛙4個(gè)群體線粒體基因選擇壓力分析Table 5 Analysis of selection pressure of 4 Rana kukunoris populations in Qinghai Province based on mitochondrial genes

2.3 群體歷史動態(tài)

4個(gè)群體中，除2 600 m群體的Fu’s檢驗(yàn)值為正值（0.877），其余3個(gè)群體均為負(fù)值，3 800 m群體的Tajima’s和Fu’s絕對值最高、2 600 m群體的最低，Tajima’s統(tǒng)計(jì)顯示2 000 m群體和3 800 m群體差異顯著，F(xiàn)u’s統(tǒng)計(jì)顯示3 800 m群體差異顯著（表6）?；诰€粒體基因聯(lián)合數(shù)據(jù)集的高原林蛙整體的Fu’s和Tajima’s均為負(fù)值，但差異都不顯著，就單個(gè)群體而言，大多數(shù)群體的Tajima’s與Fu’s均為負(fù)值，且只有少數(shù)群體出現(xiàn)顯著差異，另外，線粒體基因聯(lián)合數(shù)據(jù)集結(jié)果顯示（圖3），4個(gè)群體的錯(cuò)配分布曲線呈雙峰。

表6 青海省4個(gè)高原林蛙群體線粒體基因序列中性檢測Table 6 Neutral test of 4 Rana kukunoris populations in Qinghai Province based on mitochondrial genes

圖3 基于線粒體基因的青海省4個(gè)高原林蛙群體的錯(cuò)配分布曲線Fig.3 Mismatch distribution curve of 4 Rana kukunoris populations in Qinghai Province based on mitochondrial genes

3 討論

3.1 高原林蛙群體的遺傳多樣性

和是衡量遺傳多樣性的2個(gè)重要參數(shù)（劉紹平等，2010）。Grant和Bowen（1998）提出，和的標(biāo)準(zhǔn)臨界值分別為0.5和0.005，將2個(gè)指數(shù)的搭配分為 4 種類型：＜0.5，＜0.005；＞0.5，＜0.005；＜0.5，＞0.005；＞0.5，＞0.005?；诰€粒體4個(gè)基因聯(lián)合數(shù)據(jù)集分析的遺傳多樣性結(jié)果屬于第2種類型，這可能是高原林蛙群體曾經(jīng)歷過瓶頸效應(yīng)后，伴隨了迅速的群體擴(kuò)張與變異的積累。這與北非綠蛙的遺傳多樣性相似（Farjallah.，2012）。此外，基于線粒體基因?qū)Σ煌乩砣后w的高原林蛙的遺傳多樣性研究表明，青藏高原北部地區(qū)適宜生境多于其他地區(qū)，其遺傳多樣性也高于其他地區(qū)（Zhou.，2013）。由此推測適宜生境的數(shù)量或面積影響遺傳多樣性，本研究中線粒體4個(gè)基因聯(lián)合數(shù)據(jù)集對不同海拔的高原林蛙遺傳多樣性結(jié)果表明，2 000 m群體和2 600 m群體的和高于其他群體，推測低海拔生境適宜度整體上高于高海拔，但整體遺傳多樣性仍處于較低水平。而物種的遺傳多樣性降低將導(dǎo)致群體生存能力降低，甚至滅絕（Pauls.，2013）。目前高原林蛙被世界自然保護(hù)聯(lián)盟（IUCN）列為無危（LC），但其較低的群體遺傳多樣性提示該物種同樣值得進(jìn)一步關(guān)注。

3.2 群體遺傳結(jié)構(gòu)和遺傳分化

系統(tǒng)發(fā)育樹的結(jié)果顯示，2 000 m和2 600 m群體的多數(shù)樣本以及3 200 m群體聚為一支，表明這3個(gè)群體的遺傳關(guān)系較近。基于4個(gè)線粒體基因聯(lián)合數(shù)據(jù)集的所有單倍型中，無4個(gè)群體共享單倍型，且大部分單倍型是群體獨(dú)有，說明在進(jìn)化過程中，每個(gè)群體為了適應(yīng)外界環(huán)境形成了獨(dú)特單倍型。3 800 m群體的單倍型數(shù)量高于其他群體，可能是高海拔地區(qū)人煙稀少、人為干擾較少，自然環(huán)境中高原林蛙的基因交流較頻繁。AMOVA結(jié)果表明，群體間的差異大于群體內(nèi)，遺傳分化主要發(fā)生在群體間。Wright等（1978）建議：為0～0.05表示群體間遺傳分化很小、0.05～0.15表示存在中等程度的遺傳分化、0.15～0.25表示遺傳分化較大、＞0.25表示有很大的遺傳分化。本研究中基于4個(gè)線粒體基因聯(lián)合數(shù)據(jù)集的為0.630，表明高原林蛙群體遺傳分化程度很大。

群體間的基因流與遺傳分化負(fù)相關(guān)，當(dāng)＜1時(shí)，基因流弱，遺傳分化處于較強(qiáng)水平；當(dāng)1＜＜4時(shí)，基因流中等，遺傳分化處于中等水平；當(dāng)＞4時(shí)，基因流強(qiáng)，遺傳分化水平較弱（Slatkin.，1987；姚佩君，2020）。4個(gè)線粒體基因聯(lián)合數(shù)據(jù)集的N為2.52，說明4個(gè)群體間的基因交流中等，但除了2 000 m群體和2 600 m群體間的N為2.53外，其他兩兩群體間的N均＜1，和上述遺傳分化結(jié)果相對應(yīng)，青藏高原高原林蛙不同海拔群體間存在一定的地理隔離，這可能與兩棲動物的遷移能力較弱，對所處生境較敏感有關(guān)（李郊，2014）。

檢測正選擇壓力對于揭示生物的適應(yīng)性進(jìn)化機(jī)制具有重要意義（許瑜婷，高芳鑾，2021）。本研究中，基因上沒有檢測到受正選擇的氨基酸位點(diǎn)，說明基因未受到選擇壓力作用；而除了最低的2 000 m群體，其余3個(gè)群體中基因都檢測到了正選擇氨基酸位點(diǎn)?？紤]到基因在有氧呼吸中的重要作用，其可能在高原林蛙適應(yīng)高海拔低溫、低氧環(huán)境的過程中起到了積極的作用。在高原林蛙的核基因中也發(fā)現(xiàn)了與能量代謝相關(guān)的基因受到了正選擇（Yang.，2012）。

3.3 群體歷史動態(tài)

應(yīng)用Tajima’s與Fu’s中性檢驗(yàn)推測群體歷史時(shí)，如果均為負(fù)值，且達(dá)到顯著差異（＜0.05），則說明序列中含有比中性進(jìn)化模型更多的核苷酸位點(diǎn)變化，預(yù)示著被研究群體可能經(jīng)歷過擴(kuò)張（Tajima.，1989）。4個(gè)線粒體基因聯(lián)合數(shù)據(jù)集的Tajima’s和Fu’s均為負(fù)值，但差異都不顯著。單個(gè)群體中，大多數(shù)群體的Tajima’s與Fu’s也均為負(fù)值，且只有少數(shù)群體出現(xiàn)顯著差異；從錯(cuò)配分布曲線來看，4個(gè)群體總體曲線為雙峰。高原林蛙群體可能未經(jīng)歷過群體擴(kuò)張事件。

2 000 m和2 600 m的高原林蛙群體遺傳多樣性高于3 200 m和3 800 m群體，4個(gè)海拔的群體遺傳分化程度高，基因可能對其適應(yīng)高海拔環(huán)境起到了積極作用，群體可能未經(jīng)歷擴(kuò)張事件。本研究可為理解高原林蛙適應(yīng)高原環(huán)境及其保護(hù)提供參考。