趙　芳　陳振華　李文靖　張同作　皮　立　林恭華（中國(guó)科學(xué)院西北高原生物研究所，西寧8000；2青海中科蜜蜂研究院，西寧8000；3青海青藏華峰中蜂蜂業(yè)有限公司，貴德8700）

青海東方蜜蜂貴德種群的遺傳分化地位

趙芳1，2陳振華2，3李文靖1，2張同作1皮立1林恭華1，2
（1中國(guó)科學(xué)院西北高原生物研究所，西寧810001；2青海中科蜜蜂研究院，西寧810001；3青海青藏華峰中蜂蜂業(yè)有限公司，貴德811700）

對(duì)青海東方蜜蜂貴德種群（QHGD）進(jìn)行全基因組測(cè)序，并結(jié)合現(xiàn)有其他地區(qū)東方蜜蜂的分子序列信息，分析QHGD種群的遺傳分化地位。測(cè)序、組裝和比對(duì)后得到Cytb、tRNAIle～ND2、16S rRNA、EF1-α、COI～COII序列對(duì)應(yīng)的比對(duì)文件。單倍型分析顯示，QHGD種群的Cytb、EF1-α、COI～COII序列與國(guó)內(nèi)其他種群有重疊，但tRNAIle～ND2、16S rRNA序列則形成特有的單倍型?；趖RNAIle～ND2、16S rRNA、EF1-α合并序列的Network分析顯示，QHGD種群與四川九寨（SCJZ）、四川馬爾康（SCMEK）等地方種群成輻射狀與一個(gè)median vector節(jié)點(diǎn)（mv1）連接。本研究的結(jié)果表明，青海貴德種群是較為獨(dú)特的一個(gè)進(jìn)化分支單元，需要在未來(lái)的保護(hù)和開(kāi)發(fā)利用過(guò)程中予以關(guān)注。

東方蜜蜂；貴德種群；遺傳分化；基因組

青海東方蜜蜂是優(yōu)良的中華蜜蜂地方品種之一，主要分布于青海省東部農(nóng)業(yè)區(qū)。青海東方蜜蜂在分類(lèi)上屬于東方蜜蜂阿壩亞種（Apis cerana abansis）［1］，在形態(tài)和行為等方面都表現(xiàn)出明顯的高海拔、低氣溫的適應(yīng)特征［2，3］，對(duì)青藏高原地區(qū)養(yǎng)蜂業(yè)發(fā)展和生態(tài)平衡維持具有重要作用。長(zhǎng)期以來(lái)，青海東方蜜蜂的養(yǎng)殖多是農(nóng)戶自發(fā)的散養(yǎng)模式，沒(méi)有形成產(chǎn)業(yè)化。由于繁育體系不到位，致使種群數(shù)量不斷衰減，加上外來(lái)西方蜜蜂的威脅，青海東方蜜蜂一度處于頻臨滅絕的邊緣［4，5］。

值得一提的是，當(dāng)前對(duì)青海東方蜜蜂的研究，其對(duì)象都是青海省民和縣的種群［5］，而對(duì)其他地區(qū)的種群則很少涉及。青海省貴德縣蜂農(nóng)自1979年開(kāi)始利用當(dāng)?shù)氐囊吧淙禾剿髑嗪|方蜜蜂的科學(xué)繁育技術(shù)，并于2010年取得突破。以青海青藏華峰中蜂蜂業(yè)有限公司為平臺(tái)，現(xiàn)已繁育優(yōu)良生產(chǎn)群（命名為貴德種群）1，500群，生產(chǎn)性能表現(xiàn)優(yōu)異，其產(chǎn)品獲得第44屆APIMONDIA國(guó)際養(yǎng)蜂大會(huì)暨博覽會(huì)金牌和銀牌各一枚［6］。為了更好地保護(hù)和利用這一優(yōu)良品種資源，有必要對(duì)其遺傳特征進(jìn)行科學(xué)描述。本研究對(duì)青海東方蜜蜂貴德種群進(jìn)行基因組重測(cè)序，利用分子生物學(xué)研究方法，探討其遺傳分化地位。

