蜜蜂球囊菌菌丝和孢子中全长转录本的差异表达分析

杜 宇, 蒋海宾, 王 杰, 范小雪, 王秀娜, 冯睿蓉, 张文德,隆 琦, 熊翠玲, 郑燕珍,2, 陈大福,2, 郭 睿,2,*

(1. 福建农林大学动物科学学院(蜂学学院), 福州 350002; 2. 福建农林大学, 蜂产品加工与应用教育部工程研究中心, 福州 350002;3. 福建农林大学生命科学学院, 福州 350002; 4. 福建农林大学, 福建省病原真菌与真菌毒素重点实验室, 福州 350002)

蜜蜂球囊菌Ascosphaeraapis(简称“球囊菌”)是一种致死性真菌病原，专性侵染蜜蜂幼虫而导致白垩病。该病于1992年在我国全面暴发，并往往发生在凉爽湿润的季节，与蜂群大量繁殖的春季和秋季相吻合，可导致成年蜜蜂数量和蜂群群势的急剧下降，严重影响蜂产品的可持续发展(Aronstein and Murray, 2010; 赵红霞等, 2014)。低日龄蜜蜂幼虫易受球囊菌感染，而成年蜜蜂可作为媒介在蜂群内传播疾病。球囊菌孢子随食物进入蜜蜂幼虫中肠，因此时中肠与后肠隔绝，肠道内低氧环境使孢子处于低水平萌发状态，直至幼虫到预蛹的过渡期，中肠与后肠连通，球囊菌进入后肠接触氧气后剧烈萌发，菌丝开始大量生长，相继穿透肠壁和体壁，最终蔓延包裹幼虫全身，形成白色或棕黑色的坚硬虫尸(Jensenetal., 2013)。诸多生态因子如温度、湿度、酸碱度(pH)、二氧化碳(CO2)浓度、辐射及营养素条件等均对球囊菌的生长和白垩病暴发产生一定影响，外部温度介于(31～35)±0.5℃，相对湿度在80%以上，pH值介于5～7.8，且CO2浓度高于10%的环境条件有利于球囊菌孢子的萌发(梁勤等, 2000, 2001; Vojvodicetal., 2011; Jensenetal., 2013; 董文滨, 2014)。球囊菌侵染可引起宿主的胁迫应答。Aronstein等(2010)利用cDNA-AFLP和qRT-PCR技术研究发现1日龄西方蜜蜂Apismellifera每头幼虫摄入1×105个球囊菌孢子后，其摄食率迅速下降，并影响宿主体内免疫信号途径的激活以及抗菌肽相关基因的转录过程。王立立(2014)研究发现意大利蜜蜂Apismelliferaligustica5日龄幼虫被5×107孢子/mL浓度的球囊菌侵染后，血淋巴内酚氧化酶原级联反应为主的体液免疫被显著激活，并催化黑化反应。

此前，由于基因组基因功能注释信息的缺失，球囊菌的分子生物学研究进展缓慢。Shang等(2016)通过二代测序技术对球囊菌ARSEF 7405菌株进行了基因组测序、组装和注释，公布了scafford水平的球囊菌参考基因组(assembly AAP 1.0)，为其分子生物学及组学研究奠定了基础。近年来，以RNA-seq为代表的第二代高通量测序技术发展迅猛，为探究蜜蜂与球囊菌的互作机制提供了强大技术手段。前人利用二代测序技术对球囊菌进行的转录组研究仅有一例报道，Cornman等(2012)对来自培养基的球囊菌菌丝和来自西方蜜蜂幼虫感染组织的球囊菌菌丝进行了转录组测序，功能分析表明球囊菌的差异表达基因(differentially expressed genes, DEGs)参与了交配类型、细胞内信号转导和应激反应。近年来，笔者所在团队基于二代测序技术对球囊菌进行了一系列转录组研究，例如基于Illumina测序数据组装和注释了球囊菌的首个参考转录组，并大规模开发了球囊菌的SSR(张曌楠等, 2017)；利用二代转录组数据对球囊菌的参考基因组注释进行了补充和完善(郭睿等, 2019)；通过比较转录组分析初步揭示了球囊菌对意蜂幼虫和中蜂幼虫的侵染机制(陈大福等, 2017; 郭睿等, 2017)；并在全基因组水平鉴定和分析了球囊菌的微小RNA(microRNA, miRNA)、长链非编码RNA(long non-coding RNA, lncRNA)和环状RNA(circular RNA, circRNA)(Chenetal., 2018; Guoetal., 2018; 郭睿等, 2018)。二代测序技术虽具有通量和准确性较高的优势，但同时也存在测序读段较短(不超过300 bp)的劣势，需要通过生物信息学算法对测序得到的短读段进行拼接和组装，以得到较长的转录本(Boldogköietal., 2019)。因此，基于二代短读段数据只能进行基因的定量和差异表达分析，无法开展转录本的结构、定量和差异表达分析(Tombáczetal., 2018)。此外，基于二代测序数据的基因定量实际上是对来源于同一个基因的不同剪接异构体的表达量取平均值，因而混淆了不同剪接异构体的真实表达量。因此，在转录本水平进行定量和差异表达分析对于深入理解物种的转录组的复杂性及动态变化尤为必要。

