松乳菇子实体两个发育时期的转录组分析

宋志强，丁 祥，唐 贤，朱 淼，侯怡铃，*

(1.西华师范大学 生命科学学院，西南野生动植物资源保护教育部重点实验室，四川 南充 637009； 2.西华师范大学 环境科学与工程学院，四川 南充 637009)

食用菌子实体的形成与其自身发育的分子机制息息相关，快速发展的转录组测序技术逐渐揭示了大型真菌发育的分子机制。不同专家学者对食用菌发育分子机制的研究方法各不相同，其中de novo转录组测序[1]因其成本低、周期短、准确性高、能够不受物种是否存在参考基因组的限制等优点被广泛应用。李帆等[2]通过转录组测序分析发现刺芹侧耳菌丝期的抗氧化酶、氨基酸代谢、碳代谢等均会影响刺芹侧耳的生长发育。吴小梅等[3]通过对双孢蘑菇原基期、采收期和开伞后期3个发育时期的转录组测序分析，发现氨基酸代谢会影响双孢蘑菇子实体发育。陈志宏等[4]研究发现，草菇子实体转录组水平测序中的elnL基因编码的bZIP型转录因子，能够通过控制与草菇菌柄伸长有关的一系列基因的表达，影响草菇子实体的生长发育。申进文等[5]通过转录组测序发现真菌转录因子MHR超家族基因中的fst3基因调控糙皮侧耳子实体的生长发育，fst3基因在双核菌丝期、原基期、子实体成熟期的表达量依次增强，子实体成熟期最强。

松乳菇(Lactariusdeliciosus)隶属红菇科乳菇属[6]，子实体早期呈半球形，后渐平展呈波形，中央凹陷，呈鲜艳的黄色、橙黄色至虾仁色，光滑，湿润时黏，有绒毛状同心环带，后变浅色，主要分布于长江中下游及沿海地区，四川云南等地均有松乳菇分布，且松乳菇能与松衫、马尾松等形成菌根，是珍稀的野生食用菌[7-8]。目前，通过de novo测序探究影响松乳菇发育过程中的分子机制还未见报道。本文通过de novo转录组测序，对采自四川省小金县松乳菇子实体幼菇期和成熟期的差异表达基因进行分析，为探究松乳菇在生长发育过程中所涉及到的基因及代谢机制提供一定的理论基础与指导意义。

1 材料与方法

1.1 材料制备

本实验所研究的幼菇期(LDG-1)与成熟期(LDG-2)的松乳菇来源于四川省阿坝藏族自治州小金县不同发育时期的子实体[9]。采集完成后，取小块新鲜松乳菇子实体转移至RNA储存液中保存[10]。

1.2 样品总RNA的提取与检测

松乳菇子实体中总RNA按照OMEGA(OMEGA’s Fungal Total RNA Extraction Kit)真菌总RNA提取试剂盒的使用说明进行提取。使用琼脂糖凝胶电泳、NanoDrop、Qubit 2.0、Agilent 2100等方法分别对RNA降解与污染程度、纯度、浓度以及完整性进行精确定量的检测分析[11]。

1.3 生物信息学分析

将检测合格后的样品总RNA送北京诺禾致源生物信息科技有限公司进行cDNA文库构建，库检合格后进行测序。

为获得高质量分析数据，通过过滤筛选原始测序数据(raw data)，得到分析数据clean reads，利用Trinity[12]对分析数据clean reads进行de novo拼接组装，在此基础上选取每个基因中最长的转录本作为非重复独立基因Unigene，并在KOG/COG数据库将非重复独立基因Unigene进行功能注释。利用RSEM[13]软件进行基因表达定量分析(FPKM>0.3)。对差异基因进行GO富集分析和KEGG富集分析时采用DESeq[14]进行基因差异性表达分析，阈值为Padj<0.05。

2 结果与分析

2.1 转录组的数据质量分析

