線粒體DNA分析青藏高原東方蜜蜂(Apis cerana)種群遺傳結構

劉意秋，周丹銀，龔雪陽，馮毅楠，和紹禹*，趙文正*

(1.云南農(nóng)業(yè)大學 動物科學技術學院，云南 昆明 650201；2.云南農(nóng)業(yè)大學 東方蜜蜂研究所，云南 昆明 650201)

種群遺傳結構分析的主要內(nèi)容為：種群內(nèi)及種群間等位基因/基因型的遺傳差異和頻率分布。分析結果可為探索種質(zhì)資源的演化形成過程，定義重要進化單元[1]，規(guī)劃其保護、發(fā)展和利用提供重要基礎數(shù)據(jù)[2]。東方蜜蜂(ApisceranaFabricius，1793)是亞洲的土著蜜蜂種類，其自然分布地域遍及亞洲大陸溫帶、熱帶、亞熱帶地區(qū)及其眾多島嶼[3]?；趥€體數(shù)量龐大、適應性和抗逆性較強的特點，東方蜜蜂為其生境內(nèi)的生態(tài)系統(tǒng)提供著十分重要的傳粉服務，對于當?shù)厣鷳B(tài)系統(tǒng)多樣性的維持和穩(wěn)定具有不可替代的作用[4]。此外，在一些植被豐富但耕地稀缺的山區(qū)、半山區(qū)，飼養(yǎng)東方蜜蜂還為當?shù)鼐用竦慕?jīng)濟增收做出貢獻。

線粒體DNA(mtDNA)是真核生物的核外遺傳物質(zhì)(之一)，由于具有母系遺傳、不發(fā)生重組、種內(nèi)變異豐富[5]等特點而被廣泛應用于蜜蜂種群遺傳結構分析。利用mtDNAtRNAleu-CO2基因片段對西方蜜蜂開展了大量調(diào)查研究[6～11]，這些研究結果與Ruttner(1988)利用形態(tài)學方法劃分的西方蜜蜂(Apismellifera)種下分類體系基本一致。然而這段mtDNA片段并不十分適宜對東方蜜蜂進行分析，原因在于該片段長度相對較短[12]甚至在某些東方蜜蜂個體中完全缺失[13]。之前一些學者也利用細胞色素C氧化酶亞基1基因(COX1)[14]NADH脫氫酶亞基2基因(ND2)[15]等mtDNA片段對東方蜜蜂開展了遺傳多樣性分析。本試驗中采用的是細胞色素C氧化酶亞基1-亞基2基因片段(CO1-2)。

我國境內(nèi)的東方蜜蜂又稱中華蜜蜂或中蜂。中華蜜蜂的種群遺傳結構及遺傳多樣性狀況前人已有研究報道，然而這些研究所得結論不盡相同。楊冠煌等對中華蜜蜂資源進行了詳細的形態(tài)學和生物學調(diào)查，調(diào)查結果認為我國境內(nèi)分布著5個東方蜜蜂亞種，其中滇、川、藏省區(qū)分布有西藏亞種(Apisceranaskorikovi)、印度亞種(Apisceranaindica)、阿壩亞種(Apisceranaabansis)和中華亞種(Apisceranacerana)，海南亞種(Apisceranahainana)只分布于海南島上[16,17]。但Radloff等的形態(tài)學研究則認為我國境內(nèi)的東方蜜蜂分為“阿壩”、“華中華東”和“華南”3個形態(tài)類群[18]。在后續(xù)的研究中，一些學者利用線粒體DNA[19,20]和微衛(wèi)星DNA[21]分子標記已經(jīng)找到了支持海南亞種(形態(tài)類群)的分子證據(jù)。然而對于青藏高原地區(qū)東方蜜蜂種群遺傳結構及該區(qū)域西藏亞種和阿壩亞種(形態(tài)類群)間遺傳分化的分子數(shù)據(jù)目前鮮見報道。鑒于此，我們采集了青藏高原地區(qū)(包括滇、川、甘、青、藏等省區(qū))的東方蜜蜂樣本并進行線粒體DNA分子標記檢測和數(shù)據(jù)分析，以期調(diào)查該區(qū)域東方蜜蜂的種群遺傳結構，并為研究東方蜜蜂西藏亞種和阿壩亞種(形態(tài)類群)的遺傳分化積累分子數(shù)據(jù)。

1 材料與方法

1.1 樣本采集與DNA提取

東方蜜蜂樣本于2013年5月-2015年6月采集自中國青藏高原及其周邊5個省、區(qū)的9個采樣點，每個點采集7～30群蜜蜂樣本，共237群。同一采樣點的樣本相互間距不超過10 km，樣本詳細信息見圖1和表1。采樣時從巢脾上直接抓取50只左右的工蜂，即刻保存于100%的乙醇中并盡快帶回實驗室貯存于-75 ℃條件下。針對每群蜜蜂樣本隨機選取1只工蜂，采用經(jīng)典的酚/氯仿法[13]提取全基因組DNA。

