基于線粒體DNA 4 個基因/區(qū)域的矮小梅花鹿群體結(jié)構(gòu)與起源進(jìn)化分析

王洪亮，郭肖蘭，張然然，王 磊，李 婷，邢秀梅

(1.中國農(nóng)業(yè)科學(xué)院 特產(chǎn)研究所，長春 130112; 2.吉林農(nóng)業(yè)大學(xué) 中藥材學(xué)院，長春 130118)

梅花鹿屬哺乳綱、偶蹄目、鹿科、鹿屬。歷史上，梅花鹿在中國具有廣泛的分布，但受捕獵和社會變遷影響，種群分布與數(shù)量大幅減少。在中國分布的梅花鹿可劃分為6 個亞種：東北亞種、山西亞種、四川亞種、華北亞種、華南亞種、臺灣亞種。其中華北、山西和臺灣3 個亞種已經(jīng)被宣布在野外滅絕[1]。目前，梅花鹿已被國際自然與自然資源保護(hù)聯(lián)盟(IUCN)編寫的紅皮書列為瀕危物種，也是中國的Ⅰ級重點保護(hù)動物[2]。

線粒體DNA是動物機體細(xì)胞線粒體內(nèi)共價閉合的環(huán)狀雙鏈DNA，包括一條重鏈和一條輕鏈，能夠進(jìn)行自我復(fù)制、轉(zhuǎn)錄和翻譯，具有嚴(yán)格母系遺傳、分子結(jié)構(gòu)簡單、一級結(jié)構(gòu)進(jìn)化快的特點，因此可作為研究動物起源進(jìn)化、群體遺傳結(jié)構(gòu)的重要分子標(biāo)記。鹿科動物線粒體DNA有13個蛋白編碼基因、22個tRNA基因、2個rRNA基因和一個特殊的D-loop區(qū)構(gòu)成。D-loop區(qū)不編碼基因，不受選擇壓力影響，因此突變速率較高，并且得到積累，而蛋白編碼區(qū)受選擇壓力影響，突變速率慢于D-loop區(qū)。在鹿科動物中，13個蛋白編碼基因包括2個ATP酶亞單位(ATP6和ATP8)，1個細(xì)胞色素b(Cytb)、3個細(xì)胞色素C氧化酶的亞單位(COⅠ、COⅡ和COⅢ)和 7 個NADH還原酶復(fù)合體亞單位(ND1～6、ND4L)[3]。

矮小是自然界常見的一種現(xiàn)象，除了受營養(yǎng)條件影響外，更多地由于遺傳因素導(dǎo)致生長發(fā)育受阻，導(dǎo)致其體高矮于正常個體，目前除在人類發(fā)現(xiàn)外，在雞、馬、豬等動物中也被發(fā)現(xiàn)。矮小型動物盡管個體矮小，但往往具備獨特的生產(chǎn)性能或特點，具有一定的產(chǎn)業(yè)前景。通過對中國梅花鹿資源的調(diào)查、評價和搜集，目前已發(fā)現(xiàn)吉林省通化地區(qū)東豐型梅花鹿存在矮小種群，成年個體比正常個體矮5～10 cm，但目前尚不知導(dǎo)致其矮小的遺傳機制或功能基因。為了研究矮小梅花鹿群體母系起源和群體遺傳結(jié)構(gòu)，本研究對位于線粒體DNA的12S、ND5、Cytb基因和D-loop區(qū)序列擴(kuò)增并測序，通過系統(tǒng)發(fā)育和群體遺傳分析確定梅花鹿矮小型群體的母系起源和群體遺傳結(jié)構(gòu)，期望能夠了解其起源、資源現(xiàn)狀，為矮小型梅花鹿資源的開發(fā)和利用提供理論依據(jù)。

1 材料與方法

1.1 材 料

吉林省通化地區(qū)矮小梅花鹿群體，共31只，其中雄性17只，個體編號為奇數(shù)，雌性14只，個體編號為偶數(shù)，鹿只麻醉后，頸靜脈采血10 mL凍存，用于提取DNA。

1.2 方 法

采用全血基因組DNA提取試劑盒提取總DNA，并通過微量紫外分光光度計檢測質(zhì)量濃度，小于50 ng/μL重新提取，并用蒸餾水調(diào)整至50 ng/μL。