LDG-1与LDG-2进行测序后分别获得原始数据(raw reads)，去除带接头的reads和低质量的reads后分别获得7.44G(LDG-1)和7.21G(LDG-2)分析数据(clean reads)，其中每个样品的Q20(%)均大于97.5%，Q30(%)均大于93.5%，GC含量在54.7%左右，碱基错误率为0.01%，拼接分析数据clean reads后一共获得129 004个转录本(transcript)，70 223个不同长度的非重复独立基因(unigene)，占总转录本的54.43%。以上结果说明该测序结果可靠性较高，可用于下一步分析研究(表1、图1)。

2.2 转录组基因KOG功能注释与GO富集分析

转录组基因KOG功能注释结果表明，LDG-1和LDG-2中的基因具有全面且复杂的KOG功能类别。同时，用GO分析[15]方法对LDG-1和LDG-2的基因进行了基因功能分类，结果显示，影响松乳菇在生长发育过程中的关键基因功能分类主要集中在生物过程(biological process)中的新陈代谢过程、细胞过程、单体过程；细胞成分(cell components)中的细胞、细胞组分和分子功能(molecular function)中的结合、催化过程。

2.3 两个发育时期基因表达水平分析

通过RSEM软件对基因表达进行定量分析[16]，FPKM大于60的基因高水平表达，而FPKM为0～0.3的基因为低水平表达或根本不表达[17]。分析结果显示，松乳菇在不同生长阶段(LDG-1和LDG-2)基因表达水平存在一定的差异。在松乳菇子实体幼菇期(LDG-1)有34 091个基因表达(FPKM≥0.5)，约占基因总数的74.41%，其中，2 684个基因高表达(FPKM≥60)，约占基因总数的5.86%；而在松乳菇子实体成熟期(LDG-2)，有36 650个基因表达(FPKM≥0.5)，约占基因总数的75.67%，其中，2 607个基因高表达(FPKM≥60)，约占基因总数的5.38%。

LDG-1a、LDG-1b、LDG-1c是幼菇期(LDG-1)的3个重复样品；LDG-2a、LDG-2b、LDG-2c是成熟期(LDG-2)的3个重复样品。
LDG-1a，LDG-1b，LDG-1c were three replicates of young mushroom (LDG-1) stage. LDG-2a，LDG-2b，LDG-2c were three replicates of and mature mushroom(LDG-2) stage.

图1 拼接转录本与Unigenes序列长度分布图Fig.1 Distribution of assembled transcripts and unigenes length

对LDG-1与LDG-2中绝对表达量FPKM大于60的基因进行聚类分析(图2)，结果表明，ATP酶、多功能伴侣蛋白等家族基因为LDG-1和LDG-2主要高表达基因。进一步分析发现，松乳菇子实体幼菇期(LDG-1)有7个基因为极高表达(FPKM>2 000)(表2)，包括H3、H4、H2B、MGAM等；5个基因在松乳菇子实体成熟期(LDG-2)极高表达(表3)，包括MSI、EEF1A、ARD1等。其中，WrbA和UBE2D-E在松乳菇幼菇期(LDG-1)和成熟期(LDG-2)中均极高表达。

2.4 两个发育时期差异表达基因(DEGs)的分析

使用DESeq对LDG-1和LDG-2的差异表达基因(DEGs)进行分析，并设定显著差异表达的阈值Padj<0.05。松乳菇两个发育时期差异表达基因 (DEGs)的分析结果显示，LDG-2期相较于LDG-1期，共存在差异基因16 789个，其中7 635个基因上调和9 154个基因下调。

在|log2FoldChange|>9阈值内的差异基因中，差异倍数较大的10个基因被上调，包括MYO1、PPP1C、SDHC、CNOT3、WBP1、COX10、TAF11、DLD、EARS、MPAO(图3)； 7个基因被下调，包括AmyA、GAD、G6PD、BUD31、RFC2_4、trmB和RP-L5e等(图4)。

