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

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

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

0 引言

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

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

1.1 供試材料

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

1.2 主要試劑及儀器

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

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

1.3 方法

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

1.3.3 cDNA文庫(kù)構(gòu)建及 RNA seq 利用 RNAiso Reagent試劑盒抽提中蜂幼蟲腸道組織的總RNA，然后用RNase-free DNaseI去除基因組DNA殘留。RNA的質(zhì)量通過瓊脂糖凝膠電泳和 NanoDrop ND-2000（NanoDrop，Wilmington，DE，USA）進(jìn)行檢測(cè)。利用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對(duì)上述cDNA進(jìn)行純化，cDNA文庫(kù)通過TruSeqTMDNA Sample Prep Kit-Set A進(jìn)行構(gòu)建和TruSeq PE Cluster Kit進(jìn)行擴(kuò)增。委托廣州基迪奧生物技術(shù)有限公司對(duì)上述12個(gè)腸道樣品進(jìn)行深度測(cè)序，測(cè)序平臺(tái)為Illumina Hiseq 2500，各腸道樣品的3個(gè)生物學(xué)重復(fù)均同時(shí)進(jìn)行測(cè)序。

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

1.3.5 Unigenes注釋 利用 BLASTx（E-value＜10?5）將測(cè)序序列比對(duì)NCBI nr數(shù)據(jù)庫(kù)（http://www.ncbi. nlm.nih.gov）、Swiss-Prot 數(shù)據(jù)庫(kù)（http://www.expasy.ch/ sprot）、KOG（Clusters of orthologous groups for eukaryotic complete genomes）數(shù)據(jù)庫(kù)（ftp://ftp.ncbi.nih. gov/pub/COG/KOG/kyva）和KEGG代謝通路（pathway）數(shù)據(jù)庫(kù)（http://www.genome.jp/kegg/）。利用BLASTX將組裝出來的unigenes序列與Nr數(shù)據(jù)庫(kù)進(jìn)行比對(duì)后，取每個(gè)unigenes在Nr庫(kù)中比對(duì)結(jié)果最好（E值最低）的那一條序列為對(duì)應(yīng)同源序列（如有并列，取第一條）確定同源序列所屬物種，統(tǒng)計(jì)比對(duì)到各個(gè)物種的同源序列數(shù)量?；贜r database注釋結(jié)果，利用Blast2GO進(jìn)行unigenes的GO注釋，利用WEGO軟件對(duì)每一個(gè)轉(zhuǎn)錄本進(jì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位點(diǎn)：二核苷酸重復(fù)≥6次，三核苷酸重復(fù)≥5次，四核苷酸重復(fù)≥5次，五核苷酸重復(fù)≥5次和六核苷酸重復(fù)≥5次。根據(jù) MISA的輸出結(jié)果，利用 Primer Premier 5（PREMIER Biosofe Int.，Palo Alto，CA）對(duì)每一個(gè)含有16 bp堿基重復(fù)的SSR設(shè)計(jì)引物。

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

2 結(jié)果

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

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

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

2.2 Unigenes注釋

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

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

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

