Transcriptome Analysis of Ophiocordyceps Sinensis Under Low-Temperature Treatment

This study aims to investigate the effect of low temperature on the gene expression in Ophiocordyceps Sinensis. To set up the experiments, two groups of Ophiocordyceps Sinensis were grown at 4°C and 18°C, respectively, for 60 days. Differential gene expression was analyzed by transcriptome sequencing and the morphology of mycelium was studied under scanning electron microscope for the two groups. The results showed that there were signicant differences in hypha morphology between the two groups. Transcriptome analysis indicated that 2799 genes were differentially expressed in the low-temperature treatment group compared to the normal temperature group. Among these genes, 1489 were up-regulated, whereas 1310 were down-regulated. Further analysis revealed that 33 genes were only expressed at low temperature compared to 26 genes which were only expressed at normal temperature. GO enrichment analysis showed that the 33 genes were related to redox enzymes, which have functions in clearing reactive oxygen species and free radicals caused by low-temperature environment. KEGG enrichment analysis showed that the 33 genes were mainly involved in the cyanoamino acid metabolism, biosynthesis of amino acids, carbon metabolism, pantothenate and CoA biosynthesis, biotin metabolism, phenylpropanoid biosynthesis, apoptosis and other metabolic pathways. Our observations provide fundamental clues on the relationship between the low-temperature stimulation and the fruiting body development and further promote the research progress on Ophiocordyceps Sinensis.


Introduction
The Chinese caterpillar fungus is a parasitic complex composed of the remains of the caterpillar and the fungal sexual stroma (Lo et al., 2013;Zhou et al., 2014;Liu et al., 2015). It is only distributed in the Qinghai-Tibet Plateau and its surrounding high-altitude areas, and has been used in traditional Chinese herbal medicine (Shrestha et al., 2010). Due to the high demand in markets, it has been placed at risk of extinction due to overharvesting. According to the latest IUCN Red List of Threatened Species released in 2020, O. Sinensis has been listed as a 'Near Threatened' species (Yang, 2020).
It is widely accepted that O. Sinensis is the only asexual type of Chinese caterpillar fungus (Sung et  Generally, the transformation from vegetative growth to reproductive growth of fungi needs low temperature for stimulation. It has been shown that the primordium of Lentinula edodes culture only differentiates under appropriate low temperature, whereas it cannot be formed under constant temperature (Feng et al., 2005). The low-temperature stimulation triggers the initiation of the sexual development of fungi, makes some speci c genes differentially expressed, and induces the start of sexual development and the fruiting body development. For instance, the fruiting body of Flammulina velutips could grow under darkness and low-temperature conditions due to the expression of Pf1 and Pf3 (Sakamoto, 2010). One study showed that in Pleurotus ostreatus, the transcription levels of the fourth chitinase gene of immature Pleurotus ostreatus fruiting body were higher than that in the mature stage during development, and the transcription levels were signi cantly increased after low-temperature treatment (Nishihara et al., 2007). In Flammulina velutips, the result also showed that the JRL1 gene played signi cant roles for fruiting body formation under low temperature-treatment (Lu et al., 2016).
It has been reported that the initiation of the O. Sinensis sexual development may be related to the unique environmental factors in the Qinghai-Tibet Plateau ecosystem . O. Sinensis generally goes through the frozen time for 5 months. Only in late April that is after the end of the frozen time, the fruiting bodies start to grow. According to this phenomenon, our experiments were designed under lowtemperature and normal-temperature treatments (control group) and differential gene expression was studied under these conditions for O. Sinensis by transcriptome analysis. Our results provide fundamental clues on the relationship between the low-temperature stimulation and the fruiting body development and promote the research progress on O. Sinensis.

Experimental strains and temperature treatment
Using the GZC: (Fungal Name registration. FN570054), O. Sinensis strain was isolated from the Chinese caterpillar fungus tissues and cultivated under arti cial conditions to get its fruiting body. The strain is stored in the Biological Seedling Engineering Institute of Lanzhou Jiaotong University. After 60 days of vegetative growth at 18℃, the strains were randomly divided into two groups: the control group (Group B) incubated at 18℃ and the experimental group incubated (group C) at 4℃. After 60 days, the samples from the two groups were collected and stored at -80℃ until use. 2) After xation, the samples were rinsed with 0.1 mol/L phosphate buffer for 3 times, each for 15 min; 3) Samples dehydration: the samples were treated with different concentrations of ethanol with 15 min for each concentration. After that, the samples were treated twice with pure alcohol, each time for 20 min; 4) The samples were treated with alcohol/isoamyl acetate (volume ratio: 1:1) for 30 min, followed by 100% isoamyl acetate for 1-2h; 5) The samples were freeze-dried for 15h. 6) At the end of drying, the surface was sprayed with gold by an ion sputtering apparatus and then observed by the scanning electron microscope.
The prepared samples were observed under scanning electron microscope by colleagues from the Life Science Research and Experiment Center of Lanzhou University and the State Key Laboratory of Solid Lubrication, Lanzhou Institute of Chemical Physics, Chinese Academy of Sciences.

