中華蜜蜂幼蟲腸道參考轉錄組的de novo組裝及SSR分子標記鑒定

徐細建，郭睿，駱群，熊翠玲，梁勤，張串聯，鄭燕珍，張曌楠，黃枳腱，張璐，李汶東，陳大福

（1福建農林大學蜂學學院，福州 350002；2江西省養蜂研究所，南昌 330201）

中華蜜蜂幼蟲腸道參考轉錄組的de novo組裝及SSR分子標記鑒定

徐細建1，郭睿1，駱群2，熊翠玲1，梁勤1，張串聯2，鄭燕珍1，張曌楠1，黃枳腱1，張璐1，李汶東1，陳大福1

（1福建農林大學蜂學學院，福州 350002；2江西省養蜂研究所，南昌 330201）

【目的】利用RNA seq技術對中華蜜蜂（Apis cerana cerana，簡稱中蜂）幼蟲腸道參考轉錄組進行de novo組裝，并進行功能及代謝通路注釋，進而利用該轉錄組數據進行中蜂幼蟲的SSR分子標記鑒定。【方法】實驗室條件下飼養中蜂幼蟲，將純化的蜜蜂球囊菌（Ascosphaera apis，簡稱球囊菌）孢子飼喂3日齡幼蟲，剖取 4、5和 6日齡幼蟲腸道，液氮速凍。將健康幼蟲腸道與感染球囊菌的幼蟲腸道同時進行 Illumina測序。通過對raw reads的過濾得到clean reads，利用Trinity軟件組裝得到unigenes。通過BLASTx（E-value＜10?5）比對NCBI Nr、Swiss-Prot、KOG和KEGG數據庫，對unigenes進行功能和代謝通路注釋。利用MISA軟件對所有unigenes進行SSR搜索，并利用Primer Premier 5軟件設計特異性SSR引物，通過常規PCR對來源于北京、遼寧興城和四川成都的中蜂幼蟲腸道樣本進行SSR位點鑒定。【結果】中蜂幼蟲腸道的RNA seq共得到35 670 000條reads，de novo組裝得到43 557個unigenes，平均長度為898 nt。共有18 225個unigenes可被注釋到上述公共蛋白數據庫，單獨注釋到NCBI Nr、Swiss-Prot、KOG和KEGG數據庫的unigenes數分別為3 899、443、37和10個。KOG注釋結果顯示，11 442條unigenes分布于25個基因家族，其中注釋到RNA加工和修飾家族的基因數最多，達1 249個。9 679個unigenes的GO分類結果顯示，在生物學進程分類中，注釋到細胞進程的基因最多，達4 201個，在分子功能和細胞組分類中，注釋到結合與細胞的基因數最多，分別為4 935和2 900個。4 517個unigenes可注釋到KEGG數據庫中的216個代謝通路，注釋到核糖體的基因數最多，達385個。利用MISA軟件，可在7 763個unigenes搜索到13 448個SSR位點，隨機選取20對SSR引物對國內3個不同來源的中蜂幼蟲腸道樣本的SSR位點進行擴增，有6對引物可鑒定出SSR分子標記。【結論】成功組裝并注釋了中蜂幼蟲腸道參考轉錄組，可為中蜂及其幼蟲的分子生物學及組學研究提供重要的參考信息，也可用于補充、豐富和檢驗東方蜜蜂的參考基因組，基于此轉錄組數據開發出6個中蜂的SSR分子標記，可應用于中蜂的基因圖譜構建、基因多樣性分析、基因定位等研究，也說明利用轉錄組數據開發非模式生物SSRs的方法可行。

RNA seq；參考轉錄組；中華蜜蜂；unigene；SSR

0 引言