1　材料方法

1.1樣品采集和高通量測(cè)序

蜜蜂樣品采自青海青藏華峰中蜂蜂業(yè)有限公司養(yǎng)蜂場(chǎng)內(nèi)的健康蜂群，隨機(jī)選取東方蜜蜂蜂群2群（QHGD1、QHGD2），用50 ml離心管從蜂箱內(nèi)采集10只工蜂，放入液氮中固定保存，用干冰運(yùn)輸?shù)奖本┌龠~客生物科技有限公司進(jìn)行高通量測(cè)序。實(shí)驗(yàn)流程按照Illumina公司提供的標(biāo)準(zhǔn)protocol執(zhí)行：每群取3只工蜂樣品，分別提取基因組DNA，檢測(cè)合格后用超聲波將DNA片段化；對(duì)片段化的DNA進(jìn)行片段純化、末端修復(fù)、3’端加A、連接測(cè)序接頭；用瓊脂糖凝膠電泳進(jìn)行片段大小選擇，對(duì)序列長(zhǎng)度為500 bp的片段進(jìn)行PCR擴(kuò)增形成測(cè)序文庫(kù)；質(zhì)檢合格的文庫(kù)用Illumina HiSeqTM 2500進(jìn)行雙向測(cè)序（pair-end），計(jì)劃測(cè)序深度為40X。對(duì)測(cè)序得到的原始測(cè)序序列（raw reads）進(jìn)行數(shù)據(jù)評(píng)估、過(guò)濾（去除接頭、低質(zhì)量區(qū)域等），最終得到用于生物信息學(xué)分析的clean reads。

1.2分子標(biāo)記序列組裝

目前，對(duì)國(guó)內(nèi)東方蜜蜂進(jìn)行遺傳分化研究，有3篇文獻(xiàn)采樣范圍較廣且給出樣品的來(lái)源信息［7-9］。將貴德生產(chǎn)群的對(duì)應(yīng)序列和這些研究進(jìn)行比較，可以得到貴德種群的遺傳分化地位。這些研究所涉及的分子標(biāo)記包括4個(gè)線粒體DNA片段（Cytb、tRNAIle～ND2、16S rRNA、COI～COII）和1個(gè)核基因片段（EF1-α）。首先，從GenBank中下載已經(jīng)測(cè)出的東方蜜蜂線粒體全基因組（登記號(hào)：NC_014295.1）［10］，以此為餌，用SOAPaligner軟件［11］在clean reads數(shù)據(jù)庫(kù)中提取匹配的reads。然后，將這些reads與此線粒體全基因組合并成單個(gè)文件，用Geneious軟件（Biomatters Limited）中的Assembly功能模塊，以此線粒體全基因組為模板將這些reads組裝成長(zhǎng)片段（scaffold）。結(jié)合使用blast+［12］和MEGA［13］兩個(gè)軟件，以文獻(xiàn)報(bào)道的EF1-α序列（GenBank登記號(hào)：KF663570.1）為餌，從東方蜜蜂全基因組序列文件中［14］找到包含 EF1-α序列的 scaffold（A_cerana_Scaffold_0111），并截取EF1-α片段及其側(cè)翼500 bp的片段。以同樣方法，用SOAPaligner和Geneious軟件組裝相應(yīng)片段。

1.3遺傳分化地位分析