2.5 两个发育时期代谢通路KEGG分析

通过KEGG(kyoto encyclopedia of genes and genomes)对松乳菇子实体两个时期的代谢通路进行分析[18]。松乳菇子实体幼菇期(LDG-1)和成熟期(LDG-2)差异基因KEGG Pathway富集结果显示，有3 680个差异表达的基因成功注释在113条通路中。其中，在嘌呤代谢(90个差异基因，2.45%)、MAPK信号通路(90个差异基因，2.45%)、核糖体(325个差异基因，8.83%)、氧化磷酸化(132个差异基因，3.59%)、内吞作用(91个差异基因，2.47%)、内质网中的蛋白质加工(88个差异基因，2.39%)、剪接(82个差异基因，2.23%)、RNA转运(83个差异基因，2.26%)中，差异基因富集显著(图5)。

图3 差异基因的表达：上调(|log2-fold-change|>9)Fig.3 Differential expressed genes: up-regulated (|log2-fold-change|>9)

图4 差异基因的表达：下调(|log2-fold-change|>9)Fig.4 Differential expressed genes: down-regulated (|log2-fold-change|>9)

对KEGG富集通路分析可知，氧化磷酸化途径是松乳菇子实体生长发育时期能量代谢的重要途径。氧化磷酸化不仅具有电子传递、H+传递及氧的利用、产生H2O和ATP的功能，在代谢方面也具有重要作用[19]。富集在氧化磷酸化(oxidative phosphorylation)途径(图6)中的差异基因有35个基因上调，97个基因下调。上调基因主要包括 NADH脱氢酶(泛醌)黄素蛋白1(NDUFV1)、琥珀酸脱氢酶(泛醌)细胞色素b560亚单位(SDHC)、细胞色素c氧化酶4亚单位(COX4)和F型H+转运ATP酶亚单位F(ATPeFF)等；下调基因主要包括NADH脱氢酶(泛醌)1α亚基6(NDUFA6)、琥珀酸脱氢酶(泛醌)黄素蛋白亚基(SDHA)、细胞色素c氧化酶亚单位1(COX1)、F型H+转运ATP酶D亚单位(ATPeF0D)等。

图5 差异基因KEGG pathway富集分析Fig.5 KEGG pathway enrichment analysis of the DEGs

此外，幼菇期(LDG-1)与成熟期(LDG-2)差异基因还显著富集在剪接、内质网中的蛋白质加工、RNA转运和嘌呤代谢中，说明松乳菇子实体幼菇期(LDG-1)在蛋白质合成途径和核酸代谢途径与成熟期(LDG-2)存在着一定差异。通过对KEGG富集通路进一步分析，找到一条MAPK通路中影响子实体发育的关键通路，该通路中(图7)有4个基因上调，包括细胞分裂控制蛋白24基因(Cdc24)、芽出蛋白1基因(Bem1)、鸟嘌呤核苷酸结合蛋白亚单位γ基因(Ste18)和转录因子Ste12基因(ste12)。

3 讨论

本研究选用四川省小金县松乳菇幼菇期(LDG-1)和成熟期(LDG-2)的子实体，通过对其进行de novo测序及生物信息学分析，比较LDG-1和LDG-2在两个生长发育时期基因表达的差异，探究松乳菇在生长发育过程中的关键基因和关键代谢通路。

基因表达水平及聚类分析显示，ATP酶、多功能伴侣蛋白等家族基因为LDG-1和LDG-2主要高表达基因，氧化磷酸化途径中编码ATP酶家族基因[20]中的ATPeFF、ATPeF1G等亚基的表达水平在LDG-1与LDG-2两个时期均有不同程度的变化，提示能量代谢在松乳菇子实体发育过程中起到重要作用。WrbA和UBE2D-E在松乳菇幼菇期(LDG-1)和成熟期(LDG-2)中均极高表达。WrbA为NAD(P)H脱氢酶(醌)[21]，是电子传递链中的第一个蛋白，位于线粒体膜上，该酶极高表达说明在松乳菇幼菇期与成熟期会产生大量能量，加速细胞的分裂，促进子实体的生长发育。UBE2D-E为泛素结合酶e2d/e[22-23]，该酶的极高表达，说明在松乳菇幼菇期与成熟期发育时期光形态建成因子的正常泛素化水平较高，能够调控机体光形态建成，从而控制细胞的分化及结构、功能的改变。