【研究意義】蜜蜂是重要的授粉昆蟲和社會學模式昆蟲，因其對研究神經生物學、發育、社會行為和表觀基因組學的重要性而廣受關注[1-4]，西方蜜蜂（Apis mellifera）早在2006年就已完成基因組測序[5]，為研究蜜蜂行為、遺傳進化和基因功能提供了重要的信息和線索。與西方蜜蜂相比，東方蜜蜂（Apis cerana）更易適應極端天氣、飛行距離更長、具有更強的梳理行為和清潔行為以及群體防御能力[6-9]。中華蜜蜂（Apis cerana cerana，簡稱中蜂）長期進化適應本土環境，相比于西方蜜蜂具有抗螨害、耐寒、善利用零星蜜粉源等優點[8-12]。利用RNA seq技術對中蜂幼蟲腸道進行深度測序，de novo組裝其參考轉錄組并進行功能及代謝通路注釋，可為中蜂幼蟲的分子及組學研究提供重要參考信息，在此基礎上鑒定出的SSR分子標記可為在分子水平深入研究中蜂的重要性狀、復雜行為及遺傳進化提供寶貴信息。【前人研究進展】近年來，以RNA seq為代表的高通量測序技術發展迅猛，廣泛應用于動植物及微生物研究[13-19]，在蜜蜂研究方面也取得了一系列重要進展[20-21]。中國養蜂生產中的常見蜂種為意大利蜜蜂（Apis mellifera ligustica，簡稱意蜂）和中蜂。PARK等[22]通過對 A. cerana雄蜂的基因組測序和對A. mellifera及A. cerana工蜂多個組織的轉錄組測序，獲得238 Mb的基因組數據和10 651個基因信息，并對A. cerana特有基因進行了分析，但作者當時并未公布基因位置及功能注釋信息。SSR分子標記開發的傳統方法是通過構建基因組DNA文庫，成本昂貴且費時費力，而利用高通量測序技術的新一代SSR鑒定則更為經濟、高效[13-15]。梁勤等[23]利用6對微衛星DNA標記對福建省4個中蜂群體進行遺傳多樣性分析，評估群體內的遺傳變異和群體間的遺傳分化；徐新建等[24]應用10個微衛星DNA標記對海南島11個地點和大陸2個地點中蜂分析表明，海南中蜂多樣性豐富，島嶼和鄰近大陸種群發生了明顯的遺傳分化。目前，已開發的中蜂SSR分子標記很少[25-26]，制約了中蜂分子進化及種群遺傳學的發展。【本研究切入點】東方蜜蜂的基因組已完成測序并公布，但缺乏專一的中蜂幼蟲腸道參考轉錄組，嚴重制約中蜂幼蟲的病原-宿主互作及免疫應答研究。【擬解決的關鍵問題】利用RNA seq數據組裝并注釋中蜂幼蟲腸道參考轉錄組，并鑒定出若干SSR分子標記，解決中蜂幼蟲參考轉錄組缺失以及SSR分子標記較少的問題。

1 材料與方法

試驗于2015年12月至2016年8月在福建農林大學蜂學學院蜜蜂保護學實驗室完成。

1.1 供試材料

中蜂幼蟲取自福建農林大學蜂學學院教學蜂場，蜜蜂球囊菌（Ascosphaera apis）菌株保存于福建農林大學蜂學學院蜜蜂保護學實驗室。

1.2 主要試劑及儀器

RNase-free水購自中國上海生工生物公司；DNaseI和Oligotex mRNA Kits Midi試劑盒購自德國Qiagen公司；Dynal M280磁珠購自Invitrogen公司；高碘酸鈉購自美國Sigma公司；DNA ligase購自美國Thermo公司；RNA Reagent抽提試劑盒、Ex Taq polymerase及Superscript II reverse transcriptase均購自日本TaKaRa公司；純化cDNA的Ampure beads為美國 Agencourt產品；cDNA 文庫構建試劑盒TruSeqTMDNA Sample Prep Kit -Set A為美國Illumina公司產品。其他試劑均為國產分析純。

恒溫恒濕氣候箱購自中國寧波江南儀器廠；pH計購自中國上海儀電科學股份有限公司；超純水儀購自中國四川沃特爾水處理設備有限公司；高速冷凍離心機購自德國Eppendorf公司；倒置顯微鏡為中國上海光學儀器五廠產品；超凈工作臺為中國蘇州安泰空氣技術有限公司產品；PCR儀為美國Bio Rad公司產品；凝膠成像系統為中國上海培清科技有限公司產品；超低溫冰箱為中科美菱低溫科技股份有限公司產品。

1.3 方法

1.3.1 幼蟲的人工飼養 中蜂幼蟲的人工飼料參照王倩等[26]的方法配制并改良，將 D-果糖和 D-葡萄糖換為新鮮蜂蜜。預試驗中對照組中蜂幼蟲7日齡成活率達到70%以上。從福建農林大學蜂學學院蜂場選擇經PCR檢測為球囊菌陰性的健康蜂群。用滅菌的移蟲針挑取2日齡幼蟲，放入無菌的24孔細胞培養板（每孔對應1只幼蟲，孔內加有35℃預溫的幼蟲飼料），將24孔板放入恒溫恒濕培養箱，每隔24 h吸去舊飼料、加入新飼料。3日齡時，一組幼蟲飼喂含有球囊菌孢子（1×107孢子/mL）的人工飼料，另一組幼蟲飼喂以正常人工飼料。35℃，90% RH條件下飼喂幼蟲至7日齡，上述兩組幼蟲組均設置3個生物學重復。1.3.2 測序樣品準備 分別于4、5、6日齡剖取中蜂幼蟲腸道組織，為盡量減少腸道 RNA的降解，將從解剖取樣到樣品放入液氮速凍的時間控制在 30 s以內，每剖取一組樣品，液氮速凍后迅速放入-80℃超低溫冰箱保存。