近年来，以牛津纳米孔(Oxford Nanopore)长读段测序技术为代表的三代测序技术逐渐兴起和成熟，因具有超长读段(平均约15 kb)、较高准确性和直接读取核酸修饰等优势而成功应用于人、兔和橄榄实蝇Dacusoleae等物种的全长转录组研究(Chenetal., 2017; Bayegaetal., 2018; Workmanetal., 2019)。此外，由于三代测序得到的超长读段无需拼接便能获取转录本序列全长，使转录本的精确定量和差异表达分析成为现实(Weiratheretal., 2017)。目前，真菌的全长转录组研究十分有限(Piroonetal., 2018)；对于蜜蜂病原，迄今仅在东方蜜蜂微孢子虫Nosemaceranae中有过相关研究报道(陈华枝等, 2020a, 2021a)。近期，笔者利用Nanopore长读段测序技术对球囊菌的纯化菌丝(Aam)和纯化孢子(Aas)分别进行了转录组测序，基于混合数据构建和注释了球囊菌的高质量全长转录组(杜宇等, 2021b)；系统鉴定和分析了Aam和Aas中基因的可变剪接(alternative splicing, AS)和可变腺苷酸化(alternative polyadenylation, APA)(杜宇等, 2021a)。然而到目前为止，球囊菌全长转录本定量及差异表达转录本(differentially expressed transcripts, DETs)分析未见报道。

本研究在前期研究基础上对Aam和Aas中全长转录本进行定量分析，通过比较分析筛选出DETs并进行功能和通路注释，进一步对Aam和Aas中DET的结构进行比较，并构建和分析DETs对应蛋白的互作网络。研究结果可为阐明DETs在球囊菌的孢子萌发、菌丝生长和有性生殖等方面的分子功能提供理论依据和数据基础。

1 材料与方法

1.1 纳米孔长读段测序数据来源

前期已利用Oxford Nanopore技术对球囊菌Aam和Aas进行转录组测序，分别测得6 321 704和6 259 727条原始读段，经过滤分别得到5 669 436和6 233 159条有效读段，并分别鉴定9 859和16 795条非冗余全长转录本(杜宇等, 2021b)。

1.2 全长转录本的表达量计算

利用minimap2软件(Li, 2018)将1.1节Aam和Aas有效转录本读段与球囊菌参考基因组(assembly AAP 1.0)注释的已知转录本进行序列比对，获取全长转录本与已知转录本的对应信息。Nanopore全长转录组测序可以模拟成一个随机抽样的过程，为了让片段数目能真实地反映转录本表达水平，首先对各样品中的比对上球囊菌参考基因组的有效读段数目进行归一化处理，然后采用CPM(counts per million mapped reads)算法(Zhouetal., 2014)计算每一条全长转录本的表达量，公式如下：

CPM=比对到转录本上的读段数(reads mapped to transcript)/样本中对齐的总读段数(total reads aligned in sample)×1 000 000。