红色，上调；绿色，下调；黄色，既有上调，又有下调。Red， Up-regulation; Green， Down-regulation; Yellow， Both up- and down-regulation.图6 氧化磷酸化Fig.6 Oxidative phosphorylation

红色，上调。Red， Up-regulation.图7 MAPK信号通路Fig.7 MAPK signaling pathway

进一步分析发现，H3、H4、H2B、MGAM和H2A仅在LDG-1中极高表达，H3、H4、H2B、H2A为构成核小体的核心组蛋白[24]，仅在LDG-1中极高表达，提示核小体在松乳菇幼菇期(LDG-1)大量形成，帮助调控松乳菇幼菇期(LDG-1)染色质结构形成。这与王海英[25]的研究结果一致，该研究结果表明，草菇子实体菌柄发育成熟期组蛋白合成受阻，子实体发育成熟期染色体活动不频繁，复制能力减弱，存在着衰老的现象。MGAM是葡糖淀粉酶[26]，属于淀粉酶家族基因，该酶在松乳菇幼菇期(LDG-1)极高表达，提示在松乳菇幼菇期(LDG-1)有快速积累淀粉的过程，能够为幼菇期子实体提供充足的营养，促进生长发育，该结果与吴亚召等[27]发现茶树菇子实体淀粉酶的高峰出现在幼菇期到子实体成熟期的研究结果有一定的相似性。

在成熟期极高表达的eEF1A为翻译伸长因子1a[28]，属于G蛋白家族基因，该基因仅在成熟期(LDG-2)中极高表达，说明该基因能够促进松乳菇成熟期(LDG-2)细胞内受损或错误折叠的蛋白降解，促进转录与翻译，保证子实体正常的生长发育。本研究中G蛋白家族基因EEF1A在松乳菇子实体成熟期呈上调表达，而吴小梅等[3]发现，G蛋白基因在双孢蘑菇子实体发育幼菇期至成熟期呈下调表达，说明松乳菇子实体发育成熟期仍需G蛋白家族基因EEF1A的参与。除此之外，本研究中MSI、ARD1在松乳菇成熟期同样极高表达，MSI是RNA结合蛋白[29]，可以和多种辅助因子相互作用结合不同的靶RNA，仅在成熟期(LDG-2)中极高表达，提示该蛋白以不同的途径来调节松乳菇成熟期的各种代谢；ARD1为D-阿拉伯糖醇脱氢酶(NADP+)[30]，是一种以NADP+为受体的氧化还原酶，该酶仅在成熟期(LDG-2)中极高表达，说明该酶在成熟期(LDG-2)的CO2固定能力高于幼菇期(LDG-1)。以上筛选出的高表达基因对揭示松乳菇在两个不同生长发育期的关键基因提供了一定的科学依据。

