The golden genome annotation of Ganoderma lingzhi reveals a more complex scenario of eukaryotic gene structure and transcription activity

灵芝（Ganoderma lingzhi）的黄金基因组注释揭示了更复杂的真核基因结构和转录活性情况

摘要

背景 普遍认为，真核生物的核基因位于染色体上独立分布，并以单顺反子方式表达。然而，越来越多的证据表明基因结构和转录活动的复杂性超出了这一简单模型。灵芝（Ganoderma lingzhi）是一种具有代表性的药用真菌，但目前缺乏高质量的基因组注释，限制了其遗传学研究的开展。

结果 本研究报告了灵芝的黄金基因组注释，包含14,147个高置信度基因，这些基因经过广泛的人工校正。由此发现了基因结构和转录的新特征。特别是，非经典剪接位点占整个基因组的1.99%，其中GC-AG型（1.85%）、GT-AC型（0.05%）和GT-GG型（0.04%）是主要类型。我们发现1165对基因存在重叠的转录区，且92.19%的基因转录方向相反。我们还在13个灵芝菌株中发现了5412158个基因变异，经过人工校正的基因集提升了这些变异的功能注释。超过60%的灵芝基因表现出可变剪接现象。此外，我们发现两个或更多蛋白质编码基因（PCGs）可以转录成一个单一的RNA分子，这类基因被称为多顺反子基因。总共识别出了1272个与2815个PCGs相关的多顺反子基因。

结论 灵芝中广泛存在的多顺反子基因有力地支持了多顺反子在真核生物基因组中也存在的理论。通过这次黄金基因组注释揭示的非凡基因结构和转录活性，为灵芝以及其他真核生物的基因、基因组及相关研究提供了新的见解。

Here is the translation of the provided text into Chinese:

背景

灵芝（Ganoderma lingzhi），长期被认为与灵芝（G. lucidum）同义，是世界上最著名的药用大型真菌之一。它始终吸引着研究者的兴趣，并作为研究药用真菌的模型生物，持续有关于其药理活性的报道。从2012年首次发布基因组图谱以来，关于灵芝的基因组数据不断增多。目前虽然已有相对高质量的基因组组装数据，但基因组注释的质量仍然受到多种因素的显著影响，包括所使用的软件、重新注释的轮次、参考蛋白数据库和RNA测序（RNA-Seq）数据。准确和完整的基因注释是可靠参考基因组的重要方面，但通常被忽视。在我们之前的研究中，我们发现基于RNA-Seq数据对基因结构进行人工修正，可以提高后续生物信息学分析和基因克隆的准确性。

为了获得更准确和完整的注释，基因组通常会经历多轮重新注释。值得注意的例子包括由社区驱动的曲霉（Aspergillus nidulans）更新、2017年发布的拟南芥（Arabidopsis thaliana）基因组的第11次注释版本，以及2019年野草莓（Fragaria vesca V4）基因组的更新注释。这些更新主要依赖于丰富的RNA-Seq数据、新的注释方法或结合的分析方法，从而识别出大量新的基因。这些不断更新的注释版本极大地促进了科学研究，增强了我们对基因组景观的理解。软件重新注释往往依赖于现有的基因模型，但仍有很大一部分基因池未被充分研究，这可能意味着许多已注释的基因存在错误。此外，基因结构在物种之间甚至同种不同品种之间可能有很大差异，这会导致注释结果不准确。因此，利用转录组数据进行人工修正已成为一种相对精确且高效的基因结构修正方法。

真核生物基因结构的一个关键特征是转录本中存在外显子和内含子。外显子区域保留在最终成熟的mRNA分子中，而内含子区域在转录后加工过程中被剪接出去，这一过程涉及小核RNA和剪接体。GT-AG二核苷酸被称为经典剪接位点或经典内含子，而与GT-AG不同的剪接位点则称为非经典剪接位点。已识别出多种非经典剪接位点，包括GC-AG、AT-AC、GT-TG、GG-AG和AT-AG，其中GC-AG是最常见的非经典剪接位点，其次是AT-AC。虽然植物和人类中非经典剪接位点的分布已得到充分研究，但大型真菌中对此的关注相对较少。

可变剪接（AS），即在前体mRNA中对内含子和外显子的差异性处理，以产生每个基因的多个转录本同种异构体，是真核生物转录组多样化的最重要贡献者。AS通过产生两个或更多不同的蛋白同种异构体来增加蛋白质组的复杂性，并且通过产生截短的蛋白同种异构体有效地导致基因表达的下调。根据RNA-Seq分析，人类中多外显子基因的AS比率超过95%，水稻为42.4%，拟南芥为61%。然而，真菌通常表现出较低的AS比率，根据对23种真菌物种的研究，约6.4%的注释基因受到AS影响，相比早期报道的1.6%–3.6%有所提高。然而，深度RNA-Seq揭示了长孢子曲霉（Trichoderma longibrachiatum）中48.9%的基因表现出可变剪接，表明真菌中的AS被低估了。与较高等生物相比，我们对真菌中AS的理解仍然有限，尤其是在蘑菇型真菌中的普遍性、分子功能和调控机制方面。