1.3.3 cDNA文庫構建及 RNA seq 利用 RNAiso Reagent試劑盒抽提中蜂幼蟲腸道組織的總RNA，然后用RNase-free DNaseI去除基因組DNA殘留。RNA的質量通過瓊脂糖凝膠電泳和 NanoDrop ND-2000（NanoDrop，Wilmington，DE，USA）進行檢測。利用Oligotex mRNA Kits Midi試劑盒說明書，純化各樣品總RNA中的mRNA。以10 μg mRNA作為模板，GsuI-oligo dT作為反轉錄引物，用1 000 U Superscript II reverse transcriptase在42℃下孵育1 h合成第1鏈cDNA；隨后利用高碘酸鈉氧化mRNA的5′端帽子結構，并連接生物素；通過Dynal M280磁珠篩選連接了生物素的mRNA/cDNA，并通過堿裂解釋放第1鏈cDNA；然后通過DNA ligase在第1鏈cDNA 的5′末端加上接頭，利用Ex Taq polymerase合成第2鏈cDNA。最后，通過GsuI酶切去除polyA和5′端接頭。利用Ampure beads對上述cDNA進行純化，cDNA文庫通過TruSeqTMDNA Sample Prep Kit-Set A進行構建和TruSeq PE Cluster Kit進行擴增。委托廣州基迪奧生物技術有限公司對上述12個腸道樣品進行深度測序，測序平臺為Illumina Hiseq 2500，各腸道樣品的3個生物學重復均同時進行測序。

1.3.4 中蜂幼蟲腸道參考轉錄組的de novo組裝 首先，利用Perl腳本去除含有adaptor、未知核苷酸比例＞5%和低質量 reads，獲得 clean reads。利用軟件Trinity[27]進行中蜂轉錄組的 de novo組裝（缺省值Kmer=25）。長度短于200 bp的contigs和unigenes將被舍棄。過濾和組裝以后得到高質量的unigenes。

1.3.5 Unigenes注釋 利用 BLASTx（E-value＜10?5）將測序序列比對NCBI nr數據庫（http://www.ncbi. nlm.nih.gov）、Swiss-Prot 數據庫（http://www.expasy.ch/ sprot）、KOG（Clusters of orthologous groups for eukaryotic complete genomes）數據庫（ftp://ftp.ncbi.nih. gov/pub/COG/KOG/kyva）和KEGG代謝通路（pathway）數據庫（http://www.genome.jp/kegg/）。利用BLASTX將組裝出來的unigenes序列與Nr數據庫進行比對后，取每個unigenes在Nr庫中比對結果最好（E值最低）的那一條序列為對應同源序列（如有并列，取第一條）確定同源序列所屬物種，統計比對到各個物種的同源序列數量。基于Nr database注釋結果，利用Blast2GO進行unigenes的GO注釋，利用WEGO軟件對每一個轉錄本進行GO分類。

1.3.6 SSR分子標記開發 利用軟件 MISA（http:// pgrc.ipk-gatersleben.de/misa/）搜索unigenes的微衛星標記，按照以下標準從unigenes中查找SSR位點：二核苷酸重復≥6次，三核苷酸重復≥5次，四核苷酸重復≥5次，五核苷酸重復≥5次和六核苷酸重復≥5次。根據 MISA的輸出結果，利用 Primer Premier 5（PREMIER Biosofe Int.，Palo Alto，CA）對每一個含有16 bp堿基重復的SSR設計引物。

選取北京（B）、遼寧興城（L）、四川成都（S）3個不同來源的中蜂幼蟲腸道樣本作為模板，隨機選取20對SSR引物進行PCR擴增，PCR程序：94℃預變性5 min；94℃變性50 s，55℃退火30 s，72℃延伸30 s，共33個循環，72℃再延伸10 min。PCR產物經1%瓊脂糖凝膠電泳檢測。

2 結果

2.1 中蜂幼蟲腸道的RNA seq及參考轉錄組de novo組裝

對上述12個腸道樣品進行Illumina測序，平均得到30 584 420條原始讀段（raw reads），去除低質量和含有接頭的reads后平均獲得29 726 139條有效讀段（clean reads），總測序長度為3 715 767 396， Q20平均為98.31%，說明測序數據質量較好，可用于下一步分析。各樣品的測序詳細信息如附表1所示。