根據(jù)GenBank數(shù)據(jù)庫中梅花鹿線粒體DNA全序列設(shè)計12S、ND5、Cytb基因和D-loop區(qū)引物，將引物設(shè)計在上下游基因區(qū)域內(nèi)，確保獲得完整基因片段，其中ND5基因片段較長，因此分為 2 段擴(kuò)增。引物信息見表1。

表1 引物名稱及序列組成Table 1 The sequence of primers used in this research

PCR反應(yīng)體系為50 μL：基因組DNA 2 μL，上下引物游各1 μL，dNTPs Mixture 4 μL，10×PCR buffer(Mg2+)5 μL，TaqDNA聚合酶(5 U/μL)0.25 μL，用ddH2O補至50 μL。

擴(kuò)增條件：95 ℃預(yù)變性5 min；94 ℃變性30 s，最適溫度退火30 s(表 1)，72 ℃延伸1 min，共36個循環(huán)；最后72 ℃延伸10 min，10 ℃保存。PCR產(chǎn)物用15 g/L瓊脂糖凝膠電泳分析。檢測合格后將PCR 產(chǎn)物送上海生工生物工程公司測序。

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

1.3.1 序列分析 通過DNASTAR 軟件包和Clustal X[4]對序列進(jìn)行拼接、比對、校對和編輯，利用軟件DnaSP 5.0[5]計算各基因序列多態(tài)位點數(shù)、核苷酸多樣性(Nucleotide diversity，Pi) 、單倍型數(shù)及單倍型多樣性( Haplotype diversity，Hd) 。

1.3.2 群體結(jié)構(gòu)分析 采用Modeltest 3.7軟件[6]對12S、Cytb、ND5和D-loop基因/區(qū)域序列進(jìn)行最佳替代模型的篩選，同時計算相關(guān)參數(shù)。利用所得模型、參數(shù)結(jié)合相近梅花鹿亞種對應(yīng)基因序列(各亞種名稱及線粒體DNA基因組登錄號：東北梅花鹿1，Cervusnipponhortulorum-1，NC_013834；東北梅花鹿2，Cervusnipponhortulorum-2，HQ191428；北海道亞種，Cervusnipponyesoensis，AB210267；四川亞種，Cervusnipponsichuanicus，JN389443)，以馴鹿對應(yīng)基因序列為外群(Rangifertarandus，NC_007703)，利用Paup 4.0軟件構(gòu)建單倍型系統(tǒng)發(fā)育樹，探討各單倍型系統(tǒng)發(fā)育關(guān)系。利用Mega 7軟件[7]基于Kimura雙參數(shù)模型計算矮小群體以及各亞群之間的遺傳距離。

1.3.3 群體歷史動態(tài)分析 利用Arlequin 3.52軟件[8]，以中性檢驗和核苷酸不配對分布兩種方法來檢測矮小梅花鹿群體歷史動態(tài)。先以Tajima’s D[9]和Fu’s Fs[10]進(jìn)行中性檢驗，再進(jìn)行核苷酸不配對分布檢測。用最小方差法檢驗核苷酸不配對觀測值與群體擴(kuò)張模型預(yù)期分布之間是否一致。結(jié)合DnaSP 5.0 構(gòu)建錯配分布圖，通過可視化曲線圖分析種群歷史動態(tài)。

2 結(jié)果與分析

2.1 序列變異分析

通過在12S、ND5、Cytb和D-loop 4個基因/區(qū)域相鄰區(qū)域設(shè)計引物、擴(kuò)增，獲得完整的基因序列。矮小型梅花鹿4個基因/區(qū)域序列中A+T含量均明顯高于C+G，符合線粒體DNA序列A+T偏向性，同時，堿基突變類型也存在偏倚，其轉(zhuǎn)換數(shù)遠(yuǎn)大于顛換數(shù)，在該群體中，4 個基因共發(fā)生 80 處轉(zhuǎn)換突變，而顛換突變只發(fā)生 2 處。12S、ND5、Cytb3 個基因沒有檢測到堿基插入或缺失突變(表2)，而D-loop區(qū)存在單堿基和雙堿基插入或缺失，導(dǎo)致其序列長度產(chǎn)生變異，出現(xiàn)988和993 bp 2種(表3)。

表2 基于4個基因/區(qū)域的序列變異分析Table 2 Variation analysis of sequences of four genes/regions