1.3 DETs的筛选

将同一全长转录本在Aam和Aas中表达量的比值定义为相对变化倍数(fold change)。通过对差异显著性P值进行校正得到错误发现率(false discovery rate, FDR)。由于Nanopore全长转录组测序的差异表达分析是对大量的全长转录本表达量进行独立的统计假设检验，会存在假阳性问题，因此在进行差异表达分析过程中，采用了公认的Benjamini-Hochberg校正方法对原有假设检验得到的显著性P值进行校正，并最终采用FDR作为DETs筛选的关键指标。 根据fold change≥2且FDR<0.01的标准，采用EBSeq软件(Lengetal., 2013)从AasvsAam比较组中筛选DETs。

1.4 DETs的生物信息学分析

利用百迈克云平台(https:∥international.biocloud.net)的相关工具绘制DETs的曼哈顿图及表达量聚类热图。 利用BLAST工具将DETs比对GO数据库(https:∥www.geneontology.org)和KEGG数据库(https:∥www.genome.jp/kegg/)，获得相应的功能和通路注释。再根据DETs与KEGG通路的注释关系构建调控网络，并利用Cytoscape软件(Smootetal., 2011)进行网络可视化。参照Sylvain等(2007)的方法，通过IGV(Integrative Genomics Viewer)浏览器对DETs的结构进行可视化。

2 结果

2.1 球囊菌菌丝和孢子中全长转录本的表达量分析

Aam和Aas共有的表达转录本数量为19 966条，特有的表达转录本数量分别为1 273和2 856条。Aam中全长转录本的表达量介于0.0010～76 721.3655，其中表达量最高的为ONT.6598.1, ONT.4715.1和ONT.3612.1等(表1)。Aas中全长转录本的表达量介于0.0010～20 116.2770，其中表达量最高的为ONT.7346.1, ONT.2508.21和ONT.2150.1(表2)。此外，部分全长转录本同时在Aam和Aas中高量表达，例如ONT.705.2(CPM值分别为2 290.0785和9 827.5333), ONT.2558.1(CPM值分别为2 917.3948和6 693.9654)和ONT.4966.2(CPM值分别为6 640.7930和5 544.5046)。

表1 蜜蜂球囊菌菌丝(Aam)中CPM值前10位的全长转录本信息概要Table 1 Summary of full-length transcripts in Ascosphaeraapis mycelium (Aam) with the top 10 CPM values

2.2 球囊菌菌丝和孢子中DETs的筛选及表达量聚类

在AasvsAam比较组中共筛选到3 230条DETs，包含3 072条上调转录本和158条下调转录本。 上述DETs的log2(fold change)介于-6.5506～14.4956。上调倍数最大的全长转录本是ONT.4715.1，其次为ONT.7176.3和ONT.3781.2(表3)；下调倍数最大的全长转录本是rna5277，其次为ONT. 2275.2和ONT.7416.7(表4)。

表3 蜜蜂球囊菌Aas vs Aam比较组中前10位上调转录本的信息概要Table 3 Summary of the top 10 up-regulated transcripts in the Aas vs Aam comparison group of Ascosphaera apis

表4 蜜蜂球囊菌Aas vs Aam比较组中前10位下调转录本的信息概要Table 4 Summary of the top 10 down-regulated transcripts in the Aas vs Aam comparison group of Ascosphaera apis

2.3 球囊菌菌丝和孢子中DETs的功能和通路注释

GO数据库结果显示，上述DETs涉及生物学进程(biological process)、细胞组分(cellular component)和分子功能(molecular function) 3个大类，其中有1 251条DETs注释到生物学进程大类的16个功能条目，包括细胞进程(cellular process)、代谢进程(metabolic process)及单组织进程(single-organism process)等(图1)；有1 170条DETs注释到细胞组分大类的14个功能条目，包括细胞(cell)、细胞部分(cell part)及细胞器(organelle)等(图1)；有1 157条DETs注释到分子功能大类的12个功能条目，包括催化活性(catalytic activity)、结合(binding)和结构分子活性(structural molecule activity)等(图1)。

