3.1 Natural ambient temperature of Chinese caterpillar fungus
Figure 1 shows the ambient temperature of the natural growth of Chinese caterpillar fungus recorded by the author in Meiren Grassland, near Hezuo city, Gansu Province, China. This temperature is the average temperature of each month in a year. As concerning to figure 1, from November to mid-March of every year, there is a frozen time of about 5 months, during these months Chinese caterpillar fungus is always in the corpse state. When the temperature began to rise in April, the fungus started to grow. In late April, it gradually broke through the ground and then entered the harvesting period. Studies have found that the fruiting body of O. Sinensis can only be formed after a period of low-temperature induction (Tu et al., 2010; Guo et al., 2015), which is very similar to vernalization of plants. However, the molecular or biochemical mechanisms that regulate this process are not clear until present.
3.2 Effect of low-temperature treatment on the colony and hypha morphology of O. Sinensis
Refer to figure 2, there were significant differences in colony morphology under different temperatures conditions. It is obvious signs of growth and development of O. Sinensis under two temperature treatments. The O. Sinensis is a hypothermic fungus, which resists and continue to develop under low temperature. Genomic investigation showed that O. Sinensis has antifreeze protein and can increase the accumulation of lipids and unsaturated fatty acids to adapt to an extremely cold environment (Hu et al., 2013). Besides, the Genome of O. Sinensis contains many functional genes, adapted to low temperature (Xia et al., 2017). These reasons enable the O. Sinensis to continue its slow growth and development at low temperatures. Compared to the low-temperature environment of 4℃, 18℃ is the optimal temperature of the O. Sinensis growth. Under this condition, the growth is better and more aerial mycelium.
Figure 3 showed that the hypha morphology of O. Sinensis was a significant difference at different temperatures. At the magnification of 500 times, showed that the hypha at the 4℃, low-temperature treatment group is slimmer than at the 18℃, normal temperature control group, and there is a lot of helically curved hypha as shown in the red circle, while there is almost no such hypha in the 18℃ normal temperature control group. During the fungi growth and development, hyphae will appear in the phenomenon of bending and kinking, which will form the primordium, and the appearance of the primordium is an important sign of fungal fruit body development. After low-temperature treatment, the helical hypha bending was looked like having a certain connection of hypha kink and primordium formation. When the mycelia were magnified 10,000 times, characteristic verrucous structure, which conidia produce, appeared on the mycelia surface, showed in the red circle in the control group at 18℃, while the mycelia surface appeared smooth and no similar structure in low-temperature treatment at 4℃. It’s suggested that no conidia produce under low-temperature treatment.
3.3 Quality control results after total RNA extraction
Table2 Test results of samples
Sample
|
Concentration(μg/μl)
|
A260/280
|
A260/230
|
Volume (μl)
|
Total(μg)
|
Agilent2100
|
28S/18S
|
RIN
|
B1
|
0.5774
|
2.21
|
2.19
|
15
|
8.66
|
2.1
|
9.9
|
B2
|
1.0999
|
2.19
|
2.10
|
15
|
16.50
|
2.0
|
10.0
|
B3
|
0.6040
|
2.22
|
2.16
|
25
|
15.10
|
2.0
|
9.9
|
C1
|
0.4716
|
2.17
|
2.13
|
25
|
11.79
|
1.6
|
9.0
|
C2
|
0.6310
|
2.20
|
2.11
|
25
|
15.78
|
1.4
|
8.6
|
C3
|
0.3788
|
2.19
|
2.15
|
25
|
9.47
|
2.0
|
10.0
|
The quality control results after total RNA extraction are shown in Table 2. The quality of RNA samples for transcriptome sequencing requires higher nucleic acid purity when the OD260/OD230 is greater than 2, and RNA Integrity Number (RIN) value is between 1 and 10. The closer the value is to 10, the better the Integrity (Bolger et al., 2014). Refer to Table 1, the RNA purity, concentration, integrity, and other quality indicators of the 6 samples fully meet the requirements of subsequent transcriptome sequencing, so the cDNA library can be constructed for sequencing.
3.4 Quality assessment and analysis of sequencing data
Table 3 List of quality preprocessing results for sequencing data
Sample
|
Clean_reads
|
Clean_bases
|
Valid_bases
|
Q30
|
GC
|
B1
|
47847626
|
6930877946
|
93.73%
|
96.01%
|
60.36%
|
B2
|
47938792
|
6945590540
|
93.18%
|
96.31%
|
60.20%
|
B3
|
47970638
|
6960217947
|
93.84%
|
95.96%
|
60.45%
|
C1
|
47535504
|
6906826837
|
94.13%
|
96.23%
|
59.95%
|
C2
|
48020666
|
6964636068
|
93.78%
|
96.03%
|
60.14%
|
C3
|
47808510
|
6938280187
|
93.85%
|
96.00%
|
60.24%
|
To investigate the mechanism of changes in the low temperature conditions for O. Sinensis, the transcriptome sequencing was carried on at first, and then the original pieces been filtered , the data obtained was list in table 3. The group B (18 ℃ temperature control group, the same below) and group C (4 ℃ low-temperature treatment group, the same below) of Q30 average 96.09% and 96.08% respectively. Additionally, GC relative contents of group B was 60.34% and group C was 60.11%. Therefore, transcriptome sequencing data is of good quality and high accuracy, which can be used for subsequent analysis.
3.5 Splicing, assembly and length analysis of transcriptome
Transcript and Unigene sequence are counted, and the number of Unigene obtained after splicing was 37601, the maximum length of the sequence was 13572bp, and the N50 is 1524bp, as shown in the following Table 4.
Table 4 General statistics of splicing situation
|
Contig
|
Transcript
|
Unigene
|
Total Length (bp)
|
38531639
|
67412677
|
29820756
|
Sequence Number
|
99748
|
62682
|
37601
|
Max. Length(bp)
|
15695
|
13572
|
13572
|
Mean Length (bp)
|
386.3
|
1075.5
|
793.1
|
N50 (bp)
|
1239
|
2031
|
1524
|
N50 Sequence No
|
7372
|
10440
|
5658
|
N90 (bp)
|
133
|
394
|
274
|
N90 Sequence No
|
66313
|
37611
|
24866
|
GC%
|
56.79
|
58.71
|
58.22
|
3.6 Transcriptome functional annotation
3.6.1 Summary of annotation results
The assembly transcripts of O. Sinensis were compared in NR, GO, KEGG, eggNOG, and Swissprot databases respectively, and the statistical results were shown in Table 5.
Table 5 Unigene annotation results and ratios in each database
Annotation in Database
|
Isoform number
|
Percentage (%)
|
NR
|
18294
|
48.65
|
GO
|
11583
|
30.81
|
KEGG
|
3489
|
9.28
|
eggNOG
|
15819
|
42.07
|
Swissprot
|
14260
|
37.92
|
In all database
|
3112
|
8.2
|
3.6.2 GO annotation statistics
Using the BLAST2GO software (the default parameters) for GO annotation statistics. After the annotation, the result was completed corresponds to the GO Term and counted the number of Unigene annotated to the secondary classification. The results were shown in Figure 4. The number of Unigene annotated in GO is 11,583, accounting for 30.81% of all Unigene. Most of the Unigene successfully annotated here concentrates on the metabolic process, cellular process, single biological process, etc. in the biological process. The cellular component is mainly reflected in the cell membrane, organelle, cell part, etc. The molecular functions mainly include catalytic activity, binding protein, transport activity, and so on.
3.6.3 eggNOG annotation statistics
During the eggNOG annotation of Unigene, the eggNOG number of the optimal aligned result was assigned to the corresponding gene. Combining the relationship between eggNOG classification and the number, the gene was matched the corresponding eggNOG classification, and the results were shown in Figure 5. A total of 15,819 Unigene were annotated and classified into 25 functional categories of eggNOG, accounting for 42.07% of the total amount of Unigene. The most annotated genes were unknown genes functions, with 5,415. This may indicate that many Unigene spliced have not been studied yet. It is great significant to investigate the functions of these genes further.
3.6.4 KEGG annotated statistics
This part of the annotation was done automatically by KEGG's KAAS system. In the annotation process, the gene set selected the close relationship species, and BBH (bi-directional Best Hit) was the judgment principle of gene KO. After KO annotation was completed, the results were associated with the relevant metabolic pathways, as shown in Figure 6. A total of 3,468 Unigene were annotated into 35 metabolic pathways of 5 classifications in this pathway, the more typical of which are carbohydrate metabolism and amino acid metabolism in metabolism, translation in genetic information processing, signal transduction in environmental information processing, transport and catabolism in cellular processes, etc.
3.7 Differential gene analysis
In this experiment, the 18℃ normal temperature treatment group was used as the control group, and the gene expression difference of O. Sinensis was analyzed under 4℃ low temperature treatment (Figure 7). The results showed that the total number of differential expressed genes was 2799, among which the number of up-regulated genes was 1489, and the number of down-regulated genes was 1310. After further screening of differential expressed genes, 59 specifically expressed genes were obtained, including 33 differential expressed genes only at low temperature and 26 differential expressed genes that were not expressed at all at low temperature(S1), and cluster analysis was performed for these two types of genes (Figure 7C). The annotation results showed that 34 of these specifically expressed genes were unknown genes without any annotation, 9 were hypothetical protein genes, 1 was unknown functional protein gene, and 15 were known genes with annotation results. It can be seen that unknown genes, hypothetical protein genes and unknown functional protein genes account for 74.6% of the total number of specifically expressed genes. The function of these genes under the stimulation of low temperature should be investigated in future research in this field. Statistical analysis was performed on 15 genes with known functions, and the results were shown in Table 6.
Table 6 Specific expression genes with known functions
ID
|
Description
|
B (fpkm)
|
C (fpkm)
|
pvalue
|
TRINITY_DN9987_c0_g1
|
Chitinase 18-11
|
0
|
0.95
|
0.000765708
|
TRINITY_DN15069_c0_g5
|
Ribosomal protein L22/L17
|
0
|
6.94
|
0.001751776
|
TRINITY_DN15105_c0_g1
|
Gag polyprotein, partial
|
0
|
4.89
|
0.004255424
|
TRINITY_DN11747_c0_g1
|
Metalloprotease 1
|
0
|
0.34
|
0.017059201
|
TRINITY_DN25953_c0_g1
|
Matrix protein
|
0
|
3.22
|
0.021485746
|
TRINITY_DN14003_c1_g2
|
Triose-phosphate transporter
|
0
|
8.58
|
0.037477241
|
TRINITY_DN18536_c0_g1
|
COPII-coated vesicle protein
|
0
|
4.12
|
0.049795154
|
TRINITY_DN15548_c1_g1
|
Trichome differentiation protein GL1
|
2.8
|
0
|
8.94798E-15
|
TRINITY_DN6833_c0_g1
|
Acetamidase/Formamidase
|
0.55
|
0
|
0.001339126
|
TRINITY_DN22662_c0_g1
|
Polyketide synthase
|
2.92
|
0
|
0.003595873
|
TRINITY_DN11676_c0_g2
|
Tyrosyl-DNA phosphodiesterase 1
|
3.38
|
0
|
0.012752019
|
TRINITY_DN22850_c0_g1
|
Elastinolytic metalloproteinase Mep
|
5.7
|
0
|
0.017174642
|
TRINITY_DN14860_c0_g1
|
Glucose-methanol-choline oxidoreductase
|
1.84
|
0
|
0.022876803
|
TRINITY_DN19095_c0_g1
|
Lovastatin nonaketide synthase
|
5.22
|
0
|
0.023435381
|
TRINITY_DN1072_c0_g1
|
Beta-lactamase domain-containing protein
|
0.53
|
0
|
0.032627589
|
Analysis of 15 genes of known functions with specific expression revealed that 7 specific expression genes in the low temperature treatment group, among which there are 5 for structural protein gene, corresponding to ribosomal protein L22/L17, gag polyprotein, matrix protein, triose-phosphate transporter, COPII-coated vesicle protein and 2 for enzyme genes, respectively is chitinase 18-11 and metalloprotease 1. There are 8 functional genes known which specific expression in normal temperature group, among which there are 6 for synthetic metabolism related enzyme gene, relating to acetamidase/formamidase, polyketide synthase, tyrosyl-DNA phosphodiesterase 1, elastinolytic metalloproteinase Mep, glucose-methanol-choline oxidoreductase and lovastatin nonaketide synthase, there are 2 structural protein genes, trichome differentiation protein GL1 and beta-lactamase domain-containing protein respectively. It can be seen that the anabolism of O. Sinensis is blocked inside at low temperature, and a large number of structural proteins are expressed to prevent cell from damaged and adapt to the low temperature environment.
3.8 Functional annotation of differential expressed genes
3.8.1 GO annotation
The purpose of GO functional enrichment analysis is to obtain the GO functional entries with significant enrichment of the differential expressed genes, to demonstrate the gene functions that may exist in the samples. We first mapped all genes to each Term of the Gene Ontology database and calculated the number of differential genes in each Term. With the whole genome as the background, we used hypergeometric distribution to calculate the terms with significant enrichment of differential genes. The GO classification statistics of the intergroup differential genes in B-vs-C are shown in Figure 8. Divide the enrichment analysis results according to the three functional categories of molecular function, cellular component, and biological process, and order the results of each functional category from small to large according to Corrected P-value. The first 10 categories are selected for drawing, and different colors represent different functional classes.
In the GO classification statistics of the intergroup differential genes of B-vs-C, it was found that the number of the catalytic activity differential genes in the molecular functional classification was the largest, accounting for 58.5%. The proportion of differential genes related to the membrane was the highest among cellular components, which was 37.5%. Under the classification of biological processes, it can be seen that the number of differential genes in the oxidation-reduction process is the largest, accounting for 19.7%.
The maintenance of intracellular environmental balance is particularly important for the formation and normal development of fungal fruiting bodies. As the second messenger, reactive oxygen species (ROS) are involved in physiological processes such as cell signal response and cell proliferation, and the maintenance of its normal level and gradient distribution is very important for sexual reproduction and fruiting body development of fungi (Tudzynski et al., 2012). Studies have proved that temperature can regulate the activity of antioxidant enzymes in plants, thus affecting the tolerance of plants to temperature changes (Kang and Saltveit, 2002). The increase of ROS level in Cordyceps militaris cells is usually accompanied by decreasing the ability of fruiting body formation, and the degradation of frustum development in Cordyceps militaris, which can be inhibited by the expression of exogenous antioxidant enzyme-glutathione peroxidase A (Xiong et al., 2013). In the term of catalytic activity with the largest number of different genes in GO classification statistics, proteins encoded mainly include enzymes related to antioxidant reduction. This indicates that the low temperature also promotes the overexpression of antioxidant enzyme activity in O. Sinensis cells, thus eliminating the adverse effects of ROS and promoting the formation and development of fruiting bodies.
3.8.2 KEGG annotation analysis
Count the number of differential expressed genes in each KEGG Pathway at different levels, and then determine the metabolic pathways and signaling pathways in which the differential expressed genes are mainly involved. With the whole genome as the background, the hypergeometric distribution was used to calculate the pathway of significant enrichment of differential expressed genes. Enrichment analysis of B-vs-C showed that 238 differential expressed genes were successfully annotated into 123 pathways, selection of significant expression of the top 20 metabolic pathway enrichment results was shown in figure 9, mainly focused on the metabolism of carbon metabolism, biosynthesis of amino acids, nitrogen metabolism, cyanoamino acid metabolism, pantothenate and CoA biosynthesis, biotin metabolism, phenylpropanoid biosynthesis, apoptosis and other metabolic pathways. Among them, the expression of genes in the cyanoamino acid metabolism pathway was the most significant, with 1 up-regulated gene and 3 down-regulated genes, among which the significantly up-regulated gene was gamma-glutamyl transpeptidase(ggt). Details are shown in Table 7, metabolic pathway annotation results are shown in S2.
Table 7 The differential expressed genes in the cyanoamino acid metabolism pathway.
Gene ID
|
Gene name
|
Description
|
Regulate
|
log2FoldChange
|
TRINITY_DN10440
|
ggt
|
Gamma-glutamyl transpeptidase
|
up
|
1.0046
|
TRINITY_DN11153
|
fmda
|
Formaidase
|
down
|
-1.3560
|
TRINITY_DN15394
|
ghf3
|
Glycoside hydrolase, family 3
|
down
|
-1.4480
|
TRINITY_DN15752
|
cel3b
|
Beta-glucosidase
|
down
|
-2.4702
|
The gamma-glutamyl transpeptidase plays an important role in the metabolism of glutathione and amino acids, catalyzing the reaction of glutathione and amino acids to produce glutamine (Gln) and cysteamine glycine (Cys-Gly) (Figure 10). Studies have found (Xi et al., 2020) that reductive glutathione in Aspergillus nidulans can promote the production of its sexual structure. After low-temperature treatment, the expression of ggt gene of O. Sinensis was significantly up-regulated. This indicates that low-temperature treatment makes ggt gene expressed in large quantities to synthesize more gamma-glutamyl transpeptidase to regulate the content of glutathione inside, thus affecting the production of the sexual structure of O. Sinensis. Therefore, the ggt gene is likely to be a key gene for the response of O. Sinensis to the low-temperature environment and may be related to the initiation of fruit body development.
3.9 Transcriptome sequencing expression pattern verification was performed by RT-PCR
In this study, key differential expressed genes were obtained through transcriptomic sequencing technology. In order to further verify the accuracy of the transcriptomic sequencing data, 8 differential expressed genes were selected for Real-time PCR verification. The results are shown in Figure 11, and the verification results are consistent with the transcriptomic sequencing data shown in Table 8. This indicates that the transcriptome data obtained in this experiment have high reliability and biological significance.
Table 8 Transcriptome sequencing results of differential genes
Gene ID
|
Description
|
FPKM(B)
|
FPKM(C)
|
B-vs-C
|
TRINITY_DN10440_c0_g1
|
Gamma-glutamyltranspeptidase
|
8.43
|
18.64
|
Up
|
TRINITY_DN11153_c0_g1
|
Formamidase
|
145.15
|
59.89
|
Down
|
TRINITY_DN15752_c0_g2
|
Hypothetical protein
|
41.06
|
7.77
|
Down
|
TRINITY_DN8973_c0_g1
|
Sorbitol dehydrogenase
|
13.41
|
29.29
|
Up
|
TRINITY_DN15673_c0_g2
|
Nitrate reductase
|
25.9
|
95.54
|
Up
|
TRINITY_DN15159_c2_g2
|
NADP-specific glutamate dehydrogenase
|
38.35
|
369.31
|
Up
|
TRINITY_DN12710_c0_g2
|
Cyanate hydratase
|
46.74
|
23.55
|
Down
|
TRINITY_DN14699_c0_g2
|
Ribose 5-phosphate isomerase
|
52.38
|
296.29
|
Up
|