在12S、ND5、Cytb和D-loop 4 個基因/區(qū)域中，以D-loop區(qū)檢測到的突變最豐富，且集中于序列前部，除包含31處轉(zhuǎn)換/顛換突變外，還發(fā)現(xiàn)6處插入/缺失突變，符合線粒體DNA D-loop區(qū)低選擇壓力與突變積累特點，其中有5處為單堿基插入/缺失，1處為雙堿基插入/缺失。而12S基因僅發(fā)現(xiàn)4處突變，核苷酸多樣性在 4 個基因/區(qū)域中最低。在ND5基因 21 處突變中，3 處發(fā)生在密碼子第 1 位，1處發(fā)生在密碼子第 2 位，17 處發(fā)生在密碼子第 3 位。在Cytb基因 26 處突變中，6 處發(fā)生在密碼子第 1 位，1 處發(fā)生在密碼子第 2 位，19 處發(fā)生在密碼子第 3 位。

表3 基于 4 個基因/區(qū)域的遺傳多樣性指數(shù)分析Table 3 Analysis of genetic diversity index of four genes/regions

2.2 群體結(jié)構(gòu)分析

在矮小梅花鹿群體中，12S、ND5、Cytb和D-loop 經(jīng)DnasP 5.0分析，均包含3個單倍型(表3)，各基因單倍型多樣性相同，不同基因之間單倍型在群體中沒有交叉，通過 4 個基因/區(qū)域的單倍型都可以將該群體分為3個個體組成相同的亞群。

對4個基因/區(qū)域串聯(lián)進(jìn)行組合單倍型分析，表明整個群體可清晰劃分為 3 個單倍型H1、H2和 H3(表4)，說明整個矮小群體有 3 個母系起源。

因為該研究群體可以分為明顯的 3 個亞群，因此以亞群為單元研究群體內(nèi)結(jié)構(gòu)，并將4 個基因/區(qū)域串聯(lián)，結(jié)果表明該群體總體平均距離為0.078，而H3與H1、H2遺傳距離遠(yuǎn)大于H1和H2之間的距離(表5)，與系統(tǒng)發(fā)育樹結(jié)果一致，因為亞群內(nèi)個體間無序列變異，遺傳距離均為0。

表4 12S、Cytb、 ND5和D-loop區(qū)突變位置及組合單倍型堿基組成Table 4 Location of mutation in four genes/ regions and base composition of combined haplotype

注：Loci列數(shù)字表示突變所在基因/區(qū)域位置。

Note：The No.in Loci row stand for the location of mutation in gene or region.

表5 矮小梅花鹿群體各亞群間遺傳距離Table 5 Genetic distance between subpopulation based on combined haplotypes

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

以馴鹿基因為外群，通過對4個基因/區(qū)域進(jìn)行聚類分析，結(jié)果表明，4 個基因/區(qū)域得出的聚類關(guān)系一致，所有個體可以明顯地聚為3支(結(jié)果未顯示)，利用 4 個基因/區(qū)域組合單倍型，結(jié)合來自于NCBI其他序列構(gòu)建的系統(tǒng)發(fā)育樹，可以發(fā)現(xiàn)各分支都具有較高的置信度，絕大多數(shù)個體與東北梅花鹿具有更近的親緣關(guān)系，而其中 4 個個體與梅花鹿四川亞種聚為一類(圖1)，表明相比其他個體，這 4 個個體與四川亞種親緣關(guān)系更近，而日本北海道亞種與所有個體親緣關(guān)系較遠(yuǎn)。

2.4 群體歷史分析

通過12S、ND5、Cytb和D-loop區(qū)對矮小型梅花鹿群體進(jìn)行中性檢驗與錯配分析。中性檢驗可利用Tajima’s D 值和 Fu’s FS 值來分析種群是否經(jīng)歷過擴(kuò)張，當(dāng)Tajima’s D 值和 Fu’s FS 值接近于零時表明種群較為穩(wěn)定，為負(fù)值(P<0.05)表明種群近期經(jīng)歷擴(kuò)張，正值(P< 0.05)說明種群可能存在分化。錯配分析可通過分枝末端1～37數(shù)字為個體編號，奇數(shù)代表雄性，偶數(shù)代表雌性 The No. from 1 to 37 at branch end stands for individual in sika deer population, odd number for male and even number for femaleSSD(Sum of square deviations)和r(Harpending’s raggedness index)來計算，當(dāng)這兩個參數(shù)的統(tǒng)計檢驗不顯著(P>0.05)時，表明不能拒絕群體擴(kuò)張的假說，即符合原來群體擴(kuò)張的假說。當(dāng)通過可視化曲線圖對種群的歷史動態(tài)進(jìn)行分析時，錯配分布曲線為單峰表明群體經(jīng)歷過擴(kuò)張，為多峰說明近期未經(jīng)歷過擴(kuò)張。