普遍认为，真核生物的基因在基因组中是位置上独立的，并转录成不同的RNA分子。然而，最近的研究揭示了更复杂的基因结构和转录活动的情况。真核基因可能与相邻基因有重叠区域，并且在某些情况下，多个基因可以转录为一个单一的RNA分子。类似于原核生物的操纵子，两个或更多的基因可以作为一个多顺反子mRNA转录，或者也可以在真核生物中独立转录，这些多顺反子转录本被称为多顺反子。2015年，Gordon等人首次报告了在高级真菌基因组中多顺反子转录的普遍存在，尤其是在如白腐真菌（Trametes versicolor）和杉木腐菌（Gloeophyllum trabeum）等蘑菇型担子菌中。然而，关于大型真菌中多顺反子基因的系统研究仍然稀缺。

在此，我们报告了灵芝（G. lingzhi）GL0102菌株的两个单倍型基因组，分别命名为GL0102_8和GL0102_53，并通过手动修正所有基因结合全长转录组和RNA-Seq数据，获得了高质量的注释。我们的注释揭示了大量的非经典剪接位点、重叠基因、可变剪接和多顺反子基因。更新后的注释、基因结构的新特征和基因转录为其他大型真菌的研究提供了有用的资源，并加深了我们对大型真菌基因组如何编码遗传信息的理解。

结果

G. lingzhi基因组组装、注释与修正

我们对双核型（GL0102）和单核型（GL0102_8和GL0102_53）菌株进行了Illumina和PacBio Sequel平台的全基因组测序（附加文件1：表S1）。所有菌株的基因组大小估算为43.96–46.74 Mb，双核型菌株的基因组杂合度为1.55%（图1A）。GL0102_8和GL0102_53分别组装成13条染色体，基因组的总长度分别为46.35 Mb（N50 4.97 Mb）和48.56 Mb（N50 4.67 Mb）（图1B，附加文件1：表S2）。GL0102_8和GL0102_53的重复序列分别占13.46%和17.97%（附加文件1：表S3）。通过KAT分析未发现片段丢失（附加文件1：图S1），并且通过proovframe分析显示超过99.4%的高一致性（附加文件1：表S2）。基准单拷贝同源基因（BUSCO）的完整性为99.1%（附加文件1：表S4）。对菌丝体（M）、早期原基（Pe）和晚期原基（Pl）进行了RNA-Seq、同种异构体测序（Iso-Seq）和蛋白质组学分析（附加文件1：表S5），这些数据用于后续的基因预测与修正。GL0102_8和GL0102_53的基因注释使用EuGene [25]工具进行，然后通过Apollo [26]对GL0102_53基因组中的所有基因进行手动修正。根据RNA-Seq和Iso-Seq转录本的支持，基因结构进行了人工修正。修正结果显示，4277个完全匹配的基因保持不变，删除了26,975个外显子，新增6175个新外显子，删除23,681个内含子，新增9718个新内含子，删除6369个位点，新增990个新位点。最终，在手动修正后的GL0102_53基因组中，共注释了14,147个高置信度基因。

基因组调查、组装和基因修正 A. 对灵芝基因组调查和菌丝体的K-mer分析。 B. GL0102_53基因组的全景视图。a，染色体；b，基因密度；c，GC含量；d，重复序列含量；e，GL0102_53与GL0102_8的变异密度对比；f，同源区块（序列长度，≥5 kb）。GL0102_53基因组的平均GC含量为55.72%。值得注意的是，重复序列比例较高的区域显示出较少的基因和较低的GC含量。所有统计数据均通过50 kb不重叠窗口计算。 C. P450基因修正前后的基因结构。箭头表示转录方向。 D. P450基因扩增产物的电泳结果。 E. RNA-Seq读取映射率与GL0102_53的mRNA序列进行比对，分别为修正前和修正后的数据。

在手动修正GL0102_53基因组的过程中，我们发现即使是一些研究较为深入的基因家族，其工具预测的基因结构也存在许多错误。例如，一个P450基因通过将三个独立的基因（chr8g0151541，图1C）组合在一起进行重建。一个P450基因的3'端通过增加四个内含子和四个外显子得到修正（chr6g0087771，附加文件1：图S2A），而5'端则通过增加两个内含子和两个外显子得到修正（chr13g0208881，附加文件1：图S2B）。此外，通过PCR扩增和桑格测序验证了手动修正后的P450基因结构（图1D）。在手动修正后的GL0102_53基因组中，所有基因都进行了UniProt数据库的搜索，发现有393个基因完全没有匹配，其中只有15个有PFAM域注释。在这393个基因中，83个是通过手动修正过程中发现的新基因。

