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

（1福建農(nóng)林大學(xué)蜂學(xué)學(xué)院，福州 350002；2江西省養(yǎng)蜂研究所，南昌 330201）

中華蜜蜂幼蟲腸道參考轉(zhuǎn)錄組的de novo組裝及SSR分子標(biāo)記鑒定

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

（1福建農(nóng)林大學(xué)蜂學(xué)學(xué)院，福州 350002；2江西省養(yǎng)蜂研究所，南昌 330201）

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

RNA seq；參考轉(zhuǎn)錄組；中華蜜蜂；unigene；SSR

0 引言

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

1 材料與方法

試驗于2015年12月至2016年8月在福建農(nóng)林大學(xué)蜂學(xué)學(xué)院蜜蜂保護學(xué)實驗室完成。

1.1 供試材料

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

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產(chǎn)品；cDNA 文庫構(gòu)建試劑盒TruSeqTMDNA Sample Prep Kit -Set A為美國Illumina公司產(chǎn)品。其他試劑均為國產(chǎn)分析純。

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

1.3 方法

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

1.3.3 cDNA文庫構(gòu)建及 RNA seq 利用 RNAiso Reagent試劑盒抽提中蜂幼蟲腸道組織的總RNA，然后用RNase-free DNaseI去除基因組DNA殘留。RNA的質(zhì)量通過瓊脂糖凝膠電泳和 NanoDrop ND-2000（NanoDrop，Wilmington，DE，USA）進行檢測。利用Oligotex mRNA Kits Midi試劑盒說明書，純化各樣品總RNA中的mRNA。以10 μg mRNA作為模板，GsuI-oligo dT作為反轉(zhuǎn)錄引物，用1 000 U Superscript II reverse transcriptase在42℃下孵育1 h合成第1鏈cDNA；隨后利用高碘酸鈉氧化mRNA的5′端帽子結(jié)構(gòu)，并連接生物素；通過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進行構(gòu)建和TruSeq PE Cluster Kit進行擴增。委托廣州基迪奧生物技術(shù)有限公司對上述12個腸道樣品進行深度測序，測序平臺為Illumina Hiseq 2500，各腸道樣品的3個生物學(xué)重復(fù)均同時進行測序。

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

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

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

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

2 結(jié)果

2.1 中蜂幼蟲腸道的RNA seq及參考轉(zhuǎn)錄組de novo組裝

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

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

2.2 Unigenes注釋

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

表1 中蜂幼蟲腸道參考轉(zhuǎn)錄組組裝結(jié)果統(tǒng)計Table 1 Summary of A. c. cerana larval gut’s reference transcriptome assembled in this study

表2 公共蛋白數(shù)據(jù)庫注釋統(tǒng)計表Table 2 Summary of annotation information of all unigenes in public protein databases

注釋到Nr數(shù)據(jù)庫中unigenes的E-value分布顯示（圖1），比對到物種序列的E-value均＜10-5，其中E-value＜10-100的有49.76%，說明比對結(jié)果可信度較高。注釋基因同源序列的物種分布統(tǒng)計結(jié)果顯示前 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的基因數(shù)為5 753（31.57%），注釋到A. dorsata和A. florea的基因數(shù)分別為3 695（20.27%）和2 489（13.66%）（表3）。

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

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

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

2.3 Unigenes的Gene Ontology（GO）分類

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

2.4 Unigenes的KEGG代謝通路注釋

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

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

2.5 SSR分子標(biāo)記鑒定

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

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

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

表5 中蜂幼蟲腸道SSR位點統(tǒng)計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 不同串聯(lián)重復(fù)單元類型的SSR在總SSR中所占比例Fig.5 Frequency of SSR motif in total SSRs

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

3 討論

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

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

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

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

4 結(jié)論

成功組裝中蜂幼蟲腸道參考轉(zhuǎn)錄組并對其進行了功能及代謝通路注釋，可為中蜂幼蟲的分子及組學(xué)研究提供重要的參考信息，也可用于補充、豐富和檢驗已公布的東方蜜蜂基因組，基于該轉(zhuǎn)錄組數(shù)據(jù)開發(fā)出的6個中蜂的SSR分子標(biāo)記可應(yīng)用于中蜂的基因圖譜構(gòu)建、基因多樣性分析、基因定位等研究，同時也說明利用轉(zhuǎn)錄組數(shù)據(jù)開發(fā)非模式生物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] 周冰峰, 朱翔杰. 論中華蜜蜂種質(zhì)資源的保護//中國蜂業(yè)科技可持續(xù)發(fā)展學(xué)術(shù)研討會暨蜂業(yè)科技與生態(tài)研討會, 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] 龔一飛, 張其康. 蜜蜂分類與進化. 福州: 福建科學(xué)技術(shù)出版社, 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文庫構(gòu)建和主要觸角特異蛋白基因克隆、定位及其表達鑒定[D]. 杭州: 浙江大學(xué), 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] 梁勤, 張立卿. 利用微衛(wèi)星標(biāo)記分析福建4個中華蜜蜂群體的遺傳多樣性. 福建農(nóng)林大學(xué)學(xué)報 (自然科學(xué)版), 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] 徐新建, 周姝婧, 朱翔杰, 周冰峰. 海南島中華蜜蜂遺傳多樣性的微衛(wèi)星DNA分析. 昆蟲學(xué)報, 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)遺傳多樣性分析. 福建農(nóng)林大學(xué)學(xué)報 (自然科學(xué)版), 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] 王倩, 孫亮先, 肖培新, 劉鋒, 康明江, 胥保華. 室內(nèi)人工培育中華蜜蜂幼蟲技術(shù)研究. 山東農(nóng)業(yè)科學(xué), 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.

（責(zé)任編輯 岳梅）

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

國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系（蜜蜂）建設(shè)專項（CARS-45-KXJ7）、福建省自然科學(xué)基金（2013J01070）

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