對 clean reads進行進一步序列拼接和去冗余處理，組裝得到43 557條unigenes，平均長度達898 nt， N50為1 704 nt（表1）。統計結果顯示，unigenes的數目隨著序列長度的增加而減少，在200—299 nt長度范圍內數目最多，符合生物體序列長度分布的基本規律。長度＞1 000 nt的unigenes有10 454條，占總unigenes的24.00%。上述結果說明中蜂幼蟲腸道的組裝質量較好。轉錄組測序數據已上傳NCBI SRA數據庫，SRA號：SRA456721。

2.2 Unigenes注釋

利用 BLASTx（E-value＜10?5）將測序序列比對NCBI Nr、Swiss-Prot、KOG和KEGG pathway數據庫，結果顯示分別有17 456、12 830、11 442和9 045個unigenes能夠注釋到上述數據庫，有功能或代謝通路注釋的unigenes數目為18 225，占全部unigenes的41.84%，此外，有58.16%的unigenes無功能注釋（表2）。有29個unigenes在4大數據庫均有注釋，而僅能注釋到NCBI Nr、Swiss-Prot、KOG和KEGG pathway數據庫的unigenes分別為3 899、443、37和10個。

表1 中蜂幼蟲腸道參考轉錄組組裝結果統計Table 1 Summary of A. c. cerana larval gut’s reference transcriptome assembled in this study

表2 公共蛋白數據庫注釋統計表Table 2 Summary of annotation information of all unigenes in public protein databases

注釋到Nr數據庫中unigenes的E-value分布顯示（圖1），比對到物種序列的E-value均＜10-5，其中E-value＜10-100的有49.76%，說明比對結果可信度較高。注釋基因同源序列的物種分布統計結果顯示前 10位的物種依次為 Apis mellifera、Apis dorsata、Apis florea、Bombus impatiens、Bombus terrestris、Lasius niger、Megachile rotundata、 Harpegnathos saltator、Capsaspora owczarzaki ATCC 30864和Cerapachys biroi，注釋到A. mellifera的基因數為5 753（31.57%），注釋到A. dorsata和A. florea的基因數分別為3 695（20.27%）和2 489（13.66%）（表3）。

KOG注釋結果顯示，11 442個unigenes分布于25個基因家族（圖2）。其中，注釋基因數最多的為信號轉導機制，其次為一般功能預測和翻譯后修飾、蛋白翻轉和分子伴侶。值得注意的是，有 170條 unigenes注釋到防御機制，它們可能在中蜂幼蟲抵御病原入侵過程發揮重要作用。

圖1 E值分布Fig.1 Distribution of E-value in four databases

表3 Unigenes的物種分布統計表（前10位）Table 3 Unigenes distribution in different species (top 10 species)

2.3 Unigenes的Gene Ontology（GO）分類

對所有unigenes進行GO分類，共有9 679個unigenes具有GO功能注釋，這些基因的功能分為生物學過程、細胞組分和分子功能3類。如圖3所示，生物學進程中，注釋到行為、生物黏附、生物調控、細胞殺傷、細胞成分組織或生物合成、細胞進程、生長、免疫系統進程、定位、運動、代謝進程多組織進程、多細胞組織進程、生殖、生殖進程、應激、信號、單一有機體進程的unigenes數目分別為22、92、1 655、2、519、4 156、7、16、1 132、36、4 146、52、220、25、31、819、593和3 263個；細胞組分中，注釋到細胞、細胞連接、細胞零、細胞外基質、細胞外基質組分、胞外區、胞外區零件、大分子復合物、細胞膜、細胞膜零件、細胞膜內腔、細胞器、細胞器零件、突觸、突觸零件、病毒、病毒零件的unigenes數目分別為2 900、15、2 900、39、6、27、26、1 287、1 702、1 511、65、1 893、700、41、37、49和49；分子功能中，注釋到抗氧化活性、結合、催化活性、通道調節子活性、電子轉運活性、酶調節活性、脒基核苷酸交換因子活性、分子功能調節因子、分子轉換器活性、核酸結合轉錄因子活性、蛋白結合轉錄因子活性、結構分子活性、轉運子活性的unigenes數目分別為48、4 935、2、34、89、61、150、316、177、30、402和521。

2.4 Unigenes的KEGG代謝通路注釋