原始和修正后的基因集在多个层面上进行了评估。通过BUSCO分析评估的完整性从修正前的89.7%增加到99.1%（附加文件1：表S4）。在修正后的基因集中的可注释基因比例显著增加。具体来说，修正后的基因集中大约72.63%的基因被EggNOG注释，52.7%被Pfam注释，49.91%被UniProt注释（附加文件1：表S6）。此外，M、Pe和Pl的RNA-Seq数据被映射到原始和修正后的基因集的mRNA序列上。与原始基因集相比，修正后的基因集的映射率提高了近5%（图1E），这表明在使用修正后的基因集作为参考时，基因表达水平的计算更加准确。

在灵芝基因组中发现丰富的非经典剪接位点

剪接位点由一对供体和受体核苷酸组成，定义了外显子和内含子的边界，可以分为经典和非经典两类。在原始基因集中，经典剪接位点占98.01%，而非经典剪接位点占1.99%（GC-AG 1.26%，GT-AC 0.72%，GC-AC 0.01%）。然而，在修正后的基因集中，发现了更多类型和更高比例的非经典剪接位点。具体来说，97.56%的剪接位点是GT-AG，剩下的是94个非经典剪接位点（占2.44%），代表16种供体和16种受体类型（表1和附加文件1：表S7）。在所有非经典剪接位点中，20种类型的比例大于0.01%，其中GC-AG（1.85%）、GT-AC（0.05%）和GT-GG（0.04%）是三种最常见的非经典剪接位点（表1）。非经典剪接位点分布在13条染色体上，涉及1847个基因。含有非经典剪接位点的基因没有表现出显著的功能富集。1653个基因中识别出单对非经典剪接位点，而194个基因中识别出两对或更多对剪接位点。例如，chr5g0072011（SNF2家族N端结构域含量蛋白，附加文件1：图S3A）包含三个GC-AG剪接位点，而chr5g0071881（未表征的蛋白，附加文件1：图S3B）包含两个。

Diversity of G. lingzhi splicing sites

Type	Number	Ratio	Type	Number	Ratio
GT-AG	89,855	97.555000	GG-AG	20	0.021714
GC-AG	1706	1.852190	GT-CG	19	0.020628
GT-AC	49	0.053199	GT-AT	18	0.019542
GT-GG	38	0.041256	CT-AG	17	0.018457
GA-AG	35	0.037999	GT-GT	17	0.018457
GT-TG	22	0.023885	GT-CT	15	0.016285
GT-AA	22	0.023885	AT-AC	13	0.014114
GT-GC	21	0.022800	GT-TC	11	0.011943
GT-CC	21	0.022800	GT-TA	10	0.010857
CA-AG	21	0.022800	Others	156	0.169368
AT-AG	21	0.022800

非经典剪接信号识别不足可能会对基因结构预测和功能注释产生负面影响。 例如，在原始基因集中，由于未能识别AT-AC剪接位点，chr6g0094801的中断导致功能域不完整（图2A）。设计了一对引物，用于扩增chr6g0094801的AT-AC剪接位点的两侧区域（图2B和附加文件1：表S8）。通过桑格测序确认了AT-AC剪接位点以及一个额外的GC-AG剪接位点（图2C）。此外，随机选择了10种非经典剪接位点进行确认（附加文件1：表S8和附加文件1：图S3C和S3D）。

灵芝中的非经典剪接位点 A. chr6g0094801中的AT-AC和GC-AG剪接位点。箭头表示转录方向。zf-C4pol，C4型锌指DNA聚合酶delta。 B. 显示AT-AC和GC-AG剪接位点的PCR产物的琼脂糖凝胶电泳。 C. 包含AT-AC剪接位点的PCR产物的桑格测序。

灵芝基因组中存在显著的重叠基因

现有的基因预测工具设计时假设真核生物基因组中基因不会重叠，而新兴的证据揭示了真核生物中重叠基因的存在[21, 27]。实际上，在本研究或之前的研究中使用的基因预测工具预测的灵芝基因组中并未发现重叠基因[2-5]。然而，在手动修正过程中，通过全长转录本支持，灵芝基因组中发现了大量重叠基因。在GL0102_53基因组中，发现1165对基因具有重叠的转录区域，共涉及2266个基因，分布在13条染色体上（图3A）。这些基因大多是一对一重叠，然而63个基因与两个或三个基因重叠（图3B）。在这些重叠基因对中，1074对表现为相反的基因转录方向。大多数重叠基因在3′末端重叠（73.47%），且重叠区域相对于基因全长的平均比例为26.70%。此外，有146个基因位于其他基因内，并且方向相反。重叠区域的GC含量（54.20%）显著低于基因组中所有基因区域的GC含量（56.82%）（图3C）。这些重叠基因没有表现出显著的功能富集。通过PCR扩增和桑格测序，使用GL0102的cDNA作为扩增模板，确认了重叠基因对chr1g0000801（p450）和chr1g0000811（SET结构域蛋白）之间的重叠区域（图3D）。