KEGG数据库注释结果显示分别有36, 639和832条DETs注释到细胞进程(celluar process)、遗传信息处理(genetic information processing)和新陈代谢(metabolism)相关通路(图2)。注释DETs数量最多前5位通路分别为核糖体(ribosome)、抗生素的生物合成(biosynthesis of antibiotics)、碳代谢(carbon metabolism)、氨基酸的生物合成(biosynthesis of amino acids)及糖酵解/糖异生(glycolysis/gluconeogenesis)(图2)。

进一步分析发现64条DETs注释到糖酵解/糖异生，包括ONT.1193.1, ONT.1218.4和ONT.29.1等62条上调转录本，以及2条下调转录本(ONT.7961.1和ONT.7961.5)；ONT.5566.7, ONT.5566.3和ONT.4500.10等26条上调转录本注释到内吞作用。5条上调转录本(ONT.1358.7, ONT.1474.3, ONT.750.14, ONT.750.2和rna3824)以及2条下调转录本(ONT.5822.1和rna735)注释到MAPK信号通路；2条上调转录本(ONT.750.14和ONT.750.2)同时注释到内吞作用和MAPK信号通路(图3)。

图3 蜜蜂球囊菌Aas vs Aam比较组中差异表达转录本(DETs)与糖酵解/糖异生、内吞作用及MAPK信号通路的注释关系网络Fig. 3 Annotation relationship networks of glycolysis/gluconeogenesis, endocytosis and MAPK signaling pathway associateddifferentially expressed transcripts (DETs) in the Aas vs Aam comparison group of Ascosphaera apis

2.4 球囊菌菌丝和孢子中DETs的结构比较

对Aam和Aas中DETs的结构进行比较，发现685个DETs(ONT.3612.1, ONT.1081.2, ONT.1218.6和ONT.1534.4等)在菌丝和孢子中除具有不同的表达量外，同时也具有不同的结构，例如全长转录本ONT.3612.1由分别在Aam和Aas中发现的剪接异构体ONT.2373.1和ONT.3407.1合并得到，ONT.3612.1在Aam和Aas中的CPM值分别为33 850.23和47.22；此外，ONT.2373.1和ONT.3407.1的结构存在差异，二者均含有1个外显子，但前者的外显子明显短于后者(图4)。

图4 球囊菌菌丝(Aam)和孢子(Aas)中差异表达转录本 ONT.3612.1结构的IGV浏览器视图Fig. 4 IGV browser view of the structure of differentially expressed transcript ONT.3612.1in Ascosphaera apis mecylium (Aam) and spore (Aas)白色矩形表示scafford，上面的数字表示scafford上的碱基位置；红色矩形、蓝色矩形和绿色矩形表示外显子，箭头表示转录方向；Final代表在Aam和Aas中鉴定到的转录本的合并结果。White rectangle indicates scafford, and the numbers above scafford indicate base position. Red rectangle, blue rectangle and green rectangle indicate exons, and arrows indicate the direction of transcription. Final represents the integration of identified transcripts in Aam and Aas.

3 讨论与结论

前期研究中，笔者所在团队基于二代短读段数据对球囊菌菌丝和孢子转录组中的表达基因进行了定量和比较分析，分别检测到7 557和7 640个表达基因，筛选出7 362个共有表达基因以及195和278个特有表达基因，并鉴定到977个上调基因和687个下调基因(蒋海宾等, 2020)。由于二代测序产生的短读段限制，无法对转录本进行结构、精确定量和差异表达分析。同一个基因转录而来的前体mRNA(pre-mRNA)通过AS可形成不同的剪接异构体，最终形成不同的蛋白质而表现出不同的表型(Modrek and Lee, 2002; Wangetal., 2008)。