Total RNA extraction and mRNA separation
The total RNA was extracted from the samples with the modi ed CTAB method. NanoDrop ultraviolet spectrophotometer and Agilent2100 were used for quality check and quanti cation. After qualify detection, the mRNA was separated by using the magnetic bead for subsequent experiments.

cDNA library construction and sequencing
The library construction and sequencing work was completed by Shanghai OE Biotech. Co., Ltd. These libraries were sequenced on the Illumina sequencing platform (Illumina HiSeq X Ten) and 150 bp pairedend reads were generated.

Sequencing data analysis
Raw sequencing data was ltered to obtain high-quality Clean Data. Trinity (Grabherr et al., 2011) software was used to splicing Clean Reads according to the no-reference genome analysis strategy, and the transcript was obtained for subsequent analysis. The reason for the non-reference genome analysis is that the published genome data of O. Sinensis is controversial, and there is no widely recognized reference genome of O. Sinensis (Hu et al., 2013;Li et al., 2016;. After splicing, DESeq software (Anders and Huber, 2010) was used to analyze the differential expression of Unigene between sample groups. The screening condition was expression differential multiple |log2FoldChange|>1, signi cance P-value<0.05. GO enrichment analysis and KEGG enrichment analysis were performed for the differentially expressed genes obtained from screening to annotate and predict the functions and metabolic pathways of the related genes. Raw Illumina sequencing data of O. sinensis were submitted to NCBI as BioProject SUB8768026.
2.3.6 qRT-PCR analysis cDNA was synthesized from mRNA using Kit (Vazyme). Then the reaction was carried out on the LightCycler® 480 Real-time PCR Instrument using the QuantiFast® SYBR® Green PCR Master Mix (Qiagen, Germany). The results were analyzed by using the calculation method of 2-∆ Ct with 18SrRNA of O. sinensis as internal reference gene. The related primer sequences were analyzed by Roche LCPDS2 software and synthesized by Beijing Qingke Xinye Biotechnology Co., Ltd. Primer sequences are listed in  Figure 3 showed that the hypha morphology of O. Sinensis was a signi cant difference at different temperatures. At the magni cation 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 magni ed 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. 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. To investigate the mechanism of changes in the low temperature conditions for O. Sinensis, the transcriptome sequencing was carried on at rst, and then the original pieces been ltered , 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.

Quality assessment and analysis of sequencing data
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.

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. 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.  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 re ected in the cell membrane, organelle, cell part, etc. The molecular functions mainly include catalytic activity, binding protein, transport activity, and so on.

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 classi cation and the number, the gene was matched the corresponding eggNOG classi cation, and the results were shown in Figure 5. A total of 15,819 Unigene were annotated and classi ed 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 signi cant to investigate the functions of these genes further.

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 classi cations 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.

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 speci cally 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 speci cally 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 speci cally expressed genes. The function of these genes under the stimulation of low temperature should be investigated in future research in this eld. Statistical analysis was performed on 15 genes with known functions, and the results were shown in Table 6. 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 signi cant enrichment of the differential expressed genes, to demonstrate the gene functions that may exist in the samples. We rst 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 signi cant enrichment of differential genes. The GO classi cation 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 rst 10 categories are selected for drawing, and different colors represent different functional classes.
In the GO classi cation 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 classi cation 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 classi cation 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%.

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 signi cant 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 signi cant expression of the top 20 metabolic pathway enrichment results was shown in gure 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 signi cant, with 1 up-regulated gene and 3 down-regulated genes, among which the signi cantly up-regulated gene was gamma-glutamyl transpeptidase(ggt). Details are shown in Table 7, metabolic pathway annotation results are shown in S2. 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 signi cantly 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.

Transcriptome sequencing expression pattern veri cation 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 veri cation. The results are shown in Figure 11, and the veri cation 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 signi cance. Both the vernalization of plants and the low temperature stimulation of fungi can be regarded as a kind of adversity that low temperature makes organisms face. This kind of adversity will promote the transformation of the organism from the vegetative stage to the reproductive stage, which is usually closely related to the growth environment of the organism. The low temperature environment is exactly a kind of adversity that O. sinensis must face, so the appropriate low temperature treatment is likely to promote the O. sinensis to enter the sexual reproduction stage from the nutritional stage.

Functional analysis of unknown and known genes at low temperature
Statistical analysis of speci c differential expressed genes showed that 74.6% of speci c expressed genes were unknown genes, hypothetical protein genes and unknown functional protein genes, and the functions of these genes were not clear. In order to fully understand the hypothermia response mechanism of O. Sinensis, the function of this part of gene must be clearly understood, and the expression of this part of gene may be closely related to the initiation of sexual development of O. Sinensis. Therefore, It might be the breakthrough conten in the future on O. Sinensis research area. Analysis of speci c differential expressed genes with clear annotations showed that most of the differential expressed genes in the low-temperature treatment group were structural protein genes and a few were anabolic enzyme genes. On the contrary, most of the differential expressed genes in the normal temperature control group were anabolic enzyme genes, and a few were structural protein genes. Studies have found that the low temperature environment has a signi cant inhibitory effect on the metabolic activities of organisms (Kubien et al., 2003), and protein play an important role for organisms in the resistance to the low temperature environment (Doxey et al., 2006;Wang et al., 2012). This indicated that the anabolism of O. Sinensis was signi cantly inhibited in the low temperature environment, and other structural proteins were synthesized to adapt to the low temperature environment, enabling O. Sinensis to live through the long time in the low temperature environment.

Attributes and function of differential expressed genes
All the differential expressed genes were annotated by GO and functional cluster analysis was conducted to understand the functional classi cation of the differential expressed genes between groups. In the GO enrichment analysis between groups of B-vs-C, it was found that the differential expressed genes showed a GO Term of signi cant enrichment. Further study revealed that the differential gene coding products enriched in this functional classi cation all had redox activity. This may indirectly con rm that it is the low temperature that creates a kind of adversity, which urges O. Sinensis to enter the sexual stage to cope with this adversity. Studies have shown that the content of free radicals in plants increases signi cantly in cold environments (Hodges et al., 1997), which also leads to the expression of a large number of

What role does the gene ggt play?
KEGG is a powerful tool for metabolic analysis and metabolic regulatory network research in biological organisms, which can nd pathways with signi cant enrichment in differential expressed genes (Kanehisa et al., 2006). In this study, KEGG pathway enrichment analysis for B-vs-C showed that ggt gene was up-regulated under the in uence of temperature, and the amount of gene expression affected the content of gamma-glutamyl transpeptidase in the organism, thereby regulating the content of glutathione in the organism. Some studies have reported that the amount of glutathione may affect the development of fungal sexual structure (Thön et al., 2007). Therefore, low temperature treatment signi cantly upregulated the ggt gene expression in O. Sinensis, indicating that this gene may be a key gene for the response of O. Sinensis to low temperature environment and may be a key gene for its sexual development initiation. At present, ggt gene has been widely studied in animals. For example, in animal cell research, it has been found that ggt gene can affect mitochondrial structure maintenance, membrane integrity, cell differentiation and development (Zhang, 2005). However, there few study on fungi and plant cells until now, so intensive study research on the ggt gene in fungi would be a great signi cance. Other

Declarations
Author contributions FX designed research; QJS, LH and GZ performed experiments; FX, QH and LH analyzed data; QJS, LH, VU and YXX wrote the paper. All authors read and approved the nal manuscript.
Funding This research was funded by National Natural Science Foundation of China (31560003).

Declaration of Competing Interest
The authors declare that they have no known competing nancial interests or personal relationships that could have appeared to in uence the work reported in this paper.

Data availability statements
The data used to support the ndings of this study are available from the corresponding author upon