對所有unigenes進行KEGG代謝通路注釋，共有4 517個 unigenes注釋到 KEGG數據庫中，這些unigenes的通路信息如圖4所示。這些unigenes分布于216個已知的代謝通路中，其中富集數量最多的10個代謝通路是核糖體、碳代謝以及內質網蛋白加工、內吞、RNA轉運、嘌呤代謝、氧化磷酸化、剪接體、氨基酸生物合成和泛素介導的蛋白水解（表 4）。此外，注釋到溶酶體、MAPK信號通路、Jak-STAT信號通路、昆蟲激素生物合成、黑化作用、Ras信號通路、凋亡和嗅覺轉化上的unigenes分別為119、27、25、16、10、7、4和4個，其中富集在免疫通路上的unigenes有可能在中蜂幼蟲響應病原微生物入侵的免疫應答過程中發揮關鍵作用。

圖2 Unigenes的KOG功能分類Fig.2 KOG classification of unigenes

2.5 SSR分子標記鑒定

利用MISA軟件從43 557條 unigenes中共鑒定出13 448個SSR位點。其中二核苷酸重復最多，數目達到 7 804（58.03%），其次依次為三核苷酸、四核苷酸、五核苷酸和六核苷酸重復，數目分別為 3 797（28.23%）、1 307（9.72%）、339（2.52%）和201（1.49%）（表5）。通過對SSR基元進行分析，發現AT/AT出現的頻率最高（30.4%），其次為AG/CT（22%），不同類型的SSR在總SSR中所占的比例如圖5所示。

在上述的 13 448個 SSR位點中，利用 Primer Primer 5軟件在隨機挑選的20個SSRs序列兩側設計特異性引物，引物序列信息如附表2所示。提取4、5、6日齡中蜂幼蟲腸道總DNA，等摩爾混合作為模板進行PCR擴增。

表4 注釋到KEGG數據庫前10位代謝通路Table 4 Top 10 pathways of unigenes annotated in KEGG pathway database

表5 中蜂幼蟲腸道SSR位點統計Table 5 Characteristics of SSRs in A. c. cerana larval gut

圖3 Unigenes的GO分類Fig.3 GO classification of all unigenes

圖4 Unigenes的KEGG代謝通路注釋Fig.4 KEGG pathway annotation of all unigenes

圖5 不同串聯重復單元類型的SSR在總SSR中所占比例Fig.5 Frequency of SSR motif in total SSRs

PCR產物經1%瓊脂糖凝膠電泳檢測，結果顯示，有6對SSRs特異性引物（SSR9、SSR11、SSR14、 SSR16、SSR19、SSR20）對3個不同來源的中蜂幼蟲腸道樣品都擴增出了具有多態性的特異性條帶（圖6），說明這些SSR位點有望作為中蜂幼蟲特有的分子標記，基于轉錄組數據大規模開發SSR分子標記具有良好的前景。

3 討論

圖6 國內3個來源中蜂幼蟲腸道SSR位點鑒定Fig.6 SSR loci identification of A. c. cerana larval gut samples from three different regions in China

2015年，韓國的研究人員公布了東方蜜蜂雄蜂的基因組信息，但當時并沒有公布基因的位置及功能注釋信息[22]。WANG等[16]曾對中蜂進行過轉錄組測序，因測序組織包括 3日齡工蜂幼蟲、1日齡工蜂蛹、1日齡成年工蜂、采集蜂以及哺育蜂，故該轉錄組信息較為復雜、不夠專一。腸道是中蜂幼蟲的主要免疫器官，在抵御病原微生物入侵過程中扮演著重要角色。本研究利用RNA seq技術對中蜂腸道進行深度測序，成功組裝并注釋了專一的中蜂幼蟲腸道參考轉錄組，將有力推動中蜂及其幼蟲的分子及組學研究，如中蜂幼蟲響應球囊菌或東方蜜蜂微孢子蟲（Nosema ceranae）侵染過程中的免疫應答及分子調控研究。

養蜂生產中，意蜂幼蟲易被球囊菌感染而罹患白堊病[28]，而中蜂幼蟲具有較強的球囊菌抗性，但偶爾可見患病幼蟲。通常認為中蜂具有較強的清理行為，表現出更強的群體防御[7]，但中蜂幼蟲個體水平的免疫防御卻鮮有研究，其在中蜂幼蟲球囊菌抗性方面所發揮的作用值得深入探討。未來筆者課題組將在本研究組裝并注釋的參考轉錄組的基礎上，對病原脅迫過程中中蜂幼蟲的病原-宿主互作機制、免疫應答機制及分子調控機制進行深入系統的研究。