相对于菌丝，真菌孢子是一种休眠态，仅保持低水平的新陈代谢和生命活动(Liuetal., 2016)。本研究利用已获得的Nanopore长读段测序数据对Aam和Aas中表达的转录本进行定量，结果显示Aam中全长转录本的CPM值介于0.0010～76 721.3655(表1)，而Aas中全长转录本的CPM值介于0.0010～20 116.2770(表2)，说明球囊菌菌丝中全长转录本的整体表达水平更高，与实际情况相符。此前，我们利用二代测序技术在球囊菌孢子中鉴定到392个miRNA(陈华枝等, 2020b)、850个lncRNA(陈华枝等, 2021c)和2 995个circRNA(陈华枝等, 2021b)，并对部分ncRNA的表达进行了验证，结果说明球囊菌孢子中存在一定水平的基因转录活动。 本研究的结果进一步提供了球囊菌孢子中存在转录活动的证据。本研究发现，部分转录本在菌丝中特异性高表达，例如ONT.1954.6(CPM=1 421.3261), ONT.5520.1(CPM=6 245.3525)和ONT.7176.1(CPM=2 651.4502)等；部分转录本在孢子中特异性高表达，例如ONT.7346.1(CPM=20 116.2770), ONT.2508.21(CPM=17 482.8947)和ONT.2150.1(CPM=16 328.9977)等；且ONT.705.2, ONT.2558.1和ONT.4966.2等转录本同时在菌丝和孢子中高量表达。推测菌丝中特异性高表达的转录本与菌丝的生长和发育息息相关，孢子中特异性高表达的转录本与孢子的萌发和生长存在密切关系，而菌丝和孢子均维持高量表达的转录本在球囊菌生活史的不同阶段皆发挥必要功能。相较于孢子，球囊菌菌丝中分别有3 072和158条全长转录本上调和下调表达，上调转录本的数量远多于下调转录本，说明多数转录本在菌丝中的表达更为活跃，与实际情况相符。本研究还发现一些转录本上调和下调的幅度很大，例如ONT.4715.1[Aam中CPM=75 415.9773, Aas中CPM=11.0906, log2(fold change)=14.4956]和ONT.2430.18[Aam中CPM=1.5644, Aas中CPM=364.8805，log2(fold change)=-5.9176]等，鉴于这些转录本在球囊菌菌丝和孢子中均维持不低的表达量，值得进一步深入研究。

球囊菌在增殖过程中需要从外界摄取多种碳源、氮源、维生素和矿质元素，其中果糖和葡萄糖分别是球囊菌菌丝体生长和杂交产孢的最佳碳源(梁勤等, 2001)。真菌能通过糖酵解途径将摄取到的葡萄糖转化为丙酮酸，从而满足自身生长的能量需求(刘天明等, 2008)。因此，球囊菌在菌丝生长过程中必然进行活跃的营养摄取和能量转化活动。糖酵解与糖异生途径是方向相反的两个能量代谢途径，二者的多数酶促反应可逆。糖异生是一种从氨基酸或三羧酸(TCA)循环的中间体等非碳水化合物底物中生成葡萄糖或糖原的代谢过程，在生物体发育和适应营养饥饿方面发挥着重要作用(Marreroetal., 2010)。真菌在面对葡萄糖或其他碳源不足时能够通过糖异生途径维持自身发育所需的能量，并在侵染期间通过该途径应对摄取外源营养不足等情况(Brakhage, 2013)。Liu等(2018)研究发现糖异生途径相关的PCK1基因可启动灰霉菌Botrytiscinerea分生孢子萌发和侵染结构形成等过程，并协助灰霉菌侵染寄主细胞以及增强病原毒力。本研究中，64条DETs同时注释到糖酵解/糖异生相关通路，包括62条上调转录本(ONT.1193.1, ONT.1218.4和ONT.29.1等)以及2条下调转录本(ONT.7961.1和ONT.7961.5)(图3)，上调转录本数量远高于下调转录本的数量，暗示球囊菌菌丝可通过上调表达糖酵解/糖异生途径相关转录本，影响自身能量转化过程，以满足其旺盛的代谢活动需要。但由于糖酵解和糖异生在KEGG数据库中的pathway信息为糖酵解/糖异生，所以尚无法区分分别注释到糖酵解和糖异生的DETs。