1.2 PCR擴增與DNA測序

根據(jù)已發(fā)表的東方蜜蜂線粒體基因組全序列(GenBank登錄號：NC_014295)設計了細胞色素C氧化酶亞基1-亞基2基因片段(CO1-2)的特異性引物：CO1-2-F(5’→3’) CCTCGACGATACTCAGATTATCC，CO1-2-R(5’→3’) TCCAAGTGATGGGACAGTTCAT，由上海生工生物技術有限公司合成。PCR擴增體系及擴增條件沿用之前的研究[22]。擴增產(chǎn)物交由上海生工生物技術有限公司進行正反向測序，測序引物即PCR擴增引物。

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

利用Lasergene 7.0軟件包中的SeqMan對正、反向DNA序列進行檢查、校正、修剪并生成一致序列；從NCBI數(shù)據(jù)庫下載了24條其它地區(qū)東方蜜蜂的CO1-2片段作為參考序列(見圖3b)，利用MEGA 5.05軟件[23]對本次試驗測序獲得的CO1-2一致序列和加入24條參考序列后的全部序列分別進行Alignment計算(在默認參數(shù)下)；使用DnaSP 5.1軟件[24]進行單倍型種類、頻率、單倍型多樣度(Hd)、核苷酸多樣度(π)以及種群間遺傳距離(Fst)矩陣的統(tǒng)計計算；以GeographicDistanceMatrixGenerator v1.2.3生成種群間地理距離矩陣；使用Arlequin v3.11軟件[25]進行Tajima’s D中性檢驗和分子方差分析(AMOVA)，并對種群間遺傳距離Fst矩陣與地理距離矩陣之間實施Mantel檢驗；使用Network 5.03軟件(http://www.fluxus-engineering.com)中的Median-Joining算法對鑒別出的單倍型進行網(wǎng)絡關系構建；以蘇拉威西蜂(Apiskoschevnikovi)同源片段為外群并加入?yún)⒖夹蛄泻螅謩e利用鄰接法(Neighbour-Joining)、最大簡約法(Maximum-Parsimony)[26]和貝葉斯法(Bayesian-Algorithm)[27]對鑒別出的單倍型進行系統(tǒng)發(fā)育關系重建。鄰接法主要執(zhí)行參數(shù)為：Bootstrap Replications=10000，Substitution Model=Kimura 2-parameter model。最大簡約算法主要執(zhí)行參數(shù)為：collapse=maxbrlen，opt=acctran，hsearch start=stepwise，addseq=random，nreps=20，swap=TBR，multrees=yes，steepest=yes，bootstrap nreps=1000。貝葉斯算法主要執(zhí)行參數(shù)為：lset nst=6，rates=invgamma，prset revmat=dirichlet(1,8,1,1,8,1)，statefreq(7,2,2,9)，mcmcp ngen=5000000，nruns=2，nchains=4，temp=0.1，burninfrac=0.3。

2 結果與分析

2.1 DNA序列變異

針對9個地理種群的237只東方蜜蜂樣本，本試驗測序獲得了237條可用的線粒體DNA 序列，序列長度700～702 bp。通過與NCBI數(shù)據(jù)庫中已公布序列進行比對確認為目標片段CO1-2。序列的A+T平均含量為83.4%，與譚宏偉報道的東方蜜蜂線粒體基因組A+T含量(83.96%)[28]相近。序列比對結果共發(fā)現(xiàn)30個變異位點，其中單一多態(tài)位點4個，簡約信息位點26個，發(fā)現(xiàn)了兩處1-2個堿基的插入/缺失。30個變異位點中包括24個轉(zhuǎn)換位點5個顛換位點和1個既轉(zhuǎn)換又顛換的位點，共形成35種單倍型，依次命名為T1-T35，35種單倍型的序列變異樣式見附件。與24條參考序列合并進行分析時同源片段縮短為441 bp，序列比對發(fā)現(xiàn)59個變異位點，形成了30種單倍型。其中從樣本237條序列中探測到17種單倍型，依次命名為Tib1-Tib17；從參考序列中共發(fā)現(xiàn)16種單倍型，3種同前(Tib1、Tib8、Tib10)，13種不同的單倍型依次命名為Ref1-Ref13。

2.2 地理種群內(nèi)遺傳多樣性

9個地理種群中(除QHMH外)均發(fā)現(xiàn)一定比例(≥25%)的獨有單倍型，全部樣本的獨有單倍型率為80%，而XZLZ地理種群的獨有單倍型率為100%(表2)，表明各地理種群遺傳結構的獨特性較為顯著；核苷酸多樣度(π)指從群體中隨機抽取的兩條序列之間平均每個位點上的核苷酸差異數(shù)[29]，常用于量度種群內(nèi)包含的遺傳多樣性。當一個群體內(nèi)包含較多的相互間序列差異較大的個體時其π值相對較高，反之則π值相對較低。本次試驗中全部樣本的π值為0.003 92，比之前研究報道中廣東、廣西及湖南地區(qū)[20,22]東方蜜蜂種群的π值高，推測青藏高原地區(qū)可能蘊藏著比華南平原地區(qū)更豐富的遺傳多樣性；Tajima’s D中性檢驗結果顯示沒有任何種群呈現(xiàn)顯著或極顯著的正/負D值，各地理種群近期可能未曾經(jīng)歷明顯的瓶頸效應或種群擴張。

表2 東方蜜蜂各地理種群遺傳結構指標分析結果Table 2 Results of genetic diversity analysis for Apis cerana geographic populations

注：N樣本量；Nh單倍型數(shù)；U 獨有單倍型數(shù)；U % 獨有單倍型率；Hd 單倍型多樣度；π 核苷酸多樣度；D Tajima中性檢驗D值；*P≤0.05；**P≤ 0.01

Notes: N sample size; Nh number of haplotypes; U number of unique haplotypes; U% percent of unique haplotypes; Hd haploytpe diversity; π nucleotide diversity; D Tajima’s D test; *P≤0.05; **P≤ 0.01

2.3 青藏高原地區(qū)種群遺傳結構

方差分析結果顯示：不論如何分組，群體內(nèi)變異對總變異的貢獻率始終最大(>45.5%)且變異固定指數(shù)均達到極顯著水平(表3)。其中第5號分組方式的組間變異對總變異的貢獻率相對最高(38.42%)且變異固定指數(shù)達到顯著水平，表明按5號分組方式劃分遺傳類群相對最合理，可見不同地理種群間存在著較顯著的遺傳差異；單倍型網(wǎng)絡關系圖(圖2)顯示：共探測到7種共享單倍型，占全部單倍型種類(35種)的20%。且共享單倍型最多只出現(xiàn)于3個地理種群，未探測到廣泛分布于9個地理種群的共享單倍型。表明青藏高原東方蜜蜂各地理種群間的基因交流十分有限。同一地理種群包含的不同單倍型間呈現(xiàn)相對較大的遺傳差異，表現(xiàn)出較顯著的種群內(nèi)變異，這與方差分析的結果一致；Mantel檢驗結果顯示：9個東方蜜蜂地理種群的遺傳距離(Fst)矩陣與地理距離矩陣的相關性系數(shù)為0.113(P=0.259)，未達到顯著水平，據(jù)此推測距離隔離并不是形成青藏高原地區(qū)東方蜜蜂種群遺傳結構的主要因素。

表3 分子方差分析結果Table 3 Results of analysis of molecular variance (AMOVA)

注：分組規(guī)則，按行政區(qū)劃或地理距離遠近分組；Fct組間變異固定指數(shù)；Fsc組內(nèi)群體間變異固定指數(shù)；Fst群體內(nèi)變異固定指數(shù)；*P≤0.05；**P≤0.01。

Note: grouping rules, grouping according to administrative division or geographic distance;Fctfixation index of variation among groups;Fscfixation index of variation among populations within group;Fstfixation index of variation within populations; *P≤ 0.05;**P≤ 0.01.

圖2 35種線粒體DNA單倍型的網(wǎng)絡關系圖Fig.2 Network relationship of the 35 mtDNA haplotypes注：不同顏色代表不同地理種群；圓圈代表推測存在但未被檢測到的單倍型；方點代表突變步驟。Note：different colors represent different populations; circles represent inferred haplotypes that are detected here; square dots represent mutation steps.

與參考序列合并進行分析時獲得了30種單倍型，其系統(tǒng)發(fā)育關系如圖3所示：三種算法推斷出的系統(tǒng)發(fā)育樹拓撲結構完全一致，各大分支的支持率均大于60%。青藏高原地區(qū)東方蜜蜂與參考序列中云南、廣西、日本、韓國及一條來自印度(KF889329，Ref-6)的東方蜜蜂聚為一支，彼此間遺傳差異相對較小；與來自臺灣島、馬來西亞、文萊、印度(其它)、俄羅斯、泰國等地區(qū)的參考序列分支支持率較高，呈現(xiàn)出相對較大的遺傳差異。整體而言，單倍型間系統(tǒng)發(fā)育關系呈現(xiàn)顯著的地理分布格局。

3 討論與結論

本文采用的線粒體DNA片段長度約700 bp，從青藏高原地區(qū)東方蜜蜂樣本中共檢測出30個變異位點，占序列總長度的(30/700)4.28%。能夠較充分地反映地理種群內(nèi)/間存在的變異，從而有利于分析種群遺傳結構、評價遺傳分化程度。利用該片段探測到全部樣本的單倍型多樣度(Hd)和核苷酸多樣度(π)分別為0.939和0.00392，表明青藏高原地區(qū)東方蜜蜂種群中蘊藏著較豐富的遺傳多樣性。

圖3 a:30種單倍型的鄰接法系統(tǒng)發(fā)育樹；b:24條參考序列信息Fig.3 a:Neighbor-joining phylogenetic tree of the 30 haplotypes; b:Information about the 24 reference sequences(分支支持率依次為鄰接法、最大簡約法和貝葉斯法的自檢舉支持率)注：Tib1, Tib8, Tib10三種單倍型既出現(xiàn)于實驗樣本又出現(xiàn)于參考序列中；Ref單倍型只出現(xiàn)于參考序列中(Branch support values are bootstrap values from Neighbor-joining、Maximum parsimony and Bayesian inference in turns)Note: Tib1, Tib8 and Tib10 present in both samples of present study and reference sequences; Ref haplotypes present only in reference sequences

西藏地區(qū)XZDJ、XZLZ和XZMK采樣點檢測到的一些單倍型(T23、T25、T27、T29、T33)與四川SCMEK采樣點所檢測到的單倍型(T14和T15)相互間遺傳差異較大，這兩個區(qū)域恰是東方蜜蜂西藏亞種和阿壩亞種(形態(tài)類群)[16,18]的分布區(qū)域。這可能是兩個亞種間呈現(xiàn)遺傳分化的分子證據(jù)，對此還有待進一步開展細致的采樣與研究。

分子方差分析(表3)結果表明：種群內(nèi)變異在總變異中所占比重始終最大。Network圖(圖2)亦顯示：地理種群內(nèi)包含著相互間遺傳差異較大的單倍型，GSMX和SCMEK兩個種群尤其如此。經(jīng)與當?shù)毓ぷ魅藛T交流了解到：近年來當?shù)卦徺I了一些外地的東方蜜蜂蜂群，這些外來蜂群的混入在一定程度上增加了當?shù)貣|方蜜蜂種群的遺傳多樣度。這種人為基因交流導致的種群遺傳結構“復雜化”現(xiàn)象也不同程度地影響著其它地區(qū)的東方蜜蜂。