通过对本研究中上调与下调差异基因进一步分析，结果显示，在上调基因中，涉及氧化磷酸化中电子传递链的辅基基因SDHC，SDHC是电子传递链中琥珀酸脱氢酶复合体Ⅱ的辅基[31]，提示在松乳菇幼菇期(LDG-1)需要大量能量，氧化磷酸化是其生长发育的关键途径。MPAO为多胺氧化酶[32]，是多胺代谢途径的关键基因，提示多胺氧化酶对于预防松乳菇幼菇期(LDG-1)子实体病变起重要作用。COX10为法尼基转移酶[33]，是预防子实体病变的关键基因，能够抑制FTase的活性，阻断Ras蛋白的法尼基修饰，并有效抑制Ras蛋白的信号转导功能，抑制Ras基因的表达，从而预防松乳菇在幼菇期发生病变。PPP1C为丝氨酸/苏氨酸蛋白磷酸酶pp1催化亚单位[34]，PP1调控了糖原代谢、糖酵解、脂质代谢等生理功能，该基因上调，说明在松乳菇幼菇期在生长发育的能量代谢较为旺盛。在下调基因中，涉及戊糖磷酸途径中的关键基因G6PD，G6PD为葡萄糖-6-磷酸脱氢酶[35]，是磷酸戊糖途径中的一个关键性调控酶，能够调控松乳菇子实体细胞生长，该基因下调，提示在松乳菇成熟期糖代谢中的的磷酸戊糖途径减弱。Cdc24基因为细胞分裂控制蛋白家族基因[36]，该基因上调，说明能够参与调控松乳菇菌丝体的极性生长，这与唐贤等[17]发现，棒瑚菌子实体发育过程中上调的差异基因细胞分裂控制蛋白家族基因Cdc42基因能够调控棒瑚菌子实体菌丝的极性生长的结果一致，说明细胞分裂控制蛋白家族基因在食用菌子实体菌丝极性生长的过程中起重要作用。除此之外，荣成博等[37]发现，在白灵侧耳子实体的生长发育过程中起重要作用的差异基因为编码苯丙氨酸裂解酶的基因(CL4509)以及编码乙醇氧化酶的基因(CL2173)。施肖堃等[38]发现，编码疏水蛋白基因Hyd在双孢蘑菇发育早期至后熟期过程中发挥重要作用。但在本研究中未发现以上3种基因的差异表达情况，说明在松乳菇子实体发育过程中未涉及到以上3种基因的参与。

KEGG富集通路分析发现，氧化磷酸化途径(图6)NADH-辅酶Q等家族基因中既有上调基因也有下调基因，其中NADH脱氢酶(泛醌)黄素蛋白1(NDUFV1)、F型H+转运ATP酶亚单位F(ATPeFF)等基因表达上调；琥珀酸脱氢酶(泛醌)黄素蛋白亚基(SDHA)等基因表达下调，说明氧化磷酸化在松乳菇幼菇期(LDG-1)能量代谢较为旺盛需要大量能量供其生长发育，成熟期则能量代谢趋于平稳，较幼菇期有所降低。

对KEGG富集通路进一步分析，找到一条MAPK通路中影响子实体发育的关键通路，该通路(图7)中细胞分裂控制蛋白24(Cdc24)、Bem1基因(芽出蛋白1)、鸟嘌呤核苷酸结合蛋白亚单位γ(Ste18)和转录因子Ste12(ste12)均表达上调。该通路中信息素α与a信号分子分别与Ste2、Ste3因子受体结合，激活鸟嘌呤核苷酸结合蛋白，磷酸化细胞分裂控制蛋白42(Cdc42)，激活芽出蛋白(Bem1)和p21活化激酶1(PAK1)，再依次激活MAPK激酶激酶激酶(MAP kinase kinase kinase，MKKK)、MAPK激酶激酶(MAP kinase kinase，MKK)和MAPK激酶，进而磷酸化核融合蛋白的转录因子，促进核融合蛋白的表达，使松乳菇进入成熟期(LDG-2)，最终促进松乳菇的有性生殖。同时，细胞分裂控制蛋白42(Cdc42)激活细胞因子分解蛋白(Bin1)，调控松乳菇菌丝的极性生长。综上结果显示，MAPK信号通路在松乳菇子实体的生长发育过程中起重要作用。

综上所述，在松乳菇生长发育过程中，氧化磷酸化是松乳菇子实体生长发育以及能量代谢的重要途径。同时，MAPK信号通路在促进松乳菇有性生殖与调控菌丝体极性生长的过程中起重要作用。