基于線粒體DNA 4 個基因/區域的矮小梅花鹿群體結構與起源進化分析

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

(1.中國農業科學院 特產研究所，長春 130112; 2.吉林農業大學 中藥材學院，長春 130118)

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

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

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

1 材料與方法

1.1 材 料

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

1.2 方 法

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

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

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

PCR反應體系為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。

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

1.3 數據分析

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

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

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

2 結果與分析

2.1 序列變異分析

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

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

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

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

2.2 群體結構分析

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

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

因為該研究群體可以分為明顯的 3 個亞群，因此以亞群為單元研究群體內結構，并將4 個基因/區域串聯，結果表明該群體總體平均距離為0.078，而H3與H1、H2遺傳距離遠大于H1和H2之間的距離(表5)，與系統發育樹結果一致，因為亞群內個體間無序列變異，遺傳距離均為0。

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

注：Loci列數字表示突變所在基因/區域位置。

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 系統發育分析

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

2.4 群體歷史分析

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

圖1 基于最大似然法構建的系統發育樹Fig.1 Phylogenetic tree based on maximum likelihood method

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

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

圖2 基于 12S、Cytb、 ND5和D-loop基因/區域對矮小梅花鹿群體的錯配分布分析Fig.2 Mismatch distribution analysis of dwarf Sika deer population based on four gene/region

3 討 論

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

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

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

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

3.2 矮小型梅花鹿群體結構

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

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

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

4 結 論

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