圖1 基于最大似然法構(gòu)建的系統(tǒng)發(fā)育樹Fig.1 Phylogenetic tree based on maximum likelihood method

經(jīng)過中性檢驗，除Cytb的Tajima’s D<0，其他基因或區(qū)域Tajima’s D 值和 Fu’s FS 值均>0(P>0.05)，且Tajima’s D 值接近于0，說明種群比較穩(wěn)定，并且未經(jīng)歷過擴(kuò)張。而錯配分析顯示，除ND5、Cytb和D-loop區(qū) 的r值統(tǒng)計顯著外，其他SSD和r值統(tǒng)計均顯著(表6)，各基因錯配分布曲線均為多峰(圖2)，綜合判斷拒絕群體擴(kuò)張假設(shè)，認(rèn)為群體近期未經(jīng)歷過擴(kuò)張。

表6 基于 12S、Cytb、 ND5和D-loop基因/區(qū)域的中性檢驗與錯配分析Table 6 Neutrality test and mismatch analysis based on four genes/regions

圖2 基于 12S、Cytb、 ND5和D-loop基因/區(qū)域?qū)Π∶坊谷后w的錯配分布分析Fig.2 Mismatch distribution analysis of dwarf Sika deer population based on four gene/region

3 討 論

3.1 序列特征與群體遺傳多樣性

通過線粒體DNA12S、Cytb、ND5和D-loop 4個基因或區(qū)域的序列變異分析表明，各基因或區(qū)域均存在變異位點，且變異中轉(zhuǎn)換與顛換發(fā)生比例差異較大，與前人研究結(jié)果存在差異，孫平芳[11]對馬鹿與白唇鹿共5個群體線粒體DNA D-loop區(qū)研究發(fā)現(xiàn)49個轉(zhuǎn)換和31個顛換變異；鄧鑄疆[12]對西北馬鹿線粒體DNA D-loop區(qū)研究發(fā)現(xiàn)其變異堿基轉(zhuǎn)換/顛換比例R=2.349，而本研究中，梅花鹿群體線粒體基因變異存在較大偏向性，轉(zhuǎn)換突變數(shù)遠(yuǎn)遠(yuǎn)高于顛換，但是目前其機制和生物學(xué)意義無法確定。

本研究對不同基因核苷酸多樣性比較發(fā)現(xiàn)，其在基因間差異較大，以12S基因最低，D-loop區(qū)最高。D-loop區(qū)不編碼基因和氨基酸，因此受到的選擇壓力最小，各種突變，包括插入、缺失都被記錄并保存下來，而Cytb基因編碼氨基酸，受到的選擇壓大于D-loop區(qū)，插入、缺失突變會導(dǎo)致移碼突變，導(dǎo)致蛋白質(zhì)翻譯中斷，影響功能，這類突變危害較大，會被淘汰并清除，因此很難發(fā)現(xiàn)Cytb存在上述突變；在進(jìn)化速率方面，Cytb進(jìn)化速率慢于D-loop區(qū)，但快于12S、COI等基因，進(jìn)化速率適中。在本研究中，D-loop區(qū)存在較多插入、缺失，988/993 bp中共發(fā)現(xiàn) 31 個轉(zhuǎn)換、顛換突變，平均核苷酸差異數(shù)為 9.075，而Cytb在 1 140 bp 中共發(fā)現(xiàn) 26 個轉(zhuǎn)換、顛換突變，平均核苷酸差異數(shù)為6.318，同ND5接近。說明D-loop區(qū)和Cytb在多態(tài)性水平差異較大。

本研究共發(fā)現(xiàn) 2 種序列長度的D-loop區(qū)，分別為988和993 bp，孫平芳[11]在研究中發(fā)現(xiàn)，白唇鹿、天山馬鹿、甘肅馬鹿、塔河馬鹿、青海馬鹿的D-loop區(qū)序列長度為916～1 071 bp；涂劍鋒等[13]對25 種鹿科動物D-loop區(qū)序列分析，發(fā)現(xiàn)長度為914～1 072 bp ，其中梅花鹿臺灣亞種(EF058308)986 bp，越南梅亞種(AF291881)995 bp，參考目前NCBI數(shù)據(jù)庫中梅花鹿線粒體DNA序列，北海道亞種(AB210267)1 107 bp，屋久島亞種(AB218689)996 bp，四川亞種(JN389443)989 bp，由此可見，較高的插入/缺失多樣性導(dǎo)致各梅花鹿亞種之間D-loop區(qū)序列長度差異。