SSR分子標記的傳統開發方法是通過構建 DNA文庫進行篩選，成本高且效率低，而高通量測序技術的應用為大規模篩選 SSR分子標記帶來曙光[15]。目前，已報道的中蜂SSR分子標記非常少[24-25]，嚴重阻礙中蜂的品種鑒定及遺傳進化等研究。本研究基于中蜂幼蟲腸道的轉錄組數據預測潛在的SSR分子標記，隨機選取的20對特異性SSR引物中有6對可在北京、遼寧興城和四川成都3個不同來源的中蜂幼蟲樣品中擴增出具有多態性的片段，這些新開發的SSR分子標記有助于中蜂的基因圖譜構建、基因多樣性分析、基因定位等[29-30]研究的深入開展，說明基于轉錄組測序數據大規模開發SSR分子標記具有良好的應用前景。

4 結論

成功組裝中蜂幼蟲腸道參考轉錄組并對其進行了功能及代謝通路注釋，可為中蜂幼蟲的分子及組學研究提供重要的參考信息，也可用于補充、豐富和檢驗已公布的東方蜜蜂基因組，基于該轉錄組數據開發出的6個中蜂的SSR分子標記可應用于中蜂的基因圖譜構建、基因多樣性分析、基因定位等研究，同時也說明利用轉錄組數據開發非模式生物SSRs的方法可行。

[1] ZAYED A, ROBINSON G E. Understanding the relationship between brain gene expression and social behavior: lessons from the honey bee. Annual Review of Genetics, 2012, 46: 591-615.

[2] BEGNA D, HAN B, FENG M, FANG Y, LI J K. Differential expressions of nuclear proteomes between honeybee (Apis mellifera L.) queen and worker larvae: a deep insight into caste pathway decisions. Journal of Proteome Research, 2012, 11(2): 1317-1329.

[3] FORET S, KUCHARSKI R, PELLEGRINI M, FENG S, JACOBSEN S E, ROBINSON G E, MALESZKA R. DNA methylation dynamics, metabolic fluxes, gene splicing, and alternative phenotypes in honey bees. Proceedings of the National Academy of Sciences of the United States of America, 2012, 109(13): 4968-4973.

[4] GALIZIA G E D, GIU M. Honeybee Neurobiology and Behavior: A Tribute to Randolf Menzel. Springer, 2012: 521.

[5] WEINSTOCK G M, ROBINSON G E, GIBBS R A, WEINSTOCK G M. Insights into social insects from the genome of the honeybee Apis mellifera. Nature, 2006, 443(7114): 931-949.

[6] PENG Y S, NASR M E, LOCKE S J. Geographical races of Apis cerana Fabricius in China and their distribution. Review of recent Chinese publications and a preliminary statistical analysis. Apidologie, 1989, 20(1): 9-20.

[7] LIN Z, PAGE P, LI L, QIN Y, ZHANG Y, HU F, NEUMANN P, ZHENG H, DIETEMANN V. Go east for better honey bee health: Apis cerana is faster at hygienic behavior than A. mellifera. PLoS ONE, 2016, 11(9): e0162647.

[8] FRIES I, WEI H, SHI W, CHEN S J. Grooming behavior and damaged mites (Varroa jacobsoni) in Apis cerana cerana and Apis mellifera ligustica. Apidologie, 1996, 27(1): 3-11.

[9] 周冰峰, 朱翔杰. 論中華蜜蜂種質資源的保護//中國蜂業科技可持續發展學術研討會暨蜂業科技與生態研討會, 2008: 110-117.

ZHOU B F, ZHU X J. The protection of Apis cerana cerana//China Apicultural Conference on Sustainable Development & Apicultural Technology and Ecology, 2008: 110-117. (in Chinese)

[10] RATH W. Co-adaptation of Apis cerana Fabr. and Varroa jacobsoni Oud. Apidologie, 1999, 30(2/3): 97-110.

[11] PENG Y S, FANG Y, XU S, GE L. The resistance mechanism of the Asian honey bee, Apis cerana Fabr. to an ectoparasitic mite, Varroa jacobsoni Oudemans. Journal of Invertebrate Pathology, 1987, 49(1): 54-60.

[12] 龔一飛, 張其康. 蜜蜂分類與進化. 福州: 福建科學技術出版社, 2000: 21-26.

GONG Y F, ZHANG Q K. Classification and Evolution of Apis. Fuzhou: Fujian Science and Technology Press, 2000: 21-26. (in Chinese)

[13] 李紅亮. 中華蜜蜂頭部ESTs文庫構建和主要觸角特異蛋白基因克隆、定位及其表達鑒定[D]. 杭州: 浙江大學, 2007.

LI H L. Construction of cDNA libraries from head of bees, cloning, expression and subcellular localization of the antenna specific protein gene in the Chinese honeybee, Apis cerana cerana[D]. Hangzhou: Zhejiang University, 2007. (in Chinese)

