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

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

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

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

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

矮小是自然界常見的一種現(xiàn)象，除了受營養(yǎng)條件影響外，更多地由于遺傳因素導(dǎo)致生長發(fā)育受阻，導(dǎo)致其體高矮于正常個(gè)體，目前除在人類發(fā)現(xiàn)外，在雞、馬、豬等動(dòng)物中也被發(fā)現(xiàn)。矮小型動(dòng)物盡管個(gè)體矮小，但往往具備獨(dú)特的生產(chǎn)性能或特點(diǎn)，具有一定的產(chǎn)業(yè)前景。通過對中國梅花鹿資源的調(diào)查、評價(jià)和搜集，目前已發(fā)現(xiàn)吉林省通化地區(qū)東豐型梅花鹿存在矮小種群，成年個(gè)體比正常個(gè)體矮5～10 cm，但目前尚不知導(dǎo)致其矮小的遺傳機(jī)制或功能基因。為了研究矮小梅花鹿群體母系起源和群體遺傳結(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只，個(gè)體編號為奇數(shù)，雌性14只，個(gè)體編號為偶數(shù)，鹿只麻醉后，頸靜脈采血10 mL凍存，用于提取DNA。

1.2 方 法

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

根據(jù)GenBank數(shù)據(jù)庫中梅花鹿線粒體DNA全序列設(shè)計(jì)12S、ND5、Cytb基因和D-loop區(qū)引物，將引物設(shè)計(jì)在上下游基因區(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補(bǔ)至50 μL。

擴(kuò)增條件：95 ℃預(yù)變性5 min；94 ℃變性30 s，最適溫度退火30 s(表 1)，72 ℃延伸1 min，共36個(gè)循環(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]計(jì)算各基因序列多態(tài)位點(diǎn)數(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)行最佳替代模型的篩選，同時(shí)計(jì)算相關(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ù)模型計(jì)算矮小群體以及各亞群之間的遺傳距離。

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

2 結(jié)果與分析

2.1 序列變異分析

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

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

在12S、ND5、Cytb和D-loop 4 個(gè)基因/區(qū)域中，以D-loop區(qū)檢測到的突變最豐富，且集中于序列前部，除包含31處轉(zhuǎn)換/顛換突變外，還發(fā)現(xiàn)6處插入/缺失突變，符合線粒體DNA D-loop區(qū)低選擇壓力與突變積累特點(diǎn)，其中有5處為單堿基插入/缺失，1處為雙堿基插入/缺失。而12S基因僅發(fā)現(xiàn)4處突變，核苷酸多樣性在 4 個(gè)基因/區(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 個(gè)基因/區(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個(gè)單倍型(表3)，各基因單倍型多樣性相同，不同基因之間單倍型在群體中沒有交叉，通過 4 個(gè)基因/區(qū)域的單倍型都可以將該群體分為3個(gè)個(gè)體組成相同的亞群。

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

因?yàn)樵撗芯咳后w可以分為明顯的 3 個(gè)亞群，因此以亞群為單元研究群體內(nèi)結(jié)構(gòu)，并將4 個(gè)基因/區(qū)域串聯(lián)，結(jié)果表明該群體總體平均距離為0.078，而H3與H1、H2遺傳距離遠(yuǎn)大于H1和H2之間的距離(表5)，與系統(tǒng)發(fā)育樹結(jié)果一致，因?yàn)閬喨簝?nèi)個(gè)體間無序列變異，遺傳距離均為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ā)育分析

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

2.4 群體歷史分析

通過12S、ND5、Cytb和D-loop區(qū)對矮小型梅花鹿群體進(jìn)行中性檢驗(yàn)與錯(cuò)配分析。中性檢驗(yàn)可利用Tajima’s D 值和 Fu’s FS 值來分析種群是否經(jīng)歷過擴(kuò)張，當(dāng)Tajima’s D 值和 Fu’s FS 值接近于零時(shí)表明種群較為穩(wěn)定，為負(fù)值(P<0.05)表明種群近期經(jīng)歷擴(kuò)張，正值(P< 0.05)說明種群可能存在分化。錯(cuò)配分析可通過分枝末端1～37數(shù)字為個(gè)體編號，奇數(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)來計(jì)算，當(dāng)這兩個(gè)參數(shù)的統(tǒng)計(jì)檢驗(yàn)不顯著(P>0.05)時(shí)，表明不能拒絕群體擴(kuò)張的假說，即符合原來群體擴(kuò)張的假說。當(dāng)通過可視化曲線圖對種群的歷史動(dòng)態(tài)進(jìn)行分析時(shí)，錯(cuò)配分布曲線為單峰表明群體經(jīng)歷過擴(kuò)張，為多峰說明近期未經(jīng)歷過擴(kuò)張。

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

經(jīng)過中性檢驗(yàn)，除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ò)張。而錯(cuò)配分析顯示，除ND5、Cytb和D-loop區(qū) 的r值統(tǒng)計(jì)顯著外，其他SSD和r值統(tǒng)計(jì)均顯著(表6)，各基因錯(cuò)配分布曲線均為多峰(圖2)，綜合判斷拒絕群體擴(kuò)張假設(shè)，認(rèn)為群體近期未經(jīng)歷過擴(kuò)張。

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

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

3 討 論

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

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

本研究對不同基因核苷酸多樣性比較發(fā)現(xiàn)，其在基因間差異較大，以12S基因最低，D-loop區(qū)最高。D-loop區(qū)不編碼基因和氨基酸，因此受到的選擇壓力最小，各種突變，包括插入、缺失都被記錄并保存下來，而Cytb基因編碼氨基酸，受到的選擇壓大于D-loop區(qū)，插入、缺失突變會(huì)導(dǎo)致移碼突變，導(dǎo)致蛋白質(zhì)翻譯中斷，影響功能，這類突變危害較大，會(huì)被淘汰并清除，因此很難發(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 個(gè)轉(zhuǎn)換、顛換突變，平均核苷酸差異數(shù)為 9.075，而Cytb在 1 140 bp 中共發(fā)現(xiàn) 26 個(gè)轉(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òng)物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 個(gè)基因或區(qū)域都發(fā)現(xiàn)存在 3 種單倍型，且其包含個(gè)體均相同，通過構(gòu)建 4 個(gè)基因或區(qū)域組合單倍型并結(jié)合系統(tǒng)發(fā)育分析結(jié)果，表明該群體具有 3 個(gè)母系起源。其中 H1、H2 與梅花鹿東北亞種具有非常近的親緣關(guān)系，說明H1和H2單倍型個(gè)體起源于東北亞種，相比之下，H3與東北亞種親源關(guān)系較遠(yuǎn)，而與梅花鹿四川亞種親緣關(guān)系很近，說明H3單倍型個(gè)體母系來源于四川亞種。此外，因?yàn)榫€粒體DNA為母系遺傳，即母體細(xì)胞的線粒體DNA只能通過女兒代代相傳，當(dāng)其所產(chǎn)個(gè)體為雄性時(shí)，將會(huì)因?yàn)闊o法傳給后代而丟失，在本研究群體中，H3單倍型包括4個(gè)個(gè)體，但均為雄性，說明該H3單倍型將不會(huì)遺傳給后代，如果這一H3單倍型僅存在于該群體，則H3單倍型從此消失，因此，在保存種質(zhì)資源時(shí)，需要群體達(dá)到一定的數(shù)量規(guī)模與合理的配種方案，這樣才能最大化地保存物種的遺傳多樣性，否則會(huì)導(dǎo)致關(guān)鍵遺傳信息丟失，使某些性狀改變，甚至危及物種生存。

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

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

4 結(jié) 論

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