重叠基因的统计和验证 A. 重叠基因在不同染色体（chr）上的分布。 B. 重叠基因的分类。 C. 重叠区域和基因区域的GC含量。 "" 表示显著差异，T检验，P < 2.2e−16。 D. 重叠区域的确认。箭头表示转录方向。虚线矩形框表示重叠区域。

灵芝菌株间普遍存在的遗传变异

灵芝的菌株或单倍型在表型上有所不同。例如，当在PDA平板上培养时，GL0102_53表现出较快的生长速度和稀疏的菌丝体，而GL0102_8则表现为较慢的生长速度和更密集的菌丝体（图1A）。高质量基因组的可获得性使得通过比较基因组学进行表型与基因型关系的研究成为可能。在本研究中，12个灵芝基因组与参考基因组GL0102_53进行了比较，发现了总共5,412,158个遗传变异，其中单核苷酸多态性（SNP）最为常见，达4,919,524个（占90.90%）。此外，还有483,420个插入/缺失（indels）和9,214个结构变异（SVs）。在所有变异中，GL0102_53与GL0102_8之间发现了418,180个SNP，53,160个indels和2,257个SV，并通过PCR扩增和桑格测序确认了59个SNP，23个indels和7个SV（图4A和附加文件1：表S9）。在已识别的遗传变异中，722,123个位于或与非基因区域重叠（密度为4.42/100 bp），而4,694,470个变异位于或与基因区域重叠（密度为14.09/100 bp，图4B）。此外，2,399,162个遗传变异位于或与编码区重叠（密度为11.55/100 bp）。与SNP相比，indels和SVs在编码区的显著碱基变化较少发生。总体而言，47.34%的SNP，14.30%的indels和13.60%的SVs位于或与编码区重叠。大量（174,868个）变异位于或与编码区重叠，而在原始基因注释中这些变异位于非编码区，这表明修正后的基因集在变异功能注释方面有所改进。例如，chr7g0119171a的最后一个外显子，包含一个HSP90结构域，并在灵芝菌株之间显示出较高的SNP密度，最初被注释为位于相反链上的基因的UTR（图4C）。保守基因（变异密度<10/100 bp）在应激反应中表现出功能富集（图4D），这表明灵芝在其生存策略上相对保守。相反，对于变异密度较高的基因（变异密度>20/100 bp），未发现显著的功能富集，这表明这些基因在生存中可能具有可有可无的作用，并且与菌株的多样性相关。

灵芝中的遗传变异 A. GL0102_8与GL0102_53之间遗传变异的确认。“8”和“53”分别代表GL0102_8和GL0102_53。红色星号表示SNP位点，红色矩形框表示GL0102_53中的插入序列。 B. 基因区域的变异密度。 C. HSP90基因中的SNPs（chr7g0119171a）。箭头表示转录方向。SNPs显示在10 bp不重叠窗口中。 D. 保守基因的功能富集。计数表示在相同Gene Ontology（GO）条目下注释的基因数量，p.adjust表示统计测试指标，GeneRatio表示在一个条目中被注释的输入基因的比例。

超过一半的基因在灵芝基因组中发生了可变剪接

可变剪接（AS）极大地扩展了转录组和蛋白质组的多样性/复杂性，而准确识别AS则具有挑战性。在高质量注释的基因组和完整转录本的支持下，AS可以被全面且精确地捕捉。为了识别灵芝中的AS，应用了修正后的基因注释和M、Pe、Pl的全长转录本（图5A）。总共获得了2.36 Gb、1.88 Gb和1.97 Gb的Iso-Seq读取数据，分别用于M、Pe和Pl（附加文件1：表S5）。并从Iso-Seq读取数据中生成了137万个全长转录本。随着更多全长转录本的参与分析，能够发现更多的独特基因和独特异构体。当全长转录本的数量超过500,000时，独特基因的数量几乎没有增加（图5B）。相比之下，即使使用了100万个全长转录本，独特异构体的数量仍未达到饱和（图5B），这表明仍有许多异构体有待发现。

灵芝中的可变剪接（AS）特征 A. GL0102不同发育阶段的样本。M，菌丝体；Pe，早期原基；Pl，晚期原基。 B. 唯一基因和异构体检测的饱和测试； C. GL0102_53基因组中AS的统计； D. Gllac7基因的AS。箭头表示转录方向。 E. Gllac7的Western blot分析。

共识别出217,321个异构体，属于15,266个独特基因，分为八类。在这些类别中，不完整的剪接匹配最多，占36.86%，其次是“新型未在目录中”（29.41%）、“完全剪接匹配”（17.53%）、“新型在目录中”（10.69%）、“反义”（0.42%）、“基因间”（0.07%）、“基因组内基因”（0.04%）和其他类型（4.97%）（图5C）。在GL0102_53基因组中，5418个基因未表现出AS，而9848个基因具有两个或更多的剪接异构体，5470个基因具有超过十个的剪接异构体。整个基因组基因集的平均AS率为每个基因14.24。