[14] NOVAES E, DROST D R, FARMERIE W G, PAPPAS G J, GRATTAPAGLIA D, SEDEROFF R R, KIRST M. High-throughput gene and SNP discovery in Eucalyptus grandis, an uncharacterized genome. BMC Genomics, 2008, 9(1): 312.

[15] GUO R, WANG S, XUE R, CAO G, HU X, HUANG M, ZHANG Y, LU Y, ZHU L, CHEN F, LIANG Z, KUANG S, GONG C. The gene expression profile of resistant and susceptible Bombyx mori strains reveals cypovirus-associated variations in host gene transcript levels. Applied Microbiology and Biotechnology, 2015, 99(12): 5175-5187.

[16] WANG Z, GERSTEIN M, SNYDER M. RNA-seq: A revolutionary tool for transcriptomics. Nature Reviews Genetics, 2009, 10(1): 57-63.

[17] HAAS B J, PAPANICOLAOU A, YASSOUR M, GRABHERR M, BLOOD P D, BOWDEN J, COUGER M B, ECCLES D, LI B, LIEBER M. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nature Protocols, 2013, 8(8): 1494-1512.

[18] VAN DIJK E L, AUGER H, JASZCZYSZYN Y, THERMES C. Ten years of next-generation sequencing technology. Trends in Genetics, 2014, 30(9): 418-426.

[19] TRAPNELL C, WILLIAMS B A, PERTEA G, MORTAZAVI A, KWAN G, VAN BAREN M J, SALZBERG S L, WOLD B J, PACHTER L. Transcript assembly and quantification by RNA-seq reveals unannotated transcripts and isoform switching during celldifferentiation. Nature Biotechnology, 2010, 28(5): 511-515.

[20] WANG Z L, LIU T T, HUANG Z Y, WU X B, YAN W Y, ZENG Z J. Transcriptome analysis of the Asian honey bee Apis cerana cerana. PLoS ONE, 2012, 7(10): e47954.

[21] CORNMAN R S, BENNETT A K, MURRAY K D, EVANS J D, ELSIK C G, ARONSTEIN K. Transcriptome analysis of the honey bee fungal pathogen, Ascosphaera apis: implications for host pathogenesis. BMC Genomics, 2012, 13: 285.

[22] PARK D, JUNG J W, CHOI B S, JAYAKODI M, LEE J, LIM J, YU Y, CHOI Y S, LEE M L, PARK Y. Uncovering the novel characteristics of Asian honey bee, Apis cerana, by whole genome sequencing. BMC Genomics, 2015, 16: 1.

[23] 梁勤, 張立卿. 利用微衛星標記分析福建4個中華蜜蜂群體的遺傳多樣性. 福建農林大學學報 (自然科學版), 2009, 38(4): 388-392.

LIANG Q, ZHANG L Q. Analysis of genetic diversity of four Apis cerana cerana populations in Fujian Province with microsatellite markers. Journal of Fujian Agriculture and Forestry University (Natural Science Edition), 2009, 38(4): 388-392. (in Chinese)

[24] 徐新建, 周姝婧, 朱翔杰, 周冰峰. 海南島中華蜜蜂遺傳多樣性的微衛星DNA分析. 昆蟲學報, 2013, 56(5): 554-560.

XU X J, ZHOU S J, ZHU X J, ZHOU B F. Microsatellite DNA analysis of genetic diversity of Apis cerana cerana in Hainan Island, southern China. Acta Entomologica Sinica, 2013, 56(5): 554-560. (in Chinese)

[25] 于瀛龍, 周姝婧, 徐新建, 朱翔杰, 周冰峰. 長白山中華蜜蜂(Apis cerana cerana)遺傳多樣性分析. 福建農林大學學報 (自然科學版), 2013, 42(6): 643-647.

YU Y L, ZHOU S J, XU X J, ZHU X J, ZHOU B F. Analysis on genetic diversity of Apis cerana cerana in Changbai Mountains. Journal of Fujian Agriculture and Forestry University (Natural Science Edition), 2013, 42(6): 643-647. (in Chinese)

[26] 王倩, 孫亮先, 肖培新, 劉鋒, 康明江, 胥保華. 室內人工培育中華蜜蜂幼蟲技術研究. 山東農業科學, 2009(11): 113-116.

WANG Q, SUN L X, XIAO P X, LIU F, KANG M J, XU B H. Study on technology for indoor artificial feeding of Apis cerana cerana larvae. Shandong Agricultural Sciences, 2009(11): 113-116. (in Chinese)