3.2 矮小型梅花鹿群體結(jié)構(gòu)

本研究在矮小型梅花鹿群體中，12S、Cytb、ND5和D-loop 4 個基因或區(qū)域都發(fā)現(xiàn)存在 3 種單倍型，且其包含個體均相同，通過構(gòu)建 4 個基因或區(qū)域組合單倍型并結(jié)合系統(tǒng)發(fā)育分析結(jié)果，表明該群體具有 3 個母系起源。其中 H1、H2 與梅花鹿東北亞種具有非常近的親緣關(guān)系，說明H1和H2單倍型個體起源于東北亞種，相比之下，H3與東北亞種親源關(guān)系較遠(yuǎn)，而與梅花鹿四川亞種親緣關(guān)系很近，說明H3單倍型個體母系來源于四川亞種。此外，因為線粒體DNA為母系遺傳，即母體細(xì)胞的線粒體DNA只能通過女兒代代相傳，當(dāng)其所產(chǎn)個體為雄性時，將會因為無法傳給后代而丟失，在本研究群體中，H3單倍型包括4個個體，但均為雄性，說明該H3單倍型將不會遺傳給后代，如果這一H3單倍型僅存在于該群體，則H3單倍型從此消失，因此，在保存種質(zhì)資源時，需要群體達(dá)到一定的數(shù)量規(guī)模與合理的配種方案，這樣才能最大化地保存物種的遺傳多樣性，否則會導(dǎo)致關(guān)鍵遺傳信息丟失，使某些性狀改變，甚至危及物種生存。

另一方面，為什么吉林地區(qū)梅花鹿群體會出現(xiàn)四川亞種母系起源，其原因可能分為自然遷徙或人為遷徙，如果是自然遷徙導(dǎo)致，則該事件可能發(fā)生于人類文明出現(xiàn)之前相當(dāng)長的時期；如果是人為遷徙所致，即人類活動導(dǎo)致梅花鹿種群的被動遷移，則該事件發(fā)生與人類文明出現(xiàn)之后，尤其是梅花鹿人工馴養(yǎng)出現(xiàn)以后的引種行為。

在哺乳動物中存在一個有趣的現(xiàn)象，貝格曼規(guī)律(Bergmann’s rule)[14]，即同一物種在高緯度寒冷地區(qū)分布的個體體型會大于低緯度炎熱地區(qū)個體，在梅花鹿中也存在此種現(xiàn)象，有學(xué)者對日本 6 個梅花鹿亞種體型研究發(fā)現(xiàn)，隨著緯度由高到低，雄性梅花鹿體型變小，南部島嶼雄性梅花鹿體質(zhì)量 40 kg左右，而北部群體其雄性體質(zhì)量達(dá)100 kg[15-17]，通常，梅花鹿東北亞種的特點是體型高大，而四川亞種則體型矮小，因此有理由懷疑該矮小梅花鹿群體矮小的原因是該群體內(nèi)含有梅花鹿四川亞種血緣，該假設(shè)需要大量研究來驗證。

4 結(jié) 論

基于線粒體DNA12S、Cytb、ND5和D-loop 4 個基因或區(qū)域，對矮小型梅花鹿群體進(jìn)行群體結(jié)構(gòu)、系統(tǒng)發(fā)育及種群歷史等方面分析，得出以下結(jié)論，序列變異存在較大轉(zhuǎn)換/顛換偏向性；各基因在群體中均包括 3 種單倍型，且每種單倍型個體組成相同，因此認(rèn)為群體有 3 個母系起源，結(jié)合系統(tǒng)發(fā)育樹判斷，H1和H2單倍型個體母系起源為梅花鹿東北亞種，而H3單倍型個體母系起源可能為梅花鹿四川亞種，但該群體中僅存 4 個個體，且均為雄性，因此H3單倍型無法遺傳給后代，最終將會在該群體中消失；種群歷史分析表明，該群體沒有經(jīng)歷過擴(kuò)張。