從GenBank中下載前面提到的3篇文獻(xiàn)中所用序列［7-9］，用MEGA軟件分別將這些序列和本研究中組裝得到的貴德種群序列進(jìn)行比對(duì)，去除未匹配區(qū)段，最終得到基于不同序列片段的alignment文件。用DAMBE軟件分析變異位點(diǎn)數(shù)和單倍型組成情況，分析過(guò)程考慮插入/缺失位點(diǎn)。由于田崇浩等［8］涉及2個(gè)線粒體DNA片段（tRNAIle～ND2和16S rRNA）和1個(gè)核基因片段（EF1-α），涉及的變異位點(diǎn)信息較多，我們將3個(gè)片段合并，用DNasp軟件生成rdf文件（考慮插入/缺失位點(diǎn)）。用Network軟件計(jì)算單倍型關(guān)系（Median Joining模型），并繪制單倍型網(wǎng)絡(luò)圖。

2　結(jié)果與分析

2.1基因組測(cè)序和目標(biāo)片段組裝

對(duì)兩群蜜蜂樣品（QHGD1、QHGD2）成功進(jìn)行基因組測(cè)序，分別得到84，519，894和85，325，442條raw reads序列，兩樣本Q20值（質(zhì)量值大于等于20的堿基占總堿基數(shù)的百分比）都在90%和以上。經(jīng)過(guò)濾后分別得到81，629，314和84，873，218條clean reads，以韓國(guó)東方蜜蜂［14］為參考基因組，兩樣本mapped指數(shù)（定位到參考基因組的clean reads數(shù)占所有clean reads數(shù)的百分比）都高于78%，平均測(cè)序深度都高于40X，Coverage_ratio_5X值（5X深度及以上的堿基數(shù)占參考基因組總堿基數(shù)的比例）都高于90%（見(jiàn)表1）。以上結(jié)果顯示，兩份基因組測(cè)序質(zhì)量較好，可以滿足序列分析要求。

表1　青海東方蜜蜂貴德生產(chǎn)群基因組測(cè)序基本情況統(tǒng)計(jì)

SOAPaligner軟件比對(duì)顯示，QHGD1和QHGD2分別有939，283和716，069條clean reads與參考線粒體基因組匹配。Geneious軟件統(tǒng)計(jì)顯示，平均堿基覆蓋度分別為7，445X和5，677X，經(jīng)組裝后分別得到15，465 bp和15，462 bp長(zhǎng)度的單條scaffold，約占參考線粒體基因組長(zhǎng)度（15，897 bp）的97%。同時(shí)，QHGD1和QHGD2分別有444和361條clean reads與參考EF1-α及其側(cè)翼500 bp片段匹配；平均堿基覆蓋度分別為43X和35X，經(jīng)組裝后分別得到1，304 bp和1，286 bp長(zhǎng)度的單條scaffold，占參考EF1-α及其側(cè)翼500 bp片段長(zhǎng)度（1，306bp）的98%以上。

圖1　四川九寨（SCJZ，GenBank No:KF663582.1）群和青海貴德群（QHGD）的tRNAIle～ND2序列比對(duì)結(jié)果

圖2　四川九寨（SCJZ，GenBank No:KF572602.1）群和青海貴德群（QHGD）16S rRNA序列比對(duì)結(jié)果

2.2序列比對(duì)和遺傳分化地位分析

將本研究得到的線粒體或EF1-α及其側(cè)翼片段scaffold分別與已經(jīng)發(fā)布的 Cytb、tRNAIle～ND2、16S rRNA、EF1-α、COI～COII序列合并，比對(duì)后截取同源序列，最終得到5個(gè)alignment文件。DAMBE軟件的分析顯示，QHGD1和QHGD2兩個(gè)群的序列在上述5個(gè)基因區(qū)段完全相同，為敘述方便，將兩者視為同一個(gè)基因序列對(duì)待，統(tǒng)稱為QHGD（即“青海貴德”4個(gè)字的拼音首字母縮寫(xiě)）序列。