[27] GRABHERR M G, HAAS B J, YASSOUR M, LEVIN J Z, THOMPSON D A, AMIT I, XIAN A, FAN L, RAYCHOWDHURY R, ZENG Q. Trinity: reconstructing a full-length transcriptome without a genome from RNA-seq data. Nature Biotechnology, 2011, 29(7): 644-652.

[28] GUPTA R K, REYBROECK W. Honeybee Pathogens and Their Management, 2014: Springer, Netherlands.

[29] SOMRIDHIVEJ B, WANG S, SHA Z, LIU H, QUILANG J, XU P, LI P, HU Z, LIU Z. Characterization, polymorphism assessment, and database construction for microsatellites from bac end sequences of channel catfish (Ictalurus punctatus): A resource for integration of linkage and physical maps. Aquaculture, 2008, 275(1/4): 76-80.

[30] GAO X, HAN J, LU Z, LI Y, HE C. Retracted: characterization of the spotted seal phoca largha transcriptome using Illumina paired-end sequencing and development of SSR markers. Comparative Biochemistry and Physiology-Part D: Genomics and Proteomics, 2012, 7(3): 277-284.

（責任編輯 岳梅）

De novo Transcriptome Assembly for Apis cerana cerana Larval Gut and Identification of SSR Molecular Markers

XU XiJian1, GUO Rui1, LUO Qun2, XIONG CuiLing1, LIANG Qin1, ZHANG ChuanLian2, ZHENG YanZhen1, ZHANG ZhaoNan1, HUANG ZhiJian1, ZHANG Lu1, LI WenDong1, CHEN DaFu1

(1College of Bee Science, Fujian Agriculture and Forestry University, Fuzhou 350002;2Apiculture Institute of Jiangxi Province, Nanchang 330201)

Abstract:【Objective】 The objective of this study is to de novo assemble a reference transcriptome for Apis cerana cerana larval gut, perform gene function and pathway annotation for this transcriptome, and to identify specific SSR molecular markers for A. c. cerana larvae. 【Method】 3-day-old instar A. c. cerana larvae were fed with the purified Ascosphaera apis spores, the guts of 4-, 5- or 6-day-old honeybee larvae were sampled and used as sequencing material for RNA seq. After filtration, clean reads were obtained, and unigenes were assembled using Trinity software. BLASTX tool (E-value＜10?5) was used to search the unigenes against NCBI Nr, Swiss-Prot protein, KOG as well as KEGG databases to perform gene function and pathway annotation. MISA software was used to search microsatellite markers in the larval gut’s transcriptome. The specific primers of all SSRs were designed using Primer Premier 5 program and several pairs were used to amplify SSR loci in A. c. cerana larvae samples from 3 different regions (Beijing, Xingcheng, and Chengdu) in China by method of PCR. 【Result】 In this study, RNA seq produced 35 670 000 high quality reads, which were assembled into 43 557 unigenes with a mean length of 898 nt. 18 225 unigenes were annotated in the public protein databases. A total of 11 442 unigenes had a KOG classification and they distributed in 25 KOG categories, among them, RNA processing and modification was the largest group (1 249). 9 679 unigenes could be classified into three gene ontology (GO) categories, in which the mostly enriched ones were cellular process (4 201 unigenes), cell (2 900 unigenes) and binding (4 935 unigenes). 4 517 unigenes were annotated to 216 KEGG pathways, among them, ribosome (385 unigenes) was the largest. Finally, 13 448 SSRs were found in 7 763 unigenes, and 6 out 20 SSR loci could be successfully amplified in A. c. cerana larvae samples from 3 different regions in China using PCR. 【Conclusion】 This study assembled and annotated a reference transcriptome for A. c. cerana larval gut, which will provide a key information not only to studies on eastern honeybee and its larvae such as molecular biology and omics, but also to improve and validate the genome of A. cerana. SSR markers developed here could be applied to future investigation of A. c. cerana including gene map construction, genetic diversity analysis as well as gene location. Meanwhile, this study suggested that developing molecular markers using transcriptome data of non-model organism is a rapid and efficient method.

RNA seq; de novo assembly; Apis cerana cerana; unigenes; SSR

2016-09-29；接受日期：2016-12-12

國家現代農業產業技術體系（蜜蜂）建設專項（CARS-45-KXJ7）、福建省自然科學基金（2013J01070）

聯系方式：徐細建，Tel：18020870542；E-mail：xxjlhj2006@163.com。郭睿，Tel：15205080780；E-mail：fafu_ruiguo@126.com。徐細建和郭睿為同等貢獻作者。通信作者陳大福，Tel：0591-83726835；E-mail：dfchen826@163.com