玫烟色棒束孢转录组测序及潜在致病相关基因分析

谢梅琼,王龙江,何余容,吕利华

(1.宜春学院 生命科学与资源环境学院,江西 宜春 336000; 2.宜春学院 化学与生物工程学院,江西 宜春 336000; 3.华南农业大学 植物保护学院, 广东 广州 510642; 4.广东省农业科学院 植物保护研究所,广东 广州 510642)

虫生真菌作为一类重要的微生物杀虫剂,被广泛应用于农林害虫的防治。玫烟色棒束孢(Isariafumosorosea)作为一种有名的杀虫真菌,在国内外害虫防治中已得到广泛认可,其寄主广泛,包括小菜蛾(PlutellaxylostellaLinnaeus)、烟粉虱(BemisababaciGennadius)、棉蚜(AphisgossypiiGlover)等46个科重要的农业和果树害虫[1-3]。目前,有关玫烟色棒束孢的研究主要集中在生物学、生理生化及遗传多样性等方面[4-6],而关于其功能基因的挖掘和研究相对较少。么瑞娜[7]以小菜蛾为诱导物,利用酵母双杂交技术,构建了玫烟色棒束孢酵母双杂交cDNA文库,用于筛选玫烟色棒束孢与小菜蛾相互作用的蛋白及其对应的功能基因。黄宇[8]利用DNA步移法获得了玫烟色棒束孢Ifcdp1基因的DNA和cDNA全长,并利用基因敲除和互补的方法验证了该基因的功能。孟祥云[9]和汤强等[10]分别克隆了玫烟色棒束孢Ifuchi1 几丁质酶基因和Ifuchit2几丁质酶基因,并分析了基因的序列特征。董章勇等[11]从玫烟色棒束孢ARSEF 2679菌株全基因组的10 061个蛋白质序列中筛选到了1 112个具有信号肽的蛋白序列。目前,玫烟色棒束孢的部分功能基因虽已得到克隆、鉴定,但大量相关基因及其功能仍处于未知状态。

随着高通量测序及生物信息学领域的快速发展,人们对真菌的研究和利用范畴进一步拓宽。Gao等[12]利用高通量测序技术,对金龟子绿僵菌(Metarhiziumanisopliae)和蝗绿僵菌(Metarhiziumacridum)进行了全基因组测序,揭示两种绿僵菌的进化史,并通过比对病原菌和寄主互作基因数据库(pathogen-host interaction gene database),发现绿僵菌属真菌约有16%基因与植物病原菌致病相关基因相似,该研究结果为绿僵菌属基因的鉴定及生防菌株的基因改良奠定了基础。Zheng等[13]利用高通量测序技术对蛹虫草(Cordycepsmilitaris)进行了全基因组测序,阐明了蛹虫草的进化史,为蛹虫草的开发利用提供了基础。Xiao等[14]对球孢白僵菌(Beauveriabassiana)进行了全基因组测序,系统分析了其进化史,并发现了该物种特有的毒力基因,为球孢白僵菌的开发利用提供了基因资源。Hao等[15]利用第三代测序技术,对禾谷镰刀菌进行了全基因组的从头测序,并预测了其致病相关的碳水化合物和效应蛋白,为该菌致病基因的鉴定奠定了基础。Shang等[16]对玫烟色棒束孢ARSEF 2679菌株进行了基因组的从头测序,为玫烟色棒束孢基因组学的研究提供了参考。

本研究利用高通量测序技术对玫烟色棒束孢菌株IFCF01进行转录组测序,对获得的数据进行组装和功能注释,同时筛选其潜在的致病相关基因和简单重复序列(simple sequence repeats, SSR)分子标记,为玫烟色棒束孢功能基因的研究及菌株基因改良提供基因资源。

1 材料与方法

1.1 真菌的培养及采集