基于Cytb序列的分析顯示，QHGD序列與甘肅天水（GSTS，EF180093.1）、福建福州（FJFZ，F(xiàn)J229471.1）、四川宜賓（SCYB，F(xiàn)J229473.1）、山西沁源（SXQY，F(xiàn)J229480.1）種群的序列［7］完全一致?；趖RNAIle～ND2、16S rRNA序列的分析顯示，QHGD序列與現(xiàn)有的所有其他19種群［8］都不相同（圖1～2）?；贓F1-α序列的分析顯示，QHGD序列與絕大多數(shù)種群（JLAT、JLHD、JLHN、SDFX、SXZQ、SXTG、SXQY、SXQS、GSTS、SCMEK、YNBS、YNKM、YNXSBN、GXQZ、FJFZ、HNHK，KF663570.1）所含有的單倍型［8］完全相同?；贑OI～COII序列的分析顯示，QHGD序列與 Japan1 （AB078733.1）單倍型［9］完全相同。

基于tRNAIle～ND2、16S rRNA、EF1-α合并序列的Network分析顯示，QHGD種群與國(guó)內(nèi)其他19個(gè)地方種群沒(méi)有直接交集。相反，QHGD與四川九寨（SCJZ）、四川馬爾康（SCMEK）、山西左權(quán)（SXZQ）、山西沁源（SXQY）、山東費(fèi)縣（SDFX）共同指向一個(gè)中值向量（median vector）節(jié)點(diǎn)mv1（圖3），即這些種群為同一個(gè)祖先節(jié)點(diǎn)輻射狀分化而成。

圖3　基于tRNAIle～ND2、16S rRNA、EF1-α合并序列的Network關(guān)系圖

綠色數(shù)字表示兩兩節(jié)點(diǎn)間的差異位點(diǎn)數(shù)，未標(biāo)注的表示1個(gè)位點(diǎn)差異；紅色字母QHGD表示青海貴德生產(chǎn)群。字母代碼所表示的地點(diǎn)信息與田崇浩（2014）一致：JLAT，吉林安圖；JLHD，吉林樺甸；JLHN，吉林輝南；SDFX，山東費(fèi)縣；SXZQ，山西左權(quán)；SXTG，山西太谷；SXQY，山西沁源；SXQS，山西沁水；GSTS，甘肅天水；SCJZ，四川九寨；XZMEK，四川馬爾康；SCYB，四川宜賓；YNBS，云南保山；YNKM，云南昆明；YNXSBN，云南西雙版納；GDGZ，廣東廣州；GXQZ，廣西欽州；FJFZ，福建福州；HNHK，海南?？?。

3　討論

由于地處偏僻和產(chǎn)業(yè)化水平較低等原因，長(zhǎng)期以來(lái)，學(xué)界對(duì)青海地區(qū)的東方蜜蜂種群關(guān)注較少。楊冠煌等［1］首先提出，青海東部地區(qū)的東方蜜蜂屬于阿壩亞種，這一說(shuō)法在《中國(guó)畜禽遺傳資源志·蜜蜂志》中也得到體現(xiàn)。然而這一論述是基于青海省和四川阿壩地區(qū)地理距離和地貌特征相近而進(jìn)行的推測(cè)，并無(wú)直接研究證據(jù)支持。李華等［2］基于形態(tài)測(cè)量數(shù)據(jù)的研究表明，青海民和地區(qū)的東方蜜蜂與阿壩九寨種群有明顯不同。本研究的結(jié)果顯示，QHGD與阿壩地區(qū)的SCJZ和SCMEK種群從同一個(gè)節(jié)點(diǎn)輻射狀分化形成不同的進(jìn)化枝，換句話說(shuō)，青海貴德和阿壩地區(qū)的東方蜜蜂種群在遺傳上并無(wú)隸屬關(guān)系。以SCJZ為例，青海貴德種群在tRNAIle～ND2和16S rRNA序列上與其分別具有4 和2個(gè)位點(diǎn)差異（圖1～2）。此外，從Network關(guān)系上看，SCJZ和SCMEK與mv1節(jié)點(diǎn)僅有1個(gè)位點(diǎn)差異，而QHGD種群則與mv1節(jié)點(diǎn)間的位點(diǎn)差異達(dá)6個(gè)之多。以上結(jié)果表明，青海貴德種群可能不屬于阿壩亞種，相反，是較為獨(dú)特的一個(gè)進(jìn)化分支單元。