在一些AS检测中较为困难的案例中，全长转录本数据表现出更优的性能。例如，灵芝的一个木质素酚氧化酶（Gllac7）在原始基因预测中未检测到或注释任何AS事件[28]，而基于全长转录本发现了两个异构体（Gllac7.1和Gllac7.2），其5'端发生了变化（图5D）。由于Gllac7.1完全包含在Gllac7.2中，这两个异构体在其重叠区域具有相同的序列，使得通过短读测序难以区分。除转录本长度的差异外，这两个异构体还编码不同大小的肽，这一点通过Western blot分析得到了支持（图5E）。此外，通过PCR扩增和桑格测序确认了十个基因的AS事件（附加文件1：表S10和附加文件1：图S4）。

Iso-Seq转录本揭示了灵芝中广泛存在的多顺反子基因

普遍认为真核生物不具备多顺反子基因，但在基因修正过程中，通过Iso-Seq转录本的支持，发现了1272个位点转录成单一分子，重叠了两个或更多的蛋白质编码基因（PCGs）（图6A）。这些位点被称为多顺反子基因。或者，这些位点中的PCGs也可以独立转录，Iso-Seq转录本和蛋白质组学数据支持这一点。总的来说，这些多顺反子基因与2815个PCGs相关。具体来说，1011个多顺反子位点与两个PCGs重叠，235个与三个PCGs重叠，22个与四个PCGs重叠。多顺反子mRNA的长度范围从1028到10,459 bp，平均长度为3760.17 bp。多顺反子基因中PCGs的平均蛋白质长度（379.11 aa）显著短于单顺反子基因（449.86 aa）（图6B）。多顺反子基因对（平均距离433.98 bp）显著比其他相邻基因对（平均距离1356.41 bp）更接近（图6C）。在连接单个基因的多顺反子连接区，我们发现一些基序频繁出现。例如，富含TA的基序TACTTAYA在14.6%的分析序列中出现，而富含GC的基序CGCCGCCGCCGYYC在13.2%的序列中出现（图6D）。在几乎一半的所有基因中使用了Opal终止密码子（TGA），在多顺反子基因和单顺反子基因之间仅观察到细微差异（图6E）。我们发现，当评估上游和下游PCGs的相对阅读框时，约三分之一的多顺反子基因在框架内（图6F）。在多顺反子重叠PCGs中，1258个（44.69%）通过PFAM数据库得到了功能注释，富集了F-box-like和P450结构域（图6G）。有趣的是，有七个多顺反子基因由两个串联分布的P450基因组成。随机选择了五个多顺反子基因，并通过PCR扩增和桑格测序进行了确认（图6H和附加文件1：表S11）。

灵芝多顺反子基因的特征 A. Iso-Seq转录本支持的多顺反子基因视图。Iso-Seq转录本支持多顺反子基因和单顺反子基因的存在。肽段证据和Iso-Seq转录本表明，上游和下游基因是分别转录和翻译的。 B, C. 多顺反子基因和单顺反子基因之间蛋白质长度和基因间长度的比较分析。 D. 多顺反子基因的连接区的基序。 E, F. 多顺反子基因与单顺反子基因在终止密码子的使用频率和阅读框上的比较。 G. 多顺反子基因中的蛋白质编码基因（PCGs）的功能富集。计数表示在相同PFAM家族下注释的基因数量，p值表示统计测试指标，GeneRatio表示在PFAM家族中注释的输入基因的比例。 H. 多顺反子基因的PCR产物琼脂糖凝胶电泳。预测的目标PCR产物的大小分别为3610、2675和3019 bp。 "poly" 和 "mono" 分别表示“多顺反子基因”和“单顺反子基因”，"up" 和 "down" 分别表示“多顺反子上游基因”和“多顺反子下游基因”。"" 表示显著差异，T检验，P < 2.2e−16。

讨论

高质量的基因组注释，表现为完整且连续的基因组组装和准确的基因模型，能够极大促进生物学研究。随着测序技术和组装工具的发展，各种生物体的高质量基因组组装已逐步问世。现有多个用于基因模型预测的工具，包括AUGUSTUS [29]、GeneMark [30]、MAKER [31] 和EuGene [25]。基因注释通常依赖于现有的基因模型，但在庞大的基因池中，许多基因仍未得到充分探索，这导致了先前注释基因可能存在不准确的情况。此外，许多基因在物种之间甚至在个体间表现出遗传差异，因此很难保证基因注释的准确性。自动注释算法尚未能正确识别基因组中的所有基因，通常需要人工注释以获得准确的基因模型和基因集。由于人工注释耗时，通常只有基因组中少部分基因模型是人工注释的。近年来，基因组大规模或全基因组的手动修正受到了研究人员的广泛关注 [32]。已有多个灵芝物种的高质量基因组发布 [5]，而本研究中获得的手动修正基因集可以为这些基因组的基因预测提供出色的证据和训练/验证集。

