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

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

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

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

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

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

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

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

1　材料方法

1.1樣品采集和高通量測序

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

1.2分子標記序列組裝

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

1.3遺傳分化地位分析

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

2　結果與分析

2.1基因組測序和目標片段組裝

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

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

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

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

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

2.2序列比對和遺傳分化地位分析

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

基于Cytb序列的分析顯示，QHGD序列與甘肅天水（GSTS，EF180093.1）、福建福州（FJFZ，FJ229471.1）、四川宜賓（SCYB，FJ229473.1）、山西沁源（SXQY，FJ229480.1）種群的序列［7］完全一致。基于tRNAIle～ND2、16S rRNA序列的分析顯示，QHGD序列與現有的所有其他19種群［8］都不相同（圖1～2）?；贓F1-α序列的分析顯示，QHGD序列與絕大多數種群（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種群與國內其他19個地方種群沒有直接交集。相反，QHGD與四川九寨（SCJZ）、四川馬爾康（SCMEK）、山西左權（SXZQ）、山西沁源（SXQY）、山東費縣（SDFX）共同指向一個中值向量（median vector）節(jié)點mv1（圖3），即這些種群為同一個祖先節(jié)點輻射狀分化而成。

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

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

3　討論

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

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

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

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

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

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

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

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

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

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

［8］田崇浩.中國境內東方蜜蜂群體亞分化研究 ［D］.山西農業(yè)大學碩士學位論文，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，Fang 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）的分類研究進展.蜜蜂雜志，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

中國科學院青年創(chuàng)新促進會人才項目（NO.2015352）

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