Breeding of H. erinaceus with high intracellular polysaccharide production
The mutagenic strains HEB and HEC which screened from hundreds of strains bred through ARTP, had been identified by an antagonism test (Fig. 1A) and random amplification polymorphic DNA (RAPD) analysis to be the new strains with genetic material change compared to theoriginal strain HEA [22]. The mutant strains HEB and HEC showed excellent genetic and morphological stability through five generations of successive transfer culture. The biomass of liquid fermentation of strain HEB and HEC was higher than that of the original strain HEA, with the increase rate of 25.96% and 30.37%, respectively. The polysaccharide content in the fermented mycelia of the mutant strains HEB and HEC was also higher than HEA, increased by 23.25% and 47.45% (Fig. 1B).
In our previous study, the crude intracellular polysaccharide fractions of HEA, HEB, and HEC mycelia were extracted by hot water and isolated by gradient alcohol precipitation to 20% and 60% successively. Studies on the physicochemical characteristics of these polysaccharide fractions showed that the 20% ethanol precipitated polysaccharide fraction of the mutants had higher molecular weight than that of the original strain, and the proportion of macropolysaccharide was increased. The proportion of glucose and mannose in the polysaccharide components of the mutants significantly increased than that of the original strain [23]. Also, the activities of macrophages were enhanced by stimulation of 20% ethanol precipitated polysaccharides from HEC induced by ARTP mutagenesis compared with HEA [23]. An obvious different polysaccharide fraction X10-H3P20 between HAE and HAC was revealed by high performance size exclusion chromatography equipped with multiple angle laser light scattering and refractive index detectors (HPSEC-MALLSRI), as shown in Fig. 1C. The molecular weight of this purified polysaccharide X10-H3P20 was about 1.056 × 106 Da (Fig. 1D), and the monosaccharide composition was mainly composed of glucose with ratio of 92% (Fig. 1E) and meanwhile with a β-configuration glycosidic bonds showed by IR spectrum (Fig. 1F). This further indicated that ARTP mutagenesis resulted in synthesis of macromolecule dextran and thus improved the immune activity of polysaccharide in vitro.
Genome sequencing and general features
The H. erinaceus genome sequences from the single Molecule, Real-Time technologies were assembled into 20 scaffolds with an N50 of 258.72 kb and a total genome size of 38.16 Mb (Fig. 2 and Table 1). Prediction of the assembled genome sequence generated 9,780 gene models. The average length of coding genes was 1,355, and the ratio of the total length of the coding region to the whole genome was 34.74%. The average size of exons was 235 bp, and the average size of introns was 70 bp. The 7,137 genes encoded proteins with homologous sequences in the NCBI nr protein databases, and 6,854 genes were mappable through the KEGG pathway database (Table 1). Functional annotation analysis showed the general features, such as 5,611 conserved protein domains (containing 333 CLAN), 2,831 proteins involved in different pathways, 5,611 proteins divided into different GO terms, and 1,822 proteins assigned to different KOG classes in Table 1.
Table 1
General features of the H. erinaceus genome
General features | number | General features | number |
Size of assembled genome (Mb) | 38.16 | Pfam (genes) | 5,611 |
N50_Length (Kb) | 258.72 | Pfam (CLAN) | 333 |
GC content (%) | 53 | SwissProt | 2,267 |
Length of classified repeats (%) | 5 | KEGG alignment | 6,854 |
Number of predicted gene models | 9,780 | GO assignment | 5,611 |
Average exon size (bp) | 235 | KOG assignment | 1,822 |
Average intron size (bp) | 70 | TCDB | 277 |
Average gene length (bp) | 1,355 | DFVF | 360 |
% of Genome (Genes) | 34.74 | PHI | 406 |
% of Genome (internal) | 65.26 | P450 | 95 |
Number of tRNA genes | 204 | Secretory_Protein | 397 |
Number of Contigs | 20 | CAZy | 259 |
NR alignment | 7,137 | Secondary metabolism clusters | 19 |
The phylogenomic analysis showed that H. erinaceus had the closest evolutionary affinity with Dentipellis sp (Fig. 3A). The two species were located at the cluster of Russulales and shared a common ancestor with Polyporales. Analysis of gene family size showed that the net value was − 6919 at the node leading to Russulales and Polyporales, indicating a large number of gene loss during the evolution of Russulales and Polyporales (Fig. 3B). Venn analysis of gene families with big size showed no specific GO annotation of H. erinaceus compared to those in two species in the same cluster based on the phylogenomic tree (Fig. 3C). WEGO analysis of the amplified genes ( > = 10) of H. erinaceus showed that metabolic process, primary metabolic process, and other types of metabolic processes belong to the enriched biological process. Binding and different kinds of binding occupied the most enriched terms of molecular function (Fig. 3D).
Comparative transcriptome analysis of H. erinaceus
Transcriptome analysis showed 2068 differentially expressed genes (DEGs) in HEB_vs_HEA, and 1218 DEGs in HEC_vs_HEA (Fig. 4A and B). Venn analysis showed 768 differentially co-expressed genes among the comparison groups of HEB_vs_HEA and HEC_vs_HEA (Fig. 4C). Heatmap analysis of DEGs showed that HEB and HEC were clustered together (Fig. 4D). GO enrichment annotation of DEGs showed that HEB and HEC had the similar most enriched GO entries, such as biological process, metabolic process, single-organism metabolic process (Fig. 4E and F). The GO entries showed that oxidoreductase, catalytic activity were both enriched in HEB_vs_HEA and HEC_vs_HEA (Fig. 4E, and F), which might be closely related to the synthesis of polysaccharides according to the previous reports[24]. The GO term of cellular components was also enriched in HEC_vs_HEA, such as ribosome, ribonucleoprotein complex (Fig. 4F).
KEGG pathway enrichment using KOBAS (2.0) showed that the significantly upregulated expressed genes in HEB_vs_HEA were enriched in the starch and sucrose metabolism, carbon metabolism, and pyruvate metabolism (Additional file 1). The significantly upregulated expressed genes in HEC_vs_HEA were enriched in glycerolipid metabolism, starch and sucrose metabolism, carbon metabolism, pyruvate metabolism, glycolysis/gluconeogenesis (Additional file 2). The significantly downregulated expressed genes in HEB_vs_HEA or HEC_vs_HEA were both enriched in the ribosome (Additional file 3 and 4).
STRING enrichment showed that the significantly co-upregulated genes in HEB_vs_HEA and HEC_vs_HEA were enriched in metabolic pathways, carbon metabolism, pentose and glucuronate interconversions, pyruvate metabolism (Additional file 5A and B). The significantly co-downregulated expressed genes in HEB_vs_HEA and HEC_vs_HEA were enriched in the ribosomal pathway (Additional file 5 C and D).
These results indicated that the upregulated pathways presented in mutant strains HEB and HEC were involved in carbohydrate metabolism, and the downregulated pathways were strictly associated with protein translation.
Comparative proteomics analysis of H. erinaceus
Results of protein concentration determination using the BCA method confirmed that protein concentration decreased obviously in HEB and HEC compared to HEA (Fig. 5A), in agreement with the downregulated mRNA expression in the ribosomal pathway in HEB and HEC (Additional file 3 and 4, Additional file 5 C and D). The principal component analysis showed that HEA, HEB, and HEC had excellent repeatability and discrimination (Additional file 6A). According to the standard Score Sequest HT > 0, unique peptide ≥ 1, and the blank value was removed, 4,555 trusted proteins were screened (Additional file 7). Proteomics analysis identified 343 DEGs in HEB_vs_HEA and 266 in HEC_vs_HEA (Fig. 5B and C). The details of the DEGs could be found in Additional file 8 and 9. Venn analysis showed that 122 differentially co-expressed proteins in HEB_vs_HEA and HEC_vs_HEA (Additional file 6B).
Heatmap analysis showed the strong enrichment pathways from the significantly upregulated proteins in HEB_vs_HEA and HEC_vs_HEA, such as pyruvate metabolism, glyoxylate and dicarboxylate metabolism, and glycolytic/gluconeogenesis (Fig. 5D). The strong enrichment pathways were obtained from the significantly down-regulated proteins in HEB_vs_HEA and HEC_vs_HEA, such as longevity regulation, peroxisome, and MAPK signaling pathway (Fig. 5D). The details about all the enriched KEGG pathways in HEB_vs_HEA or HEC_vs_HEA could be found in Additional file 10–13. The 18 co-upregulated proteins from the enriched pathways in Additional file 10–13 were enriched in the pathways, such as carbon metabolism, glycolysis/gluconeogenesis, pyruvate metabolism (Fig. 5E), which conformed to the transcriptome analysis results (Additional file 1 and S2). The 22 co-downregulated proteins from the enriched pathways in Additional file 10–13 were enriched in the pathways of peroxisome (Fig. 5E).
The KEGG mapping of the pathway modules of carbohydrate metabolism in carbon metabolism in HEB_vs_HEA or HEC_vs_HEA showed the apparent upregulation of A4695 (MLS1), A6232 (MAE1), PCK1 (A5260) in the glyoxylate cycle modules (M00012) and the apparent upregulation of A8906 (PGK1) in the glycolysis module (M00001) (Fig. 5F, Additional file 14). Together with the enrichment of carbon metabolism using the significantly upregulated expressed mRNA in HEB_vs_HEA (Additional file 1) or HEC_vs_HEA (Additional file 2), these observations confirmed the upregulation activities of the pathway modules of carbohydrate metabolism. The two modules (M00012 and M00001) were linked together, leading to the production of glucose-6P, which meant that the upregulated activity of the two modules could promote the production of glucose-6P (Fig. 5F) and further provided the intermediates for polysaccharide synthesis.
Multi-omics analysis of the hypothesized mushroom polysaccharides production biosynthetic pathways
Twenty homologous genes in H. erinaceus were obtained using Blastp (1e-5) of the sequences in yeast based on the hypothesized mushroom polysaccharides biosynthetic pathways (MPBP) according to reference[14]. Heatmap analysis of the mRNA genes involved in MPBP showed a noticeable expressed difference between HEA and the two mutated strains (Fig. 6A), especially the upregulated cluster marked by purple rectangular in HEB or HEC. Among the cluster, FBP1, UGDH, GAL10, and UXS1 had the prominent upregulation mRNA expression (Fig. 6B). Only two differentially expressed proteins involved in the MPBP occurred in HEB_vs_HEA (A0648) and HEC_vs_HEA (A6180). The two genes both belonged to GAL10 (UDP-glucose-4-epimerase) involved in the synthesis of polysaccharide repeat units and also had the upregulated expression in mRNA and protein based on the transcriptome and proteomics data (Fig. 6B). These upregulated genes involved in MPBP participated in the synthesis of polysaccharides repeat units (Fig. 6C), which meant the metabolic pathway of polysaccharide in the mutated strains was increased.
The twenty homologous genes in MPBP were enriched in amino sugar and nucleotide sugar metabolism, glycolysis/gluconeogenesis, galactose metabolism, fructose, and mannose metabolism (Fig. 6D). These results agreed partially with the 18 co-upregulated proteins in HEB_vs_HEA and HEC_vs_HEA in Fig. 5E, which conformed further the MPBP activation occurred in mutant strains HEB and HEC.
The DEGs from transcriptome or proteomics were used for the enrichment analysis with STRING and DAVID to obtain the co-enriched and differentially co-expressed (CDC) genes. The CDC genes identified by the Venn analysis of HEB_vs_HEA and HEC_vs_HEA (Additional file 15) were enriched in the KEGG pathways correlated with carbohydrates metabolism, such as carbon metabolism, pyruvate metabolism, glycolysis/gluconeogenesis (Fig. 6E). This result agreed with the multi-omics analysis of the MPBP in H. erinaceus (Fig. 5D). We did the qPCR and PRM analysis to verify the expressions of the CDC genes further. The list of the primers and the peptides followed by PRM was given in Additional file 16 and Additional file 17, respectively. qPCR and PRM analysis confirmed that most of the CDC genes shared a similar expression trend with the transcriptome and proteome data (Fig. 6F), such as the upregulation of MAE1(A6232), MLS1(A4695), PCK1(A5260) in the glyoxylate cycle modules (M00012) and A8906 (PGK1) in the glycolysis module (M00001) (Fig. 5F and Additional file 14). The CDC genes shared many KEGG pathways involved in carbohydrates metabolism (Fig. 6E), and therefore, their upregulation could increase the activities of carbohydrates metabolism for providing repeat units in polysaccharides production.
Multi-omics analysis of the glucose signaling regulation dysfunction associated with β-glucan production
Considering the 92% glucose ratio in the monosaccharide composition and a β-configuration glycosidic bonds in the new polysaccharide fraction produced from HEC (Fig. 1E and F), we further did the multi-omics analysis of the β-glucan biosynthetic process (GBP) to explore the regulatory pathway of glucan production. One hundred forty-six homologous genes in H. erinaceus were obtained using Blastp (1e-5) of the sequences in the yeast GBP (GO: 0051274). Heatmap analysis of expressed mRNA genes involved in the GBP showed a noticeable expressed difference between HEA and two mutated strains (Additional file 18A). Venn analysis showed that HEB_vs_HEA had 41 DEGs involved in the GBP bigger than 14 DEGs in HEC_vs_HEA (Fig. 7A). These DEGs were enriched in the MAPK signaling pathway, longevity regulating pathway, cell cycle, meiosis, and AGE-RAGE signaling pathway in diabetic complications, which were closely associated to signal transduction pathway (Fig. 7B). The pathway of longevity regulation and MAPK signaling pathway were also found in the enrichment analysis using the significantly co-downregulated proteins in HEC_vs_HEA or HEB_vs_HEA (Fig. 5D). Multi-omics analysis suggested that the two pathways played a regulatory role in the GBP of H. erinaceus.
Venn analysis of proteomics data showed that only five DEGs involved in the GBP occurred in HEB_vs_HEA or HEC_vs_HEA (Additional file 18B). STRING enrichment of five DEGs showed that they were enriched in signal transduction (GO:0007165). The 13 DEGs involved in the GBP that occurred in HEB_vs_HEA and HEC_vs_HEA were listed in Additional file 18C. Among them, the transcriptome data and proteomics data confirmed the two downregulated genes, A8173 and A5435 (Additional file 18C). STRING annotation showed that A5435 (Sch9) was serine/threonine-protein kinase involved in ribosome biogenesis, translation initiation, and the regulation of G1 progression. The previous report indicated that Sch9 might converge with the Ras-cAMP pathway downstream of PKA on its effector Rim15 [25]. Like TOR complex 1, Sch9 was required for cytoplasmic retention of Rim15 during exponential growth on glucose-containing medium [26], indicating that Sch9 also mediated the cell growth response to glucose [27].
The correlation analysis of 122 multiple omics data showed that HEB and HEC had a stronger negative correlation with HEA in the proteomic data than in the transcriptome data (Fig. 7C), which represented the apparent difference in protein expression in HEA and the mutated strains. The cluster consisting of steroid biosynthesis, longevity regulation, MAPK signaling pathway, and peroxisome reflected the apparent change in protein expression in HEB_vs_HEA or HEC vs_HEA (Fig. 5D). Multi-omics analysis of GBP confirmed that longevity regulation was involved in the GBP (Fig. 7B), shifting our attention to this pathway. Interestingly, the KEGG pathway showed that the genes, such as Ras2 (A5450), Cyr1 (A4603), PKA (A5435), and Msn 2/4 (A6938), in the pathway of longevity regulation-yeast (Additional file 19) shared partially with those in the pathway of meiosis (Additional file 20). Venn analysis showed that 12 genes were present in longevity regulating pathway and meiosis-yeast (Fig. 7D), and most of them belonged to the genes of the RAS-cAMP-PKA pathway (Fig. 7E).
One hundred fifty homologous genes of 12 sequences in H. erinaceus were obtained using Blastp (1e-5). Heatmap analysis confirmed a noticeable decrease in the protein expression profile of the RAS-cAMP-PKA pathway in HEC compared to the HEA and HEB (Fig. 7F). As the RAS-cAMP-PKA pathway acted as a glucose signals pathway [26], the decrease in the protein expression of the RAS-cAMP-PKA pathway could lose glucose sensing in H. erinaceus. Considering the highly increased rate of 47.45% of polysaccharide content in HEC compared to 23.25% in HEB (Fig. 1B), the downregulated activity of the RAS-cAMP-PKA pathway might be the leading cause of high β-glucan production.
A putative model of high polysaccharide production in the mutated H. erinaceus
The Ras-cAMP-PKA pathway plays a prominent role in responding to glucose availability and initiating the signaling processes that promote cell growth and division [28]. Considering that the Ras-cAMP-PKA pathway eventually connects to S phase progression in the KEGG pathway of meiosis (Additional file 20), the RAS-cAMP-PKA pathway blocking might lead to the inhibition of S phase progression. Exogenous CDK1/4/9 P276-00 inhibitor experiment was carried out, and the results showed that CDK1/4/9 P276-00 increased the mycelia biomass and intracellular polysaccharide content in H. erinaceus (Fig. 8A), especially at 10 µM, which further confirmed the role of RAS-cAMP-PKA pathway in polysaccharide synthesis. Meanwhile, β-glucan content was also increased in a concentration-dependent manner (Fig. 8B), which consistent with the results of the increased glucose proportion in the monosaccharide composition and production of new β-glucan fraction in the mutant strain HEC (Fig. 1E and Fig. 1F). These observations indicated that the inhibition of S phase progression could lead to the high polysaccharide production in the mutated H. erinaceus.
Based on the above data, we proposed a putative model of high polysaccharide production in the mutated H. erinaceus (Fig. 9). This model described a disordered process of polysaccharide anabolism. The decreased expression of the RAS-cAMP-PKA pathway might cause a long delay in the glucose response, lose glucose-sensing signaling response, finally inhibiting S phase progression. The block of S phase progression might induce the activities of pathways involved in polysaccharide production, such as the elevating activity of the modules of M00012 and M00001 (Fig. 5F and Additional file 14), and the pathways involved in further synthesis of repeating units from glucose-6P (Fig. 6C). The down-regulated ribosomal protein reduced the activity of protein translation and might interact with the block of S phase progression. Previous researches indicated the role of Sch9 in the Ras-cAMP pathway downstream of PKA through its effector Rim15 [25]. Our observations suggested the putative role of Sch9 in the reduced translation initiation as well as the block of S phase progression in Fig. 9.
The conserved Ras-cAMP-PKA pathway played a central role in the regulation of many biological aspects in eukaryotic organisms [29]. For example, this pathway governs pathogenesis, morphological transitions, nutrient sensing and acquisition, sexual reproduction, and stress responses in fungi [29–32]. In S. cerevisiae, the transcriptional responses to glucose are triggered by a variety of PKA mediated pathways, alone or in combination [33]. Multi-omics analysis of our study indicated the disfunction of the RAS-cAMP-PKA pathway in the ARTP mutated strain plays an important role in the stimulating of polysaccharide, especially β-glucan production.