青藏高原東方蜜蜂不同地理種群間的共享單倍型較少，且這些共享單倍型最多只出現(xiàn)于3個地理種群(圖2)，未發(fā)現(xiàn)廣泛分布于9個地理種群的共享單倍型。表明該地區(qū)東方蜜蜂地理種群間的基因交流十分有限，各種群在近期經(jīng)歷了相對獨立的繁衍過程。青藏高原是地球上海拔最高的高原，平均海拔4 000～5 000 m。高原上絕大部分區(qū)域并不適宜中蜂生存，只在海拔相對較低且有顯花植物分布的區(qū)域生活著東方蜜蜂[16]。本試驗的樣本均采集自海拔2 200～2 700 m的區(qū)域，這些區(qū)域被連綿的高海拔山脈、高原所分隔。這在很大程度上限制了自然分蜂群的飛行范圍，阻隔了不同地區(qū)間蜂群的自然基因交流。此類地形地貌有利于不同地理種群間遺傳分化的加劇和不同亞種(形態(tài)類群)的形成，東方蜜蜂西藏亞種和阿壩亞種(形態(tài)類群)的形成正是“得益于此”。由此可見：山脈、高原阻隔是塑造青藏高原地區(qū)東方種群遺傳結構的主要作用因素。

經(jīng)與參考序列合并分析，單倍型系統(tǒng)發(fā)育關系(圖3)顯示：青藏高原地區(qū)與云南、廣西及韓國、日本的東方蜜蜂系譜關系相對較近；與臺灣島、馬來西亞、文萊、印度、泰國、俄羅斯等地區(qū)的系譜關系相對較遠。東方蜜蜂的系統(tǒng)發(fā)育關系呈現(xiàn)顯著的地理分布格局，這與我們之前研究[20]的結論是一致的。

本試驗亦存在不足之處：(1) 在與參考序列合并分析時同源片段縮短，圖3系統(tǒng)發(fā)育樹是基于441 bp線粒體DNA片段構建的，一些變異信息可能因此而被“隱藏”；(2) 數(shù)據(jù)分析中中國大陸其它地區(qū)的參考序列相對缺乏，以至于對比分析不夠全面。這些有待后續(xù)研究進一步加以完善！