本研究中，我们利用Iso-Seq和RNA-Seq转录数据，手动修正了GL0102_53的全基因组基因集，实现了灵芝基因组的超高质量注释。在基因修正过程中，Iso-Seq数据相较于RNA-Seq数据具有更大的优势，更精确地确定了基因结构、基因边界以及5′和3′-UTRs。然而，Iso-Seq数据在识别非常短的外显子方面效果不如RNA-Seq数据，可能是由于在长读对齐中短外显子的错位（例如，chr5g0072011，附加文件1：图S3）。在手动修正过程中，我们发现许多由软件预测的基因结构存在问题。即使是一些研究较为深入的基因，如P450，依赖于软件的注释仍然容易出错。因此，在手动修正过程中发现了83个新基因，这表明通过细致的基因检查可以发现更多的基因资源。此外，修正后的基因集还可以用于更准确地计算基因表达水平。这些发现强调了高质量基因组注释的重要性，这对后续分析至关重要，包括进化研究、发育研究、基因表达分析和品种鉴定。尽管已获得精确的基因结构和注释，但由于当前数据库的限制，大部分基因被注释为“功能未知基因”。未来对这些高表达或差异表达的未知功能基因进行功能研究至关重要。

本研究首次深入研究了灵芝非经典剪接位点的全景，并识别了94种非经典剪接位点。这种丰富的非经典剪接位点类型增加了基因注释的难度。GC-AG和AT-AC是已知的两种主要非经典剪接位点，基于非经典剪接位点的存在与否，将基因分为四类：GT-AG、GC-AG、AT-AC和少数非经典剪接位点基因[33]，而在本研究中，GC-AG（1.85%）、GT-AC（0.053%）和GT-GG（0.041%）是最为丰富的非经典剪接位点，AT-AC仅占0.014%。这些发现表明，与植物和人类基因组中观察到的非经典剪接位点类型和比例相比存在显著差异[11, 12]。

我们了解到，基因间编码区域的重叠在病毒基因组中是常见的。然而，真菌中也有已验证的基因重叠。例如，在白色念珠菌中，CCT8编码区与TRP1基因的编码区有13个碱基对重叠[21]。此外，哺乳动物（如小鼠和人类）中也有验证的基因重叠[34, 35]。在本研究中，通过手动修正，发现了大量的重叠基因，表明灵芝基因组的简单性。这提醒我们，在未来的基因组注释研究中不能忽视或忽略重叠基因。在原核生物中，单向（即同链）重叠最为常见，而在真核生物中，反向或反平行链重叠更为常见[34]。在本研究中，类似的结果是大多数重叠基因对表现为相反的基因转录方向（不同链重叠）。值得注意的是，如果RNA-Seq数据来自短读、链非特异性的文库，它可能无法准确量化这些重叠区域或基因。只有长读测序技术才能更精确地量化这些基因。

至今，关于丝状真菌的剪接异构体的综合分析仍然缺乏。在以往的研究中，AS通常基于短读RNA-Seq数据识别。然而，短读RNA-Seq数据在AS的识别上有缺点，因为它们不适合准确重建全长剪接异构体。近年来，PacBio Iso-Seq被用于揭示丝状真菌中的AS全景[36]。在我们的研究中，我们利用全长转录本进行AS检测，发现灵芝中64.51%的基因发生了AS，这一比例显著高于以往报道的真菌中AS的比例[20]。如果有更多的全长转录本可用，将能识别更多的AS，表明宏基真菌中的AS比率被低估了。AS显著增加了灵芝转录组的复杂性，拓展了我们对宏基真菌RNA剪接调控的理解。考虑到AS的重要性，有人建议将AS作为标准分析内容，与基因表达分析一起进行[37]。然而，AS在灵芝的生长、发育、表型复杂性和环境适应性中的作用仍需进一步研究。

历史上，人们一直认为真核生物中的每个信使RNA编码单一蛋白质。在本研究中，我们首次发现了大量的多顺反子基因，突显了它们的生物学重要性，并增强了我们对灵芝基因表达的理解。多顺反子基因的存在是否是基因表达调控的一种策略，识别能促进多顺反子表达的广谱基因间信号，需要进一步的实验表征。在体外的多顺反子转录本配对共表达可能有助于灵芝的基因工程研究和工业应用。

遗传变异的维持具有适应性意义[38]。在本研究中，我们在不同的灵芝菌株之间鉴定了大量的遗传变异，其中许多变异位于编码区。这些丰富的遗传资源和基因型变异在表型多样性和环境适应性方面对该物种有显著贡献。

以灵芝基因组为例，真菌基因组的复杂性超出了我们的预期。通过广泛的手动修正获得的高质量基因注释为灵芝基因组学相关研究提供了基础。随着高质量注释基因组和基因转录特征的获得，灵芝可以作为其他真菌物种的研究模型。

结论

总之，通过广泛的手动修正，获得了灵芝（G. lingzhi）14,147个高置信度基因的黄金注释，并且识别了基因结构和基因转录的新特征。灵芝基因组中存在丰富的非经典剪接位点和具有重叠转录区的基因。超过60%的灵芝基因表现出可变剪接现象。此外，还识别了1272个与2815个蛋白质编码基因（PCGs）相关的多顺反子基因。通过黄金注释识别出的非凡基因结构和转录活性为医学真菌研究提供了宝贵的见解。

