Overview of transcriptome and small RNA sequencing. After transcriptome sequencing and filtering, approximately 63.92 Gb of clean reads were obtained, with an average of 6.54 Gb for each sample from nine cDNA libraries, for which the Q30 base percentage was greater than 93.54% (Table S1). The clean reads from each sample were compared with the specified reference genome, for which the efficiency of alignment varied from 91.52% to 97.59%. The above results indicate that the sequencing obtained in this study were of good quality and could be used for subsequent analysis. A total of 9921 transcripts were obtained, including 1008 that mapped to predicted new genes. Gene function annotations showed that 8851 genes had significant matches in COG, GO, KEGG, KOG, Pfam, SwissPort, eggNOG, or NR databases, with 3241 (36.62%), 6129 (69.25%), 3359 (37.95%), 4490 (50.73%), 6361 (71.87%), 5298 (59.86%), 7962 (89.96%), and 8837 (99.84%) genes, respectively. In addition to the functional annotation of O. sinensis, 11.77% of the genes also had high homology with Hirsutella minnesotensis (Figure 2A).
After removing adaptors and low-quality reads, 142.96 M clean reads of small RNAs were generated, with each sample yielding greater than 11.05 M. The statistical results are shown in Table S2 with an overview of small RNA classification and annotation. The normalized clean reads were used for analysis of small RNA distribution; the length distribution map of small RNA sequences demonstrated that the length of these small RNAs was 15–35 nucleotides (nt) (Figure 2B). In general, most of the clean reads were 23–26nt in length, with reads of 25nt being the highest.
DEGs and DEMs expression analysis of O. sinensis at differential development stages. To investigate the changes in gene expression levels in O. sinensis during the development of fruiting bodies, DESeq2 software was used to compare the gene expression of samples at different stages based on clean reads. In the three comparison groups, we identified a total of 2875 DEGs. The initial stage, MC, represented a control condition; the numbers of DEGs in the ST and FB periods were 977 and 1658, respectively. 1854 DEGs were also screened between the ST and FB stages (Figure 3A). There were only 157 co-expressed genes in all three stages, and the most significant gene changes occurred during the FB stage (Figure 3B). These differential genes included Cytochrome P450 monooxygenase (gene-G6O67_005633), Catalase (gene-G6O67_006909), Glucokinase (gene-G6O67_001528) and Phosphoenolpyruvate carboxykinase (gene-G6O67_008067) (Table S3), which are key enzyme genes in many metabolic pathways.
To investigate the known and putatively novel miRNAs expressed at the three stages of O. sinensis, we first compared the known mature miRNAs and miRNA precursors in miRBase; no conserved miRNAs were identified. However, a total of 106 novel milRNAs were identified in the nine small RNA libraries using the miRDeep2 program (Table S4). A differential expression analysis of the miRNAs between these three samples was performed based on normalized read counts (TPM) for each identified miRNA. We obtained 27, 48, and 57 differentially expressed milRNAs (DEMs) in MC vs ST, ST vs FB, and MC vs FB comparisons, respectively (Figure 3C). More DEMs were downregulated during FB development. Only 12 DEMs were co-expressed in all three stages. Characterizing the differential expression of miRNAs is important in predicting the occurrence and development of fruiting bodies in O. sinensis (Figure 3D).
Functional annotation and classification of DEGs. To infer the biological functions affected by DEGs at the three stages (MC, ST, and FB), we performed GO functional analysis. In the two developmental processes, 477 and 1027 DEGs were classified into 47 terms of three major biological processes (biological processes, cellular components, and molecular functions), respectively. The enrichment results of the three major biological functions of GO are shown in Figure 4A-B (P-value ≤ 0.03). The most dominant terms included the oxidation-reduction process, integral components of membranes, oxidoreductase activity, monooxygenase activity, and iron ion binding. KEGG pathway enrichment analysis was also conducted; 97 and 220 DEGs were enriched, corresponding to 77 and 103 pathways, respectively. In the process of ossification of O. sinensis, “Starch and sucrose metabolism”, “Tryptophan metabolism”, “Tyrosine metabolism”, and “Sphingolipid metabolism” pathways were significantly enriched (Figure 4C). In the FB formation stage, the degree of enrichment of “Biosynthesis of antibiotics” and “Carbon metabolism” varied greatly (Figure 4D). These metabolic pathways may closely relate to the formation of sclerotia and fruiting bodies. All DEGs, as well as GO and KEGG analysis results are shown in Table S3. However, some DEGs encoded functionally unknown proteins, which might relate to O. sinensis growth and development; further studies will be required to verify their functionalities.
Analysis of GSEA significant enrichment. While GO and KEGG use functional annotation to analyze DEGs from functional categories, effective statistical methods were not used to analyze the overall trend of gene expression. Therefore, GSEA was adapted to analyze the enrichment of genes differentially expressed in each GO term and KEGG pathway. In this study, normalized enrichment scores were used to draw a cluster heat map (P-value < 0.05) (Figure 5). In the process of sclerotia formation, phosphorylation-related signaling (Phosphorelay signal transduction system (GO:0000160), signal transduction by protein phosphorylation (GO:0023014), phosphorelay sensor kinase activity (GO:0023014), and oxidative phosphorylation (ko00190)) and oxidation-related (response to oxidative stress (GO:0000155), cellular oxidant detoxification (GO:0098869), peroxidase activity (GO:0004601), and monooxygenase activity (GO:0004497)) signaling were significantly upregulated. During the FB growth stage, the expression trend of ribosome-related terms and pathways were increased significantly, including the nucleolus (GO:0005730), ribosome (GO:0005840, ko030100), 90S preribosome (GO:0030686), and ribosome biogenesis in eukaryotes (ko03008). However, carbohydrate transport (GO:0008643), fatty acid degradation (ko00071), glyoxylate and dicarboxylate metabolism (ko00630), and carbon metabolism related to energy metabolism were downregulated.
Some prominent pathways are listed in Table 1, such the MAPK signaling pathway. DEGs in this pathway include mitogen-activated protein kinase kinase (Pbs2), glycerol-3-phosphate dehydrogenase (GPDH), catalase (CAT), and other important genes encoding enzymes. The expression of CAT was upregulated with log2(fold change) of 2.312 and downregulated with log2(fold change) of 2.160, respectively. In the citrate cycle, the genes encoding the enzymes malate dehydrogenase (MDH), phosphoenolpyruvate carboxykinase (PEPCK), and succinyl-CoA synthetase beta chain (LSC2), which catalyze the oxidation of citric acid for energy, were highest in the ST stage, upregulated with log2(FC) of 2.237, 3.607, and 3.025, respectively, compared with the FB stage, and were slightly higher than in the MC stage.
Gene ID
|
Annotation
|
FPKM (Average value)
|
MAPK-related genes
|
MC
|
ST
|
FB
|
gene-G6O67_001149
|
ATF/CREB family transcription factor
|
159.450
|
296.989
|
310.307
|
gene-G6O67_001860
|
mitogen-activated protein kinase kinase (Pbs2)
|
38.493
|
72.839
|
117.395
|
gene-G6O67_006243
|
histidine kinase (Sln1)
|
35.787
|
77.025
|
31.850
|
gene-G6O67_006845
|
STE3
|
10.897
|
33.218
|
15.188
|
gene-G6O67_005494
|
glycerol-3-phosphate dehydrogenase (NAD+)
|
105.461
|
88.457
|
16.160*
|
gene-G6O67_007807
|
ras homolog gene family
|
7.340
|
12.113
|
5.378*
|
gene-G6O67_001713
|
Catalase (CAT)
|
290.274*
|
1624.276
|
362.997*
|
Amino acid-related genes
|
gene-G6O67_003346
|
tyrosinase
|
29.075*
|
427.065
|
1.401*
|
gene-G6O67_004157
|
acetamidase
|
36.780*
|
94.229
|
23.372*
|
gene-G6O67_001432
|
amidase signature domain protein
|
74.558
|
36.214
|
4.771*
|
gene-G6O67_006218
|
amidase
|
84.452
|
101.630
|
4.214*
|
gene-G6O67_000733
|
acetyl-CoA acetyltransferase
|
153.290
|
166.490
|
50.872*
|
gene-G6O67_001178
|
glutaryl-CoA dehydrogenase
|
57.967
|
66.856
|
33.046
|
gene-G6O67_002860
|
aldehyde dehydrogenase family
|
645.414*
|
1378.017
|
351.529*
|
gene-G6O67_005633
|
cytochrome P450 monooxygenase
|
88.415
|
163.465
|
2.668*
|
gene-G6O67_001713
|
catalase
|
290.274*
|
1624.276
|
362.997*
|
gene-G6O67_006145
|
catalase-peroxidase
|
251.388
|
239.481
|
92.286*
|
Energy-related metabolism genes
|
gene-G6O67_000383
|
isocitrate dehydrogenase
|
87.656
|
48.958
|
102.003*
|
gene-G6O67_003773
|
isocitrate lyase
|
42.254*
|
67.603
|
13.513*
|
gene-G6O67_001528
|
glucokinase
|
17.177
|
11.438
|
176.097*
|
gene-G6O67_001884
|
transaldolase
|
259.266
|
281.600
|
113.116*
|
gene-G6O67_003014
|
Acyl-CoA dehydrogenase
|
73.135
|
72.973
|
37.604
|
gene-G6O67_003548
|
citrate synthase
|
326.182
|
61.291
|
0.000
|
gene-G6O67_003742
|
NAD(P)-binding domain protein
|
140.674
|
148.164
|
67.100
|
gene-G6O67_006486
|
malate dehydrogenase
|
788.497
|
1088.906
|
236.372*
|
gene-G6O67_006510
|
isocitrate dehydrogenase NADP
|
63.750
|
30.549
|
72.962
|
gene-G6O67_007408
|
succinate dehydrogenase cytochrome b small subunit
|
269.799
|
242.459
|
111.577
|
gene-G6O67_007964
|
acyl-CoA dehydrogenase domain protein
|
77.659
|
45.968
|
17.089
|
gene-G6O67_008750
|
succinyl-CoA synthetase beta chain
|
87.758*
|
209.074
|
26.275*
|
gene-G6O67_008752
|
alpha subunit of GDP-forming succinate-CoA ligase
|
118.068
|
145.249
|
42.014*
|
gene-G6O67_002224
|
succinate dehydrogenase iron-sulfur protein
|
281.527
|
115.061
|
248.025
|
gene-G6O67_004996
|
gluconolactonase precursor
|
16.812
|
5.929
|
8.262*
|
gene-G6O67_008325
|
fructose-1,6-bisphosphatase I
|
140.215
|
107.728
|
8.231
|
gene-G6O67_008067
|
phosphoenolpyruvate carboxykinase
|
929.438
|
1282.465
|
96.630*
|
Table 1. DEGs involved in metabolic pathways. Significant differences in expression were evaluated using one-way ANOVA (ST was used as the control group, *P < 0.01).
Integrated analysis of DEGs and DEMs. To explore the regulatory relationship between milRNAs and mRNAs, 1096 potential target genes of the milRNAs were predicted, with 112 target genes obtained from the 33 DEMs in MC vs ST, and 456 target genes from the 27 DEMs in ST vs FB. In order to understand the functions of these genes targeted by DEMs, GO annotation, and KEGG enrichment were performed. Target genes were classified into cell cycle-related, cyanoamino acid metabolism, and energy metabolism related pathways (Figure 6A-B). These results indicated that milRNAs played important roles in the growth process of O. sinensis.
Using ST as a control, there were 38 and 75 DEM-DEG relationship pairs found in MC and FB, respectively, based on the above mRNA and milRNA differential analysis results (Table S5). The network regulation diagram drawn by Cytoscape of some functionally annotated target genes indicated that one DEM could regulate more than one DEG, with both positive and negative correlation. Most milRNAs had more than one possible target gene, while different milRNAs could also regulate the same targets. As miRNAs regulate gene expression mainly by promoting cleavage of the target mRNAs or regulating TFs, we focused on negatively correlated pairs. According to the target regulation map in Figure 6C-D, key enzyme genes in β oxidation gene-G6O67_007081 (3-hydroxyacyl-CoA dehydrogenase, targeted by n_os_milR90) and ecological adapting-related gene gene-G6O67_007081 (tat pathway signal sequence, targeted by n_os_milR16) were upregulated. From the ST to FB stage, gene-G6O67_006617 (ABC transporter) and gene-G6O67_008466 (SET domain protein) were significantly downregulated by n_os_milR34, with a log2(fold change) of 5.106 and 3.096, respectively. According to the target gene annotation and regulatory network, n_os_milR16, n_os_milR21, n_os_milR34, and n_os_milR90 represent candidate milRNAs to affect FB development.
Validation of the DEGs and DEMs by RT-qPCR. To confirm the reliability of the sequencing data, a total of eight DEGs and four DEMs were randomly selected to validate the RNA-Seq and small RNA expression profiles. As expected, qRT-PCR results showed that most of these mRNAs and miRNAs shared similar expression with those from the sequencing data. Pearson correlation also showed that most of the relative expression levels were strongly correlated with FPKM/TPM (Figure 7), which indicates that our data was accurate and reliable in the subsequent analyses.