玫烟色棒束孢菌株IFCF01由华南农业大学植物保护学院生物防治研究室提供,保存于-80 ℃。将菌株接种于铺有玻璃纸的马铃薯葡萄糖琼脂(potato dextrose agar,PDA)培养基平板上,置于(25±1)℃、光照14 h、相对湿度≥90%条件下培养至其产生孢子时(7 d左右),收集菌丝及孢子,经液氮速冻后,保存于-80 ℃中。

1.2 转录组测序及数据组装

将收集好的样品用冰桶送至深圳华大基因公司(BGI)进行RNA提取、文库构建及转录组的测序和序列拼装。采用 Agilent 2100 Bioanalyzer 和 ABI StepOnePlus Real-Time PCR System 进行质量检测,采用Illumina Hiseq2000系统进行转录组测序。对测序获得的原始序列(raw data)进行过滤,去掉包含接头、N比例大于5%及重复和低质量(Q≤10的碱基数占整个片段20%以上)的片段后,获得所需片段,用于后续的数据分析。使用Trinity软件[17]对所需片段进行从头组装,将具有一定长度重叠的序列组合成更长的连续序列,然后与所需片段重新比对,通过序列方向测定确定来自同一转录本的不同的重复序列以及这些序列之间的距离,将来自同一转录本的序列连在一起,得到两端不能再延长的序列,使用Tgicl软件对其进行去冗余后进一步拼接,再进行同源转录本聚类,得到最终的序列,并将其分为两部分:一部分为单基因簇聚类,即cluster,同一个cluster里面有若干条相似度大于70%的单基因(unigene)(以CL开头,后接基因家族的编号);其余的为单基因簇,即singletons(以unigene开头),代表单独的unigene。

1.3 玫烟色棒束孢转录组的功能注释

将玫烟色棒束孢转录组数据比对蛋白数据库non-redundant(Nr)、Swiss-prot、KEGG(Kyoto Encyclopedia of Genes and Genomes)和COG(E值<1×10-5),并通过核苷酸序列比对(blastn)将转录组的unigene比对核酸数据库Nt(E值<1×10-5)。通过比对获得跟unigene具有最高相似性的蛋白及unigene的蛋白功能注释信息。根据Nr注释信息,采用Blast2GO软件对所得unigene进行Gene Ontology(GO)注释,利用WEGO软件对所有unigene进行GO功能分类。根据KEGG注释信息获得unigene的通路(pathway)注释结果。

1.4 玫烟色棒束孢致病相关基因的挖掘

通过将unigene与病原真菌-寄主相互作用数据库(www.phi-base.org)进行比对分析[18],筛选玫烟色棒束孢潜在致病相关基因(E值<1×10-5)。

1.5 转录组SSR位点分析

使用SSR软件的MicroSAtellite(MISA)对玫烟色棒束孢的所有unigene序列进行SSR位点搜索。

2 结果与分析

2.1 玫烟色棒束孢转录组序列特征

玫烟色棒束孢转录组测序后共获得57 913 442条原始读数(raw reads),过滤后共获得540 896 444条干净读数(clean reads),总长度为4 868 067 960 nt,Q20值达96.35%,测序质量较高。初步组装后,共获得 44 031条重叠群(contig),总长度为26 084 190 bp,平均长度为592 bp,N50 值为1782。再次聚类和组装后,共获得 36 466 条unigene,总长度为34 551 468 bp,平均长度高达947 bp,N50值(用来衡量基因组的拼接质量)为1 676,优于前人报道[19],说明序列组装理想,可满足后续分析需求。contig和unigene序列长度在100～500 bp的最多,分别为32 556条和16 840条;在>1 500～2 000 bp的数量最少,分别为1 642和2 946条(图1)。玫烟色棒束孢原始序列已上传到NCBI数据库(检索号为SRR13236723)。

图1 玫烟色棒束孢转录组中重叠群contigs和单基因unigene序列长度分布Fig.1 Length distributions of contigs and unigene in Isaria fumosorosea transcriptome

2.2 玫烟色棒束孢转录组序列注释结果