方法

菌株、培养和样本采集

双核型灵芝（G. lingzhi）菌株GL0102（与福建农林大学真菌学研究中心的“Zhi 102”菌株相同）和由其衍生的两个单核型菌株（GL0102_8和GL0102_53），其交配类型相反，均保存在4°C的土豆葡萄糖琼脂（PDA）平板上。此外，本研究还使用了五个由原生质体单核性生长衍生的单核型灵芝菌株：GL0001_P5、GL0002_P2、GL0002_P3、GL0003_P1、GL0004_P6和GL0005_P3。所有菌株均在28°C的PDA平板上培养。

在PDA平板上培养5天的GL0102菌丝体（M）被收集并迅速冷冻在液氮中，其余菌丝体被接种到含有1 kg培养基的培养袋中。培养基由10%的橡木屑、70%的甘蔗渣、19%的麦麸、1%的石膏和最终含水量为60%的水组成。培养袋在28°C、50% ± 5%湿度下黑暗条件下孵育，待菌丝完全覆盖后转移至结果室。在结果室内，温度保持在28°C ± 5°C，室内湿度为85% ± 10%，并采用12小时光/暗周期，光照强度为300 lx。收集了三份Pe和Pl样本并迅速冷冻在液氮中。

DNA提取、基因组测序、RNA-Seq与Iso-Seq

用液氮研磨的样本进行DNA提取，使用FineOut DNA试剂盒（GENFINE Biotech (Beijing) CO., LTD, O301）并按照制造商的说明进行操作。使用NanoDrop 8000、Qubit和Femto Pulse评估DNA的纯度、浓度和完整性。保留浓度≥80 ng/μL、总量≥12 μg、Femto Pulse上的主要条带大小≥30 kb的基因组DNA（gDNA）。对于GL0102_8和GL0102_53，8 μg的gDNA与16-bp条形码序列结合，然后构建20 kb插入大小的文库。文库在PacBio Sequel II平台的一个SMRT芯片上进行测序，由Annoroad Gene Technology Co., Ltd，Beijing提供。此外，来自GL0102、GL0102_8、GL0102_53、GL0001_P5、GL0002_P2、GL0002_P3、GL0003_P1、GL0004_P6和GL0005_P3的每个菌株的10 μg gDNA用于构建插入大小为300 bp的双端文库，文库在Illumina NovaSeq平台上进行测序，由Annoroad Gene Technology Co., Ltd，Beijing提供。

M、Pe和Pl样本用于RNA-Seq和Iso-Seq（M的RNA-Seq数据来自我们的先前研究[28]）。每个样本的总RNA提取和质量控制使用先前报告的方法进行[39]。简而言之，液氮研磨的样本使用RNA试剂盒（OMEGA，R6827-01）进行RNA提取，并按照制造商的说明操作，保留RNA完整性编号（RIN）≥7.5的RNA样本。RNA-Seq文库构建和测序按照MGI测序平台的协议进行，每个样本生成至少6 Gb的150-bp双端读取。Iso-Seq文库使用SMRTbell prep试剂盒按照制造商的协议构建。简而言之，合成和扩增cDNA，为每个样本构建SMRTbell文库。所有SMRTbell文库被混合并在一个SMRT芯片上进行测序，使用PacBio Sequel II平台，由Annoroad Gene Technology Co., Ltd，Beijing提供。

基因组组装

使用FastQC（Babraham Bioinformatics - FastQC A Quality Control tool for High Throughput Sequence Data）和SeqKit（SeqKit - Ultrafast FASTA/Q kit）评估原始读取质量。低质量的碱基或读取被Skewer [40]过滤，采用以下标准：修剪3′端碱基以达到质量>30，排除长度<100 bp或平均质量<30的读取。使用GenomeScope 2.0 [41]估算基因组大小。PacBio数据通过Canu v1.8 [42]进行组装，并使用Racon [43]和Pilon [44]进行修饰。组装基因组的准确性和完整性通过K-mer分析工具包（KAT）[45]、proovframe [46]和BUSCO分析[47]进行评估，使用真菌odb10数据库。

重复序列

通过结合de novo预测和基于同源性的搜索方法，检测DNA水平的分散重复序列。使用RepeatModeler v2.0.1（http://www.repeatmasker.org/RepeatModeler/）构建de novo重复库，然后将de novo库与Repbase（真核生物重复元素数据库）混合，使用RepeatMasker v4.1.0（RepeatMasker Home Page）进行重复搜索。

基因注释

