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

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

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

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

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

1 材料與方法

1.1 樣本采集與DNA提取

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

1.2 PCR擴(kuò)增與DNA測(cè)序

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

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

利用Lasergene 7.0軟件包中的SeqMan對(duì)正、反向DNA序列進(jìn)行檢查、校正、修剪并生成一致序列；從NCBI數(shù)據(jù)庫(kù)下載了24條其它地區(qū)東方蜜蜂的CO1-2片段作為參考序列(見圖3b)，利用MEGA 5.05軟件[23]對(duì)本次試驗(yàn)測(cè)序獲得的CO1-2一致序列和加入24條參考序列后的全部序列分別進(jìn)行Alignment計(jì)算(在默認(rèn)參數(shù)下)；使用DnaSP 5.1軟件[24]進(jìn)行單倍型種類、頻率、單倍型多樣度(Hd)、核苷酸多樣度(π)以及種群間遺傳距離(Fst)矩陣的統(tǒng)計(jì)計(jì)算；以GeographicDistanceMatrixGenerator v1.2.3生成種群間地理距離矩陣；使用Arlequin v3.11軟件[25]進(jìn)行Tajima’s D中性檢驗(yàn)和分子方差分析(AMOVA)，并對(duì)種群間遺傳距離Fst矩陣與地理距離矩陣之間實(shí)施Mantel檢驗(yàn)；使用Network 5.03軟件(http://www.fluxus-engineering.com)中的Median-Joining算法對(duì)鑒別出的單倍型進(jìn)行網(wǎng)絡(luò)關(guān)系構(gòu)建；以蘇拉威西蜂(Apiskoschevnikovi)同源片段為外群并加入?yún)⒖夹蛄泻?，分別利用鄰接法(Neighbour-Joining)、最大簡(jiǎn)約法(Maximum-Parsimony)[26]和貝葉斯法(Bayesian-Algorithm)[27]對(duì)鑒別出的單倍型進(jìn)行系統(tǒng)發(fā)育關(guān)系重建。鄰接法主要執(zhí)行參數(shù)為：Bootstrap Replications=10000，Substitution Model=Kimura 2-parameter model。最大簡(jiǎn)約算法主要執(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 結(jié)果與分析

2.1 DNA序列變異

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

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

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

表2 東方蜜蜂各地理種群遺傳結(jié)構(gòu)指標(biāo)分析結(jié)果Table 2 Results of genetic diversity analysis for Apis cerana geographic populations

注：N樣本量；Nh單倍型數(shù)；U 獨(dú)有單倍型數(shù)；U % 獨(dú)有單倍型率；Hd 單倍型多樣度；π 核苷酸多樣度；D Tajima中性檢驗(yàn)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ū)種群遺傳結(jié)構(gòu)

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

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

注：分組規(guī)則，按行政區(qū)劃或地理距離遠(yuǎn)近分組；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)絡(luò)關(guān)系圖Fig.2 Network relationship of the 35 mtDNA haplotypes注：不同顏色代表不同地理種群；圓圈代表推測(cè)存在但未被檢測(cè)到的單倍型；方點(diǎn)代表突變步驟。Note：different colors represent different populations; circles represent inferred haplotypes that are detected here; square dots represent mutation steps.

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

3 討論與結(jié)論

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

圖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(分支支持率依次為鄰接法、最大簡(jiǎn)約法和貝葉斯法的自檢舉支持率)注：Tib1, Tib8, Tib10三種單倍型既出現(xiàn)于實(shí)驗(yà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采樣點(diǎn)檢測(cè)到的一些單倍型(T23、T25、T27、T29、T33)與四川SCMEK采樣點(diǎn)所檢測(cè)到的單倍型(T14和T15)相互間遺傳差異較大，這兩個(gè)區(qū)域恰是東方蜜蜂西藏亞種和阿壩亞種(形態(tài)類群)[16,18]的分布區(qū)域。這可能是兩個(gè)亞種間呈現(xiàn)遺傳分化的分子證據(jù)，對(duì)此還有待進(jìn)一步開展細(xì)致的采樣與研究。

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

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

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

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