注釋到Nr數(shù)據(jù)庫(kù)中unigenes的E-value分布顯示（圖1），比對(duì)到物種序列的E-value均＜10-5，其中E-value＜10-100的有49.76%，說明比對(duì)結(jié)果可信度較高。注釋基因同源序列的物種分布統(tǒng)計(jì)結(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個(gè)unigenes分布于25個(gè)基因家族（圖2）。其中，注釋基因數(shù)最多的為信號(hào)轉(zhuǎn)導(dǎo)機(jī)制，其次為一般功能預(yù)測(cè)和翻譯后修飾、蛋白翻轉(zhuǎn)和分子伴侶。值得注意的是，有 170條 unigenes注釋到防御機(jī)制，它們可能在中蜂幼蟲抵御病原入侵過程發(fā)揮重要作用。

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

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

2.3 Unigenes的Gene Ontology（GO）分類

對(duì)所有unigenes進(jìn)行GO分類，共有9 679個(gè)unigenes具有GO功能注釋，這些基因的功能分為生物學(xué)過程、細(xì)胞組分和分子功能3類。如圖3所示，生物學(xué)進(jìn)程中，注釋到行為、生物黏附、生物調(diào)控、細(xì)胞殺傷、細(xì)胞成分組織或生物合成、細(xì)胞進(jìn)程、生長(zhǎng)、免疫系統(tǒng)進(jìn)程、定位、運(yùn)動(dòng)、代謝進(jìn)程多組織進(jìn)程、多細(xì)胞組織進(jìn)程、生殖、生殖進(jìn)程、應(yīng)激、信號(hào)、單一有機(jī)體進(jìn)程的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個(gè)；細(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)運(yùn)活性、酶調(diào)節(jié)活性、脒基核苷酸交換因子活性、分子功能調(diào)節(jié)因子、分子轉(zhuǎn)換器活性、核酸結(jié)合轉(zhuǎn)錄因子活性、蛋白結(jié)合轉(zhuǎn)錄因子活性、結(jié)構(gòu)分子活性、轉(zhuǎn)運(yùn)子活性的unigenes數(shù)目分別為48、4 935、2、34、89、61、150、316、177、30、402和521。

2.4 Unigenes的KEGG代謝通路注釋

對(duì)所有unigenes進(jìn)行KEGG代謝通路注釋，共有4 517個(gè) unigenes注釋到 KEGG數(shù)據(jù)庫(kù)中，這些unigenes的通路信息如圖4所示。這些unigenes分布于216個(gè)已知的代謝通路中，其中富集數(shù)量最多的10個(gè)代謝通路是核糖體、碳代謝以及內(nèi)質(zhì)網(wǎng)蛋白加工、內(nèi)吞、RNA轉(zhuǎn)運(yùn)、嘌呤代謝、氧化磷酸化、剪接體、氨基酸生物合成和泛素介導(dǎo)的蛋白水解（表 4）。此外，注釋到溶酶體、MAPK信號(hào)通路、Jak-STAT信號(hào)通路、昆蟲激素生物合成、黑化作用、Ras信號(hào)通路、凋亡和嗅覺轉(zhuǎn)化上的unigenes分別為119、27、25、16、10、7、4和4個(gè)，其中富集在免疫通路上的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個(gè)SSR位點(diǎn)。其中二核苷酸重復(fù)最多，數(shù)目達(dá)到 7 804（58.03%），其次依次為三核苷酸、四核苷酸、五核苷酸和六核苷酸重復(fù)，數(shù)目分別為 3 797（28.23%）、1 307（9.72%）、339（2.52%）和201（1.49%）（表5）。通過對(duì)SSR基元進(jìn)行分析，發(fā)現(xiàn)AT/AT出現(xiàn)的頻率最高（30.4%），其次為AG/CT（22%），不同類型的SSR在總SSR中所占的比例如圖5所示。

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

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

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

3 討論

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

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

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

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

4 結(jié)論

成功組裝中蜂幼蟲腸道參考轉(zhuǎn)錄組并對(duì)其進(jìn)行了功能及代謝通路注釋，可為中蜂幼蟲的分子及組學(xué)研究提供重要的參考信息，也可用于補(bǔ)充、豐富和檢驗(yàn)已公布的東方蜜蜂基因組，基于該轉(zhuǎn)錄組數(shù)據(jù)開發(fā)出的6個(gè)中蜂的SSR分子標(biāo)記可應(yīng)用于中蜂的基因圖譜構(gòu)建、基因多樣性分析、基因定位等研究，同時(shí)也說明利用轉(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ì)資源的保護(hù)//中國(guó)蜂業(yè)科技可持續(xù)發(fā)展學(xué)術(shù)研討會(huì)暨蜂業(yè)科技與生態(tài)研討會(huì), 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] 龔一飛, 張其康. 蜜蜂分類與進(jìn)化. 福州: 福建科學(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文庫(kù)構(gòu)建和主要觸角特異蛋白基因克隆、定位及其表達(dá)鑒定[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個(gè)中華蜜蜂群體的遺傳多樣性. 福建農(nóng)林大學(xué)學(xué)報(bào) (自然科學(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é)報(bào), 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] 于瀛龍, 周姝婧, 徐新建, 朱翔杰, 周冰峰. 長(zhǎng)白山中華蜜蜂(Apis cerana cerana)遺傳多樣性分析. 福建農(nóng)林大學(xué)學(xué)報(bào) (自然科學(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

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

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