GL0102_8和GL0102_53基因组的初步注释使用EuGene 4.2 [25]进行。随后，在Apollo [26]中对GL0102_53基因组的整个基因集进行手动修正，基于Iso-Seq和RNA-Seq数据与基因组序列的比对。具体来说，根据RNA-Seq和Iso-Seq转录本的支持，调整基因边界和剪接位点，使用包括删除、合并、分割、创建、扩展、收缩等操作（附加文件1：图S5）。使用Gffcompare v0.10.4（GffCompare）比较手动修正前后的基因注释。通过搜索以下数据库对PCGs进行功能注释：Pfam [48]、UniProt（UniProt）、eggNOG [49]和InterProScan [50]。

遗传变异的识别：灵芝菌株间和菌株内的变异

公共数据库中已有八个以灵芝（G. lingzhi）或灵芝（G. lucidum）命名的基因组，其中五个基因组与GL0102_53（基于ITS2片段）共享超过95%的相似性（附加文件1：表S12）。此外，本研究还包括了六个重新测序的灵芝单核型菌株（GL0001_P5、GL0002_P2、GL0002_P3、GL0003_P1、GL0004_P6和GL0005_P3）用于遗传变异分析。基因组全基因组比较使用Minimap2（GitHub - lh3/minimap2: A versatile pairwise aligner for genomic and spliced nucleotide sequences）进行，GL0102_53作为参考基因组，变异使用paftools调用。识别了SNPs、插入/缺失（indels，长度小于50 bp）和结构变异（SVs，长度大于或等于50 bp）。

可变剪接（AS）和多顺反子基因的识别

通过Iso-Seq管道（GitHub - PacificBiosciences/IsoSeq: Iso-Seq - Scalable De Novo Isoform Discovery from Single-Molecule PacBio Reads）从原始Iso-Seq读取数据生成全长转录本。基于来自GL0102的全长转录本数据，使用SQANTI3 [51]对GL0102_53的全基因组AS进行识别。分类类别包括：完全剪接匹配（完全匹配所有剪接连接）、不完全剪接匹配（仅部分匹配参考剪接连接）、目录中新型（具有新组合的已知剪接位点的新异构体）、目录外新型（至少具有新剪接位点的异构体）、反义（与注释基因的反向基因）、基因内内含子（位于内含子内）、基因内基因组（重叠内含子和外显子）和基因间（位于基因间区域）。多顺反子基因，重叠两个或多个PCGs，基于全长转录本数据在Apollo [26]中手动构建。通过MEME [52]识别了多顺反子基因连接区富集的基序。

非经典剪接位点、重叠基因、遗传变异、AS和多顺反子基因的验证

随机选择包含非经典剪接位点、重叠基因、遗传变异、AS和多顺反子基因的区域进行分析。根据非经典剪接位点两侧的序列设计了18对引物，用于扩增24个非经典剪接位点，分为10类（附加文件1：表S8）。每个重叠区域内的两个基因序列设计了四对引物，扩增了两对重叠基因（附加文件1：表S13）。共设计了23对引物，用于扩增59个SNP、24个indels和8个SV（附加文件1：表S9）。此外，设计了11对引物，扩增了11个表现出AS的基因（附加文件1：表S10）。设计了五对引物，用于扩增五个多顺反子基因（附加文件1：表S11）。PCR扩增反应体系包括1 μL cDNA或gDNA、10 μL 2× Taq Master Mix（Vazyme Biotech Co., Ltd）、0.4 μL 10 μmol/L的正向和反向引物以及8.2 μL ddH2O。PCR反应包括初步变性步骤，95°C下3分钟，然后进行34个循环：95°C下30秒、55-60°C下30秒、72°C下30-90秒，最后以72°C下5分钟结束。对于非经典剪接位点的检测，使用GL0102的gDNA和cDNA作为模板进行扩增。相比之下，重叠基因、AS和多顺反子基因的扩增仅使用GL0102的cDNA作为模板，而GL0102_53和GL0102_8的gDNA用于检测遗传变异。所有PCR扩增产物通过琼脂糖凝胶电泳分析，并随后由广州IGE生物技术有限公司进行测序。

Western blot分析

使用Total Protein Extraction Kit（Solarbio，EX1100）按制造商的说明提取M样本的蛋白质。使用Bradford蛋白定量试剂盒（Solarbio，PC0010）检测蛋白浓度，将8 μg蛋白加载到10% SDS-PAGE凝胶上。电泳后，蛋白质转移到聚偏二氟乙烯（PVDF）膜上，并根据免疫印迹协议使用抗Gllac7抗体进行Western blot分析。

基因功能富集分析

Pfam结构域和GO富集分析使用clusterProfiler [53]进行，p值<1e−3的富集结果被保留。

蛋白质组学

本研究中使用的蛋白质组学文库是从M、Pe和Pl的总蛋白提取物通过胰蛋白酶消化制备的。消化后的肽段通过Pierce高pH反相分级分馏试剂盒（Thermo Scientific）分为六个分级。然后，在Q Exactive质谱仪（Thermo Scientific）上进行LC–MS/MS分析。蛋白质组学数据使用MaxQuant v2.4.14.0进行分析，鉴定了11,161个肽段，分别对应2826个蛋白质。