Tan et al［9］是目前已發(fā)表的唯一涉及青海東方蜜蜂（民和種群）的報(bào)道。其研究結(jié)果顯示，在COI～COII序列上，青海民和種群的單倍型（Qinghai，GenBank No：KT174438）和四川的部分種群?jiǎn)伪缎拖嗤?。有趣的是，本研究中青海貴德的COI～COII序列與Tan et al［9］文中的 Japan1（AB078733） 完全相同，與 Qinghai （KT174438）反而有1個(gè)堿基位點(diǎn)的差異。這是否意味著青海民和和貴德兩個(gè)東方蜜蜂種群有本質(zhì)上的差異，有待進(jìn)一步研究。

由于各種條件的限制，東方蜜蜂遺傳分化的研究有些混亂。形態(tài)學(xué)研究過(guò)程中，會(huì)帶入相對(duì)較多的主觀意見(jiàn)，而結(jié)合分子生物學(xué)研究手段將大大降低這種主觀性［15］。當(dāng)前，絕大多數(shù)涉及東方蜜蜂遺傳分化的研究，僅使用1個(gè)或少數(shù)幾個(gè)序列片段，由于缺乏足夠的變異信息，不同地理種群之間單倍型交疊，導(dǎo)致結(jié)果不明朗。以本研究為例，青海貴德種群在3個(gè)基因序列（Cytb、EF1-α、COI～COII）上都分別與多個(gè)其他地理種群有單倍型重疊，無(wú)法據(jù)此確定彼此的遺傳分化地位。相反，由于tRNAIle～ND2、16S rRNA、EF1-α屬于同一批樣品，將其合并后分析，可以得到相對(duì)明晰的結(jié)論，足見(jiàn)多個(gè)基因標(biāo)記方法的優(yōu)越性。本研究對(duì)青海東方蜜蜂貴德種群進(jìn)行基因組測(cè)序，得到巨量的序列信息，可以與現(xiàn)有其他地區(qū)的遺傳分化研究數(shù)據(jù)做到較好的對(duì)接，為闡明其遺傳分化地位并服務(wù)于其保護(hù)和利用提供重要的科學(xué)依據(jù)。

［1］楊冠煌，許少玉，匡邦郁，等.東方蜜蜂Apis cerana Fab.在我國(guó)的分布及其亞種分化 ［J］.云南農(nóng)業(yè)大學(xué)學(xué)報(bào)，1986，1（1）∶89-92.

［2］李華，張祖蕓，譚墾.青海東方蜜蜂的形態(tài)學(xué)研究 ［J］.蜜蜂雜志，2008，12∶5-6.

［3］李良德，馬少榮，白斌.青海東方蜜蜂遺傳資源調(diào)查 ［J］.中國(guó)蜂業(yè)，2011，62（1）∶26-27.

［4］張學(xué)功，張學(xué)成，馬少榮.青海中蜂養(yǎng)殖現(xiàn)狀及保護(hù)措施［J］.中國(guó)蜂業(yè)，2010，6（9）∶34-35.

［5］冶兆平.再議青海東方蜜蜂養(yǎng)殖現(xiàn)狀與保護(hù)措施 ［J］.中國(guó)蜂業(yè)，2014，65∶64-66.

［6］王建梅，陳黎紅，胡玭玭.中國(guó)養(yǎng)蜂學(xué)會(huì)蜂業(yè)代表團(tuán)在第44 屆APIMONDIA國(guó)際養(yǎng)蜂大會(huì)暨博覽會(huì)榮獲多項(xiàng)獎(jiǎng)牌［J］.蜜蜂雜志，2015，11∶45.