内吞作用涉及细胞的生长和分化以及分子信号传递。真菌菌丝可通过内吞作用将酶类与质膜等关键分子运送至菌丝顶端和其他部位以促进菌丝生长和极性延伸(Wangetal., 2019)。罗伯茨绿僵菌Metarhiziumrobertsii体内的MrArk1作为内吞作用的关键调控因子参与分生孢子和毒力水平的调节过程(Wangetal., 2019)。稻瘟病菌Magnaportheoryzae内吞作用基因MoARK1和MoPAN1可显著影响其细胞壁完整性、产孢量以及致病性(王佳妹, 2012)。但内吞作用在球囊菌增殖过程中的生物学功能目前仍未明确。在球囊菌侵染蜜蜂幼虫肠道过程中，其菌丝分泌各类水解酶抑制或降解宿主的防御因子，并吸收利用宿主的营养物质以促进生长过程(Jensenetal., 2013)。本研究发现ONT.5566.7, ONT.5566.3, ONT.4500.10和ONT.1179.1等26条全长转录本注释到内吞作用(图3)，且在Aam中均为上调表达。推测上述内吞作用相关全长转录本在促进菌丝生长方面发挥一定作用，并与病原毒力水平存在潜在关联。

丝裂原活化蛋白激酶(mitogen-activated protein kinase, MAPK)级联反应可引起下游靶基因的表达变化，参与真菌细胞壁完整性、交配繁殖、孢子形成以及致病性等方面的调控过程(Zhaoetal., 2007; Igbariaetal., 2008)。新月弯孢菌Curvularialunata的MAPK信号通路相关基因Clk1和Clm1分别与病原孢子形成以及细胞壁形成密切相关(Nietal., 2018)。笔者所在团队前期利用二代测序技术发现球囊菌菌丝和孢子中的lncRNA和circRNA差异表达，并分别通过顺式(cis)作用和调控来源基因转录，影响MAPK信号通路相关基因的表达(陈华枝等, 2021b, 2021c)。本研究发现5条上调转录本(ONT.1358.7, ONT.1474.3, ONT.750.14, ONT.750.2和rna3824)以及2条下调转录本(ONT.5822.1和rna735)皆可注释到MAPK信号通路(图3)。其中，ONT.750.14和ONT.750.2可同时注释到内吞作用和MAPK信号通路(图3)。上述结果表明MAPK信号通路相关全长转录本可能与球囊菌的孢子形成、细胞壁完整性维持以及杂交产孢密切相关。

三代测序技术为深入研究转录本结构提供了强大工具(Boldogköietal., 2018, 2019)。基于二代短读段测序数据只能进行基因表达量的计算和差异表达分析，但基于三代长读段测序数据不仅能够同时进行基因和转录本表达量的计算和差异表达分析，还能对基因(转录本)的结构进行精确的AS和APA分析(Abdel-Ghanyetal., 2016; Moldovánetal., 2018)。相对于较早出现的PacBio SMRT测序技术，Nanopore长读段测序技术由于诞生较晚，在转录本结构方面相关研究还较为有限。本研究发现，部分来源于同一基因的剪接异构体在球囊菌菌丝和孢子中具有不同的结构，例如全长转录本ONT.3612.1在菌丝和孢子中经过AS形成两条表达量和结构均不同的剪接异构体ONT.2373.1和ONT.3407.1(图4)，表明球囊菌在生长和发育阶段不仅伴随着转录本表达量的变化，同时也伴随着转录本结构的改变；同一基因来源的全长转录本在菌丝和孢子中经历不同的AS，形成具有不同结构的剪接异构体。这为深入理解球囊菌的菌丝生长、孢子萌发和有性生殖提供了一个崭新的维度，也为其他物种的转录本结构研究提供了可供参考的思路和方法。

球囊菌的转基因操作技术平台缺失导致其基因(剪接异构体)功能研究严重滞后。近期，Tauber等(2020)通过体外转录合成Ras家族编码基因和β-葡聚糖合成蛋白编码基因的双链RNA(dsRNA)并处理球囊菌，发现球囊菌孢子可能在萌发早期吸收上述dsRNA，孢子萌发率明显降低。该研究为深入开展球囊菌的基因(剪接异构体)的功能研究提供了方法借鉴。