通过与Nr、Nt、SwissProt、KEGG、COG和GO这6大数据库比对,共有22 741个unigene得到注释,占所有unigene 的62.36%,有接近40%的unigene未得到注释。注释到Nr、SwissProt、KEGG和COG数据库的unigene数量分别为21 737、14 243、14 364和10 888,此外,有2 274个unigene同时被这4个数据库注释,但大多数unigene仅被部分数据库注释或被单独注释(图2、表1)。

表1 玫烟色棒束孢转录组注释情况汇总表

图2 韦恩图解分析注释到Nr、KEGG、SwissPort和COG数据库的基因分布情况Fig.2 Venn diagram of gene distribution in Nr, KEGG, SwissPort and COG databases

2.3 玫烟色棒束孢转录组Nr功能注释

通过blastx对36 466条unigene进行Nr数据库注释(E值<1×10-5)。结果表明,59.61%的unigene(21 737条)在Nr中至少匹配上一个已知序列,40.39%的序列(9 185条)未被注释(图3-A)。在Nr中,E值为0的基因占比最多,为26.0%;其次为E值在>0～1×10-100的基因,占19.5%;其余E值>1×10-100～<1×10-5的合计占54.5%,说明玫烟色棒束孢转录组较大部分序列与Nr数据库序列具备非常高的同源性(图3-B)。对注释上Nr库的序列进行相似度分析发现,相似度在80%～95%的占比最大,为44.64%,其次为60%～<80%的序列,占30.43%(图3-C)。Nr功能注释到球孢白僵菌B.bassianaARSEF 2860的序列最多,为44.3%;其次为蛹虫草C.militarisCM01,占41.5%,注释到金龟子绿僵菌M.anisopliaeARSEF的占1%,另有13.2% 的基因比对上其他物种,包括蝗绿僵菌(M.acridum)CQMa 102、绿木霉TrichodermavirensGv29-8、尖孢镰刀菌FusariumoxysporumFo5176等。总的来看,在Nr功能注释中,玫烟色棒束孢很大一部分序列主要注释到球孢白僵菌和蛹虫草上(图3-D)。

A,已知序列与未知序列的分布;B, E值在1×10-5以下序列的blastx的比对结果分布;C,同源序列相似度分布;D,注释到Nr库的unigene所属物种分布。A, Distribution of known and unknown sequences; B, Distribution of blast hits for each unique sequence with a cut-off E-value below 1×10-5; C, Distribution of homologous sequences similarities; D, Species distribution of unigene annotated in the Nr database.图3 转录组序列比对Nr数据库的特征Fig.3 Characteristics of transcriptome sequence alignment in Nr database

2.4 玫烟色棒束孢转录组COG功能注释

本研究中,共有10 888个(29.86%) unigene注释到COGs数据库的25个类别,unigene富集程度最高的是主要功能预测(general function predicted only),匹配上4 060个unigene,其次为未知功能(function unknown),有2 801个unigene比对上该模块,而RNA加工修饰(RNA processing and modification) 匹配上的序列最少,只有59个unigene(图4)。这说明,玫烟色棒束孢中仍有一部分基因属于预测功能基因或不能被分类,推测可能跟微生物功能基因数据量相对较少有关,具体原因有待进一步分析。

图4 玫烟色棒束孢unigene的COG功能分类图Fig.4 COG functional classification of Isaria fumosorosea unigenes

2.5 玫烟色棒束孢转录组GO的功能注释

玫烟色棒束孢转录组共有9 617个(26.37%)unigene注释到GO数据库的3个功能模块中,包括生物过程、细胞组分和分子功能。生物过程含21个亚类, 其中的代谢过程(metabolic process)占有的unigene数量最多,为5 237条,其次为细胞过程cellular process,为5 109条,而免疫系统过程(immune system process)及生物附着(biological adhesion)匹配上的序列则比较少,分别为2条和10条。细胞组分包含11个亚类,其中细胞部分(cell part)和细胞的unigene最多,均为2 863条,细胞连接(cell junction)具有最少的匹配序列,仅为1条。分子功能包含14个亚类,其中的催化活性(catalytic activity)和结合(binding) 功能模块属于优势类别,分别有5 322条和4 343条unigene匹配上,蛋白质标签(protein tag)和金属伴侣活性(metallochaperone activity)匹配上的序列非常少,分别为2条和4条(图5)。