［7］高鵬飛，趙慧婷，張春香，等.基于線粒體Cyt b基因部分序列的中國(guó)東方蜜蜂不同地理種群的系統(tǒng)發(fā)育 ［J］.動(dòng)物學(xué)報(bào)，2008，54（6）∶1005-1013.

［8］田崇浩.中國(guó)境內(nèi)東方蜜蜂群體亞分化研究 ［D］.山西農(nóng)業(yè)大學(xué)碩士學(xué)位論文，2014.

［9］Tan K，Qu Y，Wang Z，et al.Haplotype diversity and genetic similarity among populations of the Eastern honey bee from Himalaya-Southwest China and Nepal（Hymenoptera∶Apidae）［J］.Apidologie，2015，DOI∶10.1007/s13592-015-0390-x.

［10］Tan HW，Liu GH，Dong X，et al.The complete mitochondrial genome of the Asiatic cavity-nesting honeybee Apis cerana（Hymenoptera∶Apidae）［J］.PLoS One，2011，6（8）∶e23008.

［11］Gu S，F(xiàn)ang L，Xu X.Using SOAPaligner for short reads alignment［J］.Current Protocols in Bioinformatics，2013，44∶11.11.1-11.11.17.

［12］Zhang Z，Schwartz S，Wagner L，et al.A greedy algorithm for aligning DNA sequences［J］.Journal of Computational Biology，2000，7（1-2）∶203-214.

［13］Tamura K，Stecher G，Peterson D，et al.MEGA6∶molecular evolutionary genetics analysis version 6.0［J］.Molecular Biology and Evolution，2013，30∶2725-2729.

［14］Park D，Jung JW，Choi BS，et al.Uncovering the novel characteristics of Asian honey bee，Apiscerana，by whole genome sequencing［J］.BMC Genomics，2015，16∶1.

［15］楊潔，和紹禹，苗永旺，等.東方蜜蜂（Apis cerana）的分類(lèi)研究進(jìn)展.蜜蜂雜志，2011，3∶13-15.

Genetic variation status of Apis cerana population from Guide，Qinghai

Zhao Fang1，2Chen Zhenhua2，3Li Wenjing1，2Zhang Tongzuo1Pi Li1Lin Gonghua1，2
（1 Northwest Institute of Plateau Biology，Chinese Academy of Sciences，Xining 810001；2 Qinghai-CAS Institute of Apicultural Research，Xining 810001；3 Qinghai Qingzang-Huafeng Eastern Honeybee Industry Co.，Ltd，Guide 811700）

We sequenced the wholegenome of Apis cerana from Guide population（QHGD）and，combining with sequence data on A.cerana populations from other regions，we analyzed the genetic variation status of the QHGD population.After sequencing，assembling，and aligning，we obtained alignments for five DNA fragments including Cytb，tRNAIle～ND2，16S rRNA，EF1-α，and COI～COII.Haplotype analysis showed that the QHGD population shared haplotypes with other populations on Cytb，EF1-α and COI～COII，in contrast，the QHGD population formed a unique haplotype on tRNAIle～ND2 and 16S rRNA，respectively.Network analysis based on a combination of tRNAIle～ND2，16S rRNA and EF1-αshowed that the populations of QHGD，SCJZ，SCMEK，etc.formed a radialshape from a median vector （mv1）.Our study indicated that the QHGD population was a very unique genetic branch，and people should pay attention to their special nature in the future conservation and developing processes.

Apis cerana；Guide population；genetic variation；genome

中國(guó)科學(xué)院青年創(chuàng)新促進(jìn)會(huì)人才項(xiàng)目（NO.2015352）

林恭華，副研究員，碩導(dǎo)，主要從事青藏高原動(dòng)物生態(tài)學(xué)研究，E-mail:lingonghua@nwipb.cas.cn