2.6 玫烟色棒束孢转录组KEGG代谢通路分析

对转录组的unigene参与的代谢通路进行分析,发现共有14 364个(39.39%) unigene被注释到108个KEGG pathway上。对KEGG通路中unigene富集最多的前20个通路进行分析发现,注释到代谢途径的unigene最多,为4 584个,占总数的32%;其次为次生代谢产物生物合成途径,含有1 897个unigene,占总数的13%;参与淀粉和蔗糖代谢的unigene有1 727个,占12.02%;参与氨基糖和核苷酸糖代谢的unigene有833个,占6%;参与RNA转运的unigene为829个,占6%。综上,玫烟色棒束孢代谢活动比较活跃(表2)。

表2 玫烟色棒束孢KEGG代谢通路中unigene富集最多的前20个通路

2.7 玫烟色棒束孢潜在致病相关基因的筛选

将玫烟色棒束孢36 466条ungigene与病原菌-寄主互作数据库进行blastx比对,分析发现86条unigene与数据库的基因具有较高相似性(E值<1×10-5)。通过COGs数据库注释分析,86条unigene被注释到多个功能模块中,主要包括ABC-类的多种药物运输系统、腺苷三磷酸酶和透性酶组分24%,预测的多酮氧化酶类21%,RNA结合蛋白14%,RTX毒素和钙结合蛋白13%,醌还原酶和锌依赖的氧化还原酶9%,转氨酶和酰胺酶系统5%,染色质修饰相关的转录因子2%,参与大肠杆菌素吸收的膜蛋白1%,翻译抑制因子1%及部分未知功能的基因(7%)(表3)。上述基因很可能参与调控玫烟色棒束孢侵染寄主的过程,并与其本身的各种适应机制相关,如转运蛋白吸收胞外营养、代谢途径调整细胞生理活动、分泌毒素对抗寄主免疫或直接杀死寄主,碳水化合物运输提供营养物质等,基因的具体功能有待进一步研究。

2.8 玫烟色棒束孢转录组SSR位点分析

利用MISA软件对玫烟色棒束孢的36 466条unigene总长为34 551 468 bp进行SSR位点分析,结果共检索到8 572个SSR位点,共有6 629条unigene含有SSR位点,发生频率为18.18%,平均每0.25 kb出现一个SSR位点。有1 409条unigene含有1个以上的SSR位点,化合物形成中存在的SSR数量有466个。

对SSR的重复类型进行统计,结果显示,三核苷酸重复的SSR最多,为4 077条(47.56%),其次为单核苷酸重复的SSR,有2 053条(23.95%),五核苷酸重复的SSR最少,仅为151条(1.76%)。玫烟色棒束孢SSR重复次数在4～24次,以5～7次的重复最多,此区间共有5 300个SSR位点,占61.84%。其中以5次重复的最多,为2 062个SSR位点(24.06%),5次重复中又以三核苷酸重复的SSR数量最多(1890个);其次为6次重复,共有1 925个SSR位点(22.46%),6次重复中同样以三核苷酸重复的SSR数量最多(1 223个);4次重复及8～24次重复的SSR位点共有3 272个,占38.17%。在此类重复中,以单核苷酸重复的SSR最多(2 053个),其次为二核苷酸重复的SSR(823个),四重复的SSR最少,为0(表4)。

表4 玫烟色棒束孢转录组的 SSR类型、数量及分布频率

玫烟色棒束孢的6种SSR重复类型中包含148种重复基序,其中,六核苷酸重复含有的重复基序类型最多,为76种,其次为五核苷酸重复,含38种重复基序,单核苷酸重复含的重复基序最少,仅2种;三、四核苷酸重复类型分别含有10和17种重复基序。在单核苷酸重复基序中,以A/T类型为主,有1 810个,占88.16%,C/G类型有243个,占11.84%;二核苷酸基序以AG/CT和AC/GT为主,分别占55.81%和35.69%;三核苷酸基序以CCG/CGG、AGC/CTG、AAG/CTT和ACG/CGT为主,分别占28.43%、19.48%、13.15%和11.60%;四核苷酸重复基序以AAAG/CTTT(14.85%)、AGGC/CCTG(14.41%)、ACAG/CTGT(13.10%)、ATCC/ATGG (12.66%)和ACCT/AGGT(9.17%)为主;五核苷酸重复基序以AAAAG/CTTTT(22.52%)、AAAAT/ATTTT(13.25%)、AGAGG/CCTCT(9.93%)、AGGGC/CCCTG(7.95%)为主,其他重复基序数量都≤6个;六核苷酸重复基序类型最多(76种),但每一种重复基序类型数量都不多,最多的为ACGGCG/CCGTCG(8.04%),其次为AGCGGC/CCGCTG(4.52%),其他的占比均在4%以下(表5)。

表5 玫烟色棒束孢的SSR重复基序类型及频率

3 结论与讨论

玫烟色棒束抱(Isariafumosorosea)是一类有名的虫生真菌,是生物防治中极具开发应用潜力的种类[20]。目前,关于玫烟色棒束孢的研究主要集中在生物学[21-22]、生理生化[1,23]及发酵工艺[3]等方面的研究,对其相关功能基因挖掘方面的研究较少。汤强[24]首次从玫烟色棒束孢中克隆了类枯草杆菌蛋白酶Ifupr1基因和几丁质酶Ifuchit1基因的全长cDNA序列。孟祥云[9]对玫烟色棒束孢几丁质酶基因Ifuchit2及上游调控序列进行克隆与分析。Shang等[16]对玫烟色棒束孢ARSEF 2679进行了基因组从头测序,为玫烟色棒束孢基因组学的研究提供了参考。董章勇[11]从玫烟色棒束孢ARSEF 2679菌株全基因组的10 061个蛋白质序列中筛选到了1 112个具有信号肽的蛋白序列。

本研究以对小菜蛾具有高致病性的对玫烟色棒束孢IFCF01菌株为研究对象[25],利用Illumina/Solexa测序平台对其进行转录组从头测序分析,共获得4 868 067 960 nt数据,经组装后获得36 466条总长34 551 468 bp的unigene,unigene平均长度947 bp,N50达到1 676 bp,说明序列组装理想,可满足转录组后续分析的基本要求。然而,当以Shang等[16]组装的玫烟色棒束孢ARSEF 2679基因组序列为参考基因进行有参转录组分析时,发现玫烟色棒束孢IFCF01转录组序列只有9.54%比对上参考基因组,比对率极低,其原因可能为本菌株与玫烟色棒束孢ARSEF 2679菌株存在较大遗传差异,具体原因有待进一步分析。因此,本研究对玫烟色棒束孢IFCF01采用无参转录组分析。已有报道采用无参转录组分析对越橘果实[26]、山茶[27]和瓠瓜[28]等进行物种特异性研究,获得了目标基因。本研究将玫烟色棒束孢unigene 序列与6个公共数据库进行比对,共有62.36%的基因被注释,所获得的序列信息极大丰富了玫烟色棒束孢的序列资源,也为其功能基因研究提供丰富的数据信息。玫烟色棒束孢仍有40% 的基因未被注释,该现象在其他物种[29]中也有存在,原因可能为公共数据库中玫烟色棒束孢同源序列信息相对较少,或这部分序列比较短未得到有效的比对,抑或是玫烟色棒束孢新的功能基因序列,随着高通量测序技术的发展及其在更多真菌物种基因组及转录组研究中的广泛应用,注释到的基因数量会越来越多。此外,在GO数据库中玫烟色棒束孢菌株IFCF01只有9 617条unigene被注释,占总unigene的26.37%,注释率相对比较低,可能与GO数据库不够完善有关[28]。注释上GO的46个亚类中,代谢过程、细胞部分和细胞等功能模块的基因表达比较活跃,这些基因主要为生物体进行各种生命活动所必需的大量表达的持家基因;而免疫系统过程,运动、生长和生物附着等功能模块则包含了生物体中具有特殊功能的低丰度表达基因,推测这些功能模块中基因的表达情况很可能反映玫烟色棒束孢侵染寄主过程中的基因表达时空特异性。KEGG通路注释分析结果发现,参与次生代谢物的生物合成、淀粉和蔗糖代谢、氨基糖和核苷酸糖代谢等途径的基因最多,该发现为玫烟色棒束孢致病相关基因的研究提供了依据。

病原菌-寄主互作基因数据库(pathogen-host interaction gene database,PHI数据库)收集了大量经过实验证实的细菌、真菌、卵菌等病原菌的致病相关基因,主要应用于同源分析各种病原菌潜在的毒力因子[18]。通过将玫烟色棒束孢转录组比对PHI数据库以挖掘玫烟色棒束孢潜在致病相关基因。本次比对发现玫烟色棒束孢共有86个基因比对上PHI数据库,对86个unigene进行COG功能注释发现,这些基因分布在不同的功能模块,还有一些属于未知功能的基因,这些基因是否参与玫烟色棒束孢致病过程中仍有待进一步研究。

目前未见有关玫烟色棒束孢SSR位点的分析报道。本研究从玫烟色棒束孢的36 466个unigenes序列中筛选出8 572个SSR位点,发生频率为18.18%,高于黑木耳(4.75%)[30]和杏鲍菇(8.31%)[31],低于蛹虫草的34.99%[32],其原因可能与转录组的测序与组装方法、SSR筛选的标准等有关。玫烟色棒束孢检测到的 SSR位点中,以三核苷酸重复所占比例相对最大,为47.56%,研究结果与黑木耳和毛木耳[30]以三核苷酸重复为主导相似,而与瓠瓜[28]、石蒜属[33]及杏鲍菇[31]等以单核甘酸重复为主导不同,说明不同的物种之间其主导的重复单元存在差异。据报道,三核苷酸基序在植物表现某种抗逆性时具有积极的促进作用,本研究中玫烟色棒束孢存在接近一半的三核苷酸基序,是否与其具有某种抗逆性有关有待进一步研究。在所有的基序类型中玫烟色棒束孢以A/T含量最为丰富,该结果与金群力等[34]、李强等[31]和许端祥等[28]的研究结果相似。

本研究利用Illumina Hiseq 平台对玫烟色棒束孢IFCF01菌株进行无参转录组测序分析,共获得36 466条unigene,共有62.36%的unigene在Nr、Nt、SwissProt、KEGG、COG和GO这6大数据库中得到注释。其中21 737个unigenes在Nr数据库中得到注释,大部分序列注释到虫生真菌物种中;共有10 888个unigene注释到COGs数据库的25个功能类别中;GO数据库中共有9 617个unigene注释到生物过程、细胞组分和分子功能3大类功能模块中,包含46个亚类;共有14 364个unigene被注释到KEGG的108个代谢途径中;通过比对PHI数据库及根据转录组注释结果,筛选了玫烟色棒束孢潜在致病相关基因;玫烟色棒束孢转录组共有8 572个SSR位点,发生频率为18.18%。本研究获得的玫烟色棒束孢转录组数据、潜在致病相关基因及SSR序列特征,为后续玫烟色棒束孢基因资源挖掘、基因功能研究及SSR分子标记奠定了基础。