Combining transcriptomics and metabolomics to characterize ergosterol biosynthesis in Flammulina velutipes during the fruiting process

Background: Flammulina velutipes (F. velutipes) is one of the most important mushrooms in Japan and China because of its potential medicinal and nutritional value. Ergosterol is an important precursor of vitamin D2, progesterone, hydrocortisone, brassinolide, and novel anticancer and anti-HIV drugs. One of the main methods for obtaining ergosterol is extraction of its natural products from fungi. However, because of the low production of ergosterol due to its inefficient biosynthesis, its approach cannot meet human needs. Therefore, increasing the ergosterol level in fungi is an important goal to be pursued. Although genes encoding key enzymes in the biosynthesis pathway of ergosterol have been identified in Saccharomyces cerevisiae, the mechanism and regulation of ergosterol biosynthesis in F. velutipes remain unclear, and are the focus of this study. Results: In this work, nine cDNA libraries produced from the three stages of F. velutipes were sequenced by using the Illumina HiSeqTM 4000 platform, resulting in at least 6.63 Gb of clean reads from each library. De novo sequence assemblers were used to generate 220,523 transcripts and 28,330 unigenes, which were annotated according to seven protein databases. Here, we combined transcriptomics and metabolomics approaches to investigate the biosynthesis of bioactive ergosterol from F. velutipes. The contents of 16 intermediates or metabolites from different stages (young fruiting body stage and mycelium stage) were found to be significantly different. Then, transcriptional profiling of the genes involved in ergosterol biosynthesis was carried out, and the transcripts of key genes involved in the metabolism were analysed. A total of 51 key unigenes (12 upregulated unigenes and 39 downregulated unigenes) were identified as differentially expressed. Furthermore, four genes (i.e., ERG10s, ERG1s, ERG11s and ERG26s) were identified as the most important genes that played roles in the regulation of the pathway flux towards ergosterol synthesis in F. velutipes. genes were confirmed in the transcriptome of F. velutipes . The results revealed 19 DEGs related to ergosterol biosynthesis of F. velutipes at different developmental stages. After transcriptomic and metabolomic analysis of the ergosterol biosynthesis, four genes ( ERG1s , ERG10s , ERG11s , ERG26s ) were identified as biosynthesis. Transcriptome analysis of extracts of F. velutipes in three different developmental stages resulted in 28,330 unigenes. Comparing FrII against FrI group, we identified 51 (12 up-regulated and 39 down-regulated) unigenes that encoded enzymes related to ergosterol biosynthesis. The results showed that ERG1s , ERG10s , ERG11s and ERG26s were the most important four genes in the regulation of the pathway flux towards ergosterol biosynthesis in F. velutipes . We used RT-qPCR method to check the reliability of transcriptomics of F. velutipes . The profiles of sterols in F. velutipes from three different development stages were listed. In this article, we performed metabolomics with LC-MS approach to characterize the ergosterol biosynthesis of F. velutipes , and identified 17 important differential intermediates or metabolites related to ergosterol biosynthesis. Combining transcriptomics and metabolomics, we explored a regulation relationship between ergosterol biosynthesis genes and metabolites in F. velutipes . To sum up, this study is instrumental for further research of the biosynthesis of sterol-related intermediates or metabolites in F. velutipes

metabolomics, we explored the regulatory relationship between ergosterol biosynthesis genes and metabolites in F. velutipes. Conclusions: The data obtained in this work provide useful information for understanding the biosynthesis, metabolism and regulation of sterols in F. velutipes and for channelling the metabolic flux towards ergosterol at different growth stages.
Ergosterol is the most abundant sterol in F. velutipes , so its content could serve as an important indicator for assessing the quality of sterol-related products of F. velutipes [14,15].
Ergosterol (C 28 H 43 OH) is a typical sterol of fungi and it is an important constituent of various membrane structures of fungal cells. It displayed multifaceted physiological functions in cells, such as ensuring cell viability, membrane fluidity, membrane integrity, and cellular material transport [16]. Therefore, once ergosterol is lacking, it would cause abnormal cell membrane function and even cell rupture [17][18][19]. In recent years, a variety of fungicides, collectively called sterol biosynthesis inhibitors (SBIs), have been successfully developed to target some enzymes or end products of ergosterol biosynthesis pathway, and have been widely used in medicine and agricultural production [20].
Ergosterol level is also one of quality indicators for some traditional Chinese medicines and Chinese medication preparations. More importantly, ergosterol and some of its biosynthetic intermediates are important metabolites of great economic value. In the pharmaceutical industry, ergosterol is an important precursor for vitamin D2, progesterone, hydrocortisone, and brassinolide, and almost all steps of its biosynthetic process represent potential drugs [21][22][23][24][25].
Ergosterol and its derivatives are mainly obtained by the methods of chemical synthesis, genetic engineering and metabolic engineering [26]. Because of the various steps, long route, low efficiency and high cost, chemical synthesis of ergosterol and its derivatives is not the preferable way to obtain them. One of the main approaches to producing ergosterol and its derivatives is through metabolic engineering of yeast, but the common problem with the strains used in production is the low content of ergosterol in cells [27].
The biosynthesis of ergosterol is an extremely complicated process. Transcriptional regulation of expression of related genes is one of the main means to control ergosterol biosynthesis, orfeedback regulationcanplay an important role in regulating ergosterol production [28,29]. Sterol regulatory element-binding proteins (SREBPs) are transcriptional factors which bind to thesterol regulatory elementDNA sequence. Therefore, manipulation of biosynthesis genes by genetic engineering may be an effective way to change sterol biosynthesis and intracellular sterol components. Although metabolic engineering or genetic engineering has shed light on synthetic pathways in Saccharomyces cerevisiae, how the ergosterol biosynthesis genes regulate ergosterol production has not yet been fully understood.
Nowadays multiomics has been a common biological approach for systemic genomes [35][36][37][38]. In this paper, we prepared the first RNA-seq libraries from F. velutipes samples at three developmental stages, that is, mycelium liquid culture stage (FrI), young fruiting body stage (FrII) and mature fruiting body stage (FrIII). The transcriptome technique was used to identify the DEGs involved in the biosynthesis of ergosterol in F. velutipes at different developmental stages. Afterwards, the differential metabolites were completely scanned by non-targeted metabolomic techniques. Differentially expressed genes (DEGs) combined with differential metabolites were used to assist in identifying genes associated with ergosterol biosynthesis. The results showed that a large number of ergosterol biosynthesis-related genes were confirmed in the transcriptome of F. velutipes. The results revealed 19 DEGs related to ergosterol biosynthesis of F. velutipes at different developmental stages. After transcriptomic and metabolomic analysis of the ergosterol biosynthesis, four genes ( ERG1s, ERG10s, ERG11s, ERG26s) were identified as the most important genes in the regulation of the pathway flux towards ergosterol biosynthesis in F. velutipes. And the results showed that there might be 17 differential metabolites from the ergosterol biosynthesis in F. velutipes between different developmental stages. All the data and results had a vital significance for understanding the metabolic pathway of ergosterol biosynthesis in F. velutipes.

Materials and sample collection
This article used F. velutipes samples at three different developmental stages as experimental materials. They were mycelium liquid culture stage (FrI), young fruiting body stage (FrII) and mature fruiting body stage (FrIII) (Fig. 2). Fresh mycelium, young fruiting body and mature fruiting body were acquired from a F. velutipes farm in Xianyang City, Shaanxi Province, China in June 2017. Mycelium fluid culture and the culture conditions could be consulted in literature [3]. These F. velutipes samples at three developmental stages were collected under sterile conditions, immediately frozen and stored in liquid nitrogen at -80 ℃ before the metabolic content determination and RNA isolation [2].

RNA extraction, library preparation and mRNA-Seq
Total RNA was extracted with RNA pure Plant kit (Tiangen Biotech Co., Ltd., Beijing, China). The total RNA concentration, RIN value, 28S/18S and fragment size were determined using an Agilent 2100 Bioanalyzer (Agilent Technologies Co. Ltd., Santa Clara, CA, USA) with Agilent RNA 6000 Nano Kit. The purity of the samples was measured using an ultraviolet spectrophotometer NanoDrop TM . After the isolation and fragmentation of the total RNA, the eukaryotic mRNA was enriched by using Oligo (dT) coupled to magnetic beads. The preparation of cDNA libraries was based on the method of PCR described in the following procedure. The first strand cDNA was manufactured from the fully spliced mRNA as a template, and then a second strand reaction system was used to make doublestranded cDNA. The second strand cDNA was obtained and purified by a system kit. When the cDNA ends were repaired, an 'A' nucleobase hybridized to the 3᾽ ends of the cDNA.
Then the cDNA fragments of different sizes were selected and then amplified. Then the cDNA fragments of different sizes were selected and then amplified. The qualities of amplified cDNA libraries were checked and tested. At last, the sequencing was run on an Illumina HiSeq™ 4000 platform.

De novo assembly and functional annotation
First, we filtered out low-quality sequences, contaminated adaptors and unknown N base (N percentage > 5%). The resultant data is called clean reads. We used Trinity platform (http://trinityrnaseq.github.io/) for de novo transcriptome assembly to obtain the clean reads (removing PCR repeats to improve assembly efficiency). Then we assembled the clean reads using Tgicl to get unigenes, and performed functional annotation on unigenes. FPKM calculation and DESeq2 differential expression analysis We used alignment package Bowtie2 (v 2.2.5) to map clean reads back to unigenes. Then the gene expression levels of every sample were calculated by RSEM (v1.2.12) (RNA-Seq by expectation maximization) [39,40] which is a software package for RNA-seq reads to calculate the expression of its genes and transcript isoforms. This paper used the hclust function in the R software for hierarchical clustering analysis. Fragments Per Kilobase of transcript per Million (FPKM) values were calculated by mapping the reads to fragments, which could be used to quantify the abundance of the transcripts and so to analyze differential gene expression. Differential expression analyses between genes at different developmental stages were performed with DESeq2 package [41]. Differentially expressed genes (DEGs) were revealed by adjusted value (AP ≤ 0.05) and fold change analysis (FC ≥ 2.00).

GO and KEGG enrichment analysis of DEGs
According to the official classification categories and based on the GO and KEGG annotation results, the DEGs are annotated to functional and biological pathways, and the phyper of Rstatisticspackage isused for enrichment analysis. A p-value of 0.05 ( p ≤ 0.05) implies that the accepted false discovery rate (FDR) is 5%.

RT-qPCR validation and DEGs analysis
In order to validate internal control genes for expression analysis, ABI StepOnePlusTM RT-PCR System (Applied Biosystems, USA) was used to perform RT-qPCR and analyze the expression of DEGs related to ergosterol biosynthesis. Primer Premier 5.0 software was adopted to build sequencing primers, presented in Additional file 1: Table S1. The acquisition of cDNA and RT-qPCR were carried out by PrimeScript TM RT reagent Kit and SYBR Green TM Premix Ex Taq TM П (Takara Biotech Co., Ltd., Japan). GAPDH was used as an internal control. Three independent technical replicatesand three biological replicatesfor each sample were run to measure and assess the performance of RT-qPCR.

Measurement of metabolites
In this article, we also performed a metabolomic analysis. The tests were performed on six biological replicates for each metabolome. Tissue samples stored at -80 ℃ were placed in a freezer at -20 ℃ for 30 min and then thawed in a 4 ℃ freezer. 25 mg tissue was weighed, and placed in the EP tube; 800 ul of a cooled solution of methanol/water (1:1) was added to each EP tube and two small steel balls were frozen and put to each EP tube; the sample was put in the tissue lyser and the parameters was set to 35 Hz for 4min. After grinding, the ball was removed and the tube was placed in a -20 ℃ refrigerator for 2 h. the sample was centrifuged at 30,000 g for 15 min at 4 ℃. The EP tube was carefully removed from the centrifuge and 550 ul of each sample was transferred into a new EP tube; the sample was placed on the rack of the centrifuge and a picture was taken according to the order of samples on the task sheet; The sample will be handed over to be analyzed by LC-MS, and the remaining original sample will be handed over to the sample manager for storage.
All samples were made and processed by following protocols for LC-MS data analysis.
Firstly, all chromatographic isolation was performed by using UPLC system (Waters, UK).
The column type, column temperature, mobile phase flow rate, mobile phase solvent ratio and program settings are the same as the reference [42]. The injection volume was 10 μl.
The Q-TOF was run with the capillary tube and sampling voltages set at 3 kV and 40 V in the positive ion mode and at 1 kV and 40V in the negative ion mode. MS data was generated by Xevo G2 XS QTOF with TOF mass ranging from 50 to 1200 Da and scanned 0.2 s. The MS/MS was conducted to select, separate and detect precursor ions using 20 and 40 eV in different steps and with scan time 0.2 s. When processing the data, the LE signal was received at a 3 s interval to measure the mass accuracy. And a control sample was picked up from every 10 samples for evaluation of the performance of the LC-MS on the data acquisition.

Results
De novo assembly and analysis of RNA-seq data Assignments of KOG were used to predict and classify the possible functions of the unique sequences, and describe gene evolutionary processes. In this study, 12,427 unigenes (43.87% of all assembled unigenes) had KOG annotation functions and the KOG-annotated putative proteins were classified into 25 functional groups. 24.23% of the proteins were annotated to General functional prediction only the largest group, which is thought to be responsible for processing the processors of storage proteins into mature proteins, followed by signal transduction mechanisms (18.29%) and posttranslational modification, protein turnover, chaperones (13.24%). Just a small part of proteins was assigned to cell motility (0.22%). It is noteworthy that a great deal of proteins was annotated to secondary metabolites biosynthesis, transport and catabolism (8.57%), carbohydrate transport and metabolism (9.99%), transcription (7.44%), lipid transport and metabolism (7.23%) (Fig.   3).

GO enrichment analysis and DEGs annotation and assignment
This article also dealt with enrichment analysis on gene sets, and investigated how GO terms were represented using annotations for that gene set. The GO enrichment analysis tool Blast2GO package (v 2.5.0) was used to perform enrichment analysis. The results were interpreted and displayed. It was determined that a corrected FDR was the threshold.
De novo transcriptome assembly and analysis revealed 10,266 DEGs between FrII and FrI group, 10,467 DEGs between FrIII and FrI group and 2,677 DEGs between FrIII and FrII group (Additional file 4: Fig. S2). 9,265, 10,254 and 2,384 genes from the three DEG sets were assigned at least one GO term. 49 significant shared terms were displayed in Fig. 4.
The results showed that metabolic process, cellular process, single-organism process, localization, biological regulation, cellular component organization or biogenesis and regulation of biological process were significantly shared GO terms in biological process category. Cell, cell part, membrane, organelle, membrane part, macromolecular complex and organelle part were the most shared in cellular component category. Catalytic activity, binding, transporter activity, structural molecule activity and enzyme regulator activity were markedly shared in molecular function category.

KEGG function enrichment of DEGs
The genes were divided into 7 categories according to the KEGG Orthology classification for KEGG pathway maps. Metabolic pathway analysis provided insight into the biological functions and gene interactions and we found the most represented category was metabolism (Additional file 5: Fig. S3). From the bubble map of the DEGs pathway enrichment analysis (only the top 20 metabolic pathways shown) (Fig. 5), we found that the two metabolic pathways related to ergosterol biosynthesis with significant enrichment were terpenoid backbone biosynthesis pathway (ko00900, 32 DEGs, 59 orthologs annotated) and steroid biosynthesis pathway (ko00100, 43 DEGs, 78 orthologs annotated).

DEGs related to ergosterol
We analyzed the DEGs involved in terpenoid backbone biosynthesis pathway and steroid biosynthesis pathway in F. velutipes, which were related to ergosterol biosynthesis. Venn diagram and heatmap showed the differentially expressed genes of F. velutipes between three different developmental stages (Fig. 6). The results between FrII group and FrI group showed that 51 key unigenes (12 up-regulated and 39 down-regulated) were significantly differentially expressed. In addition, we analyzed transcriptional profiles of the differential genes between FrIII group and FrI group, to find a total of 56 key significantly differentially expressed unigenes (14 up-regulated and 42 down-regulated) shown in Table   2 and Fig. 6C. The results for the analyses showed strong similarity in differential expression between FrIII -FrI DEGs and FrII -FrI DEGs. Likewise, we found that the both metabolic pathways had only 8 DEGs (4 up-regulated and 4 down-regulated) between FrIII and FrII group (Fig 6A and 6B). These results indicated that during the formation of fruiting body, the difference of gene expression between the first and second stage was very distinct in ergosterol biosynthesis, but was quite small between the second and third stage. The distinctive gene expression difference may be due to the dramatic morphological changes in the first two stages of F. velutipes, while the gene expression similarity could be explained that change in size did not significantly affect gene expression at all. Among them, more than one unigenes were identified through annotation as the same enzyme. The results in Fig. 6C demonstrated that the four genes ( ERG10s, ERG1s, ERG11s, ERG26s) could be the most important genes for ergosterol biosynthesis in F. velutipes. The information of ergosterol biosynthesis in F. velutipes concerning the two pathways in the three different developmental stages was obtained.
Then, we could learn more about the ergosterol biosynthesis pathway in F. velutipes.
To validate the reliability of transcriptome sequencing data, the sequences of 15 core differentially expressed unigenes were analyzed by RT-qPCR primers. The results of RT-qPCR analysis showed a close similarity to the results of RNA-Seq, as shown in Fig. 7.
Analysis of biosynthesis intermediates and ergosterol from F. velutipes at the three different stages In this article, we performed metabolome sequencing of F. velutipes at three different developmental stages. PCA analysis of F. velutipes metabolome at different developmental stages was shown in Fig. 8A and 8B. Metabolic pathway analysis was based on the KEGG database. To understand the ergosterol biosynthesis, we compared the metabolic profiles of F. velutipes groups at different developmental stages (Fig 8C). Differential ions and identification results for different developmental groups were shown in Additional file 7: Table S3. LC-MS revealed 1,742 differential metabolites between FrII and FrI group, of which 732 were found up-regulated respectively and 1,010 down-regulated in the positive ion modes; and revealed 2,154 differential metabolites between FrII and FrI group, of which 925 were found up-regulated respectively and 1229 down-regulated in the negative ion modes. Similarly, between FrIII and FrII group, revealed were 751 and 944 differential metabolites respectively in positive and negative ion modes, of which 246 and 197 were up-regulated, and 505 and 747 were down-regulated.
To determine the metabolomic performance, we measured the end product ergosterol in the ergosterol biosynthesis. The results were shown in Additional file 10: Fig. S5. It was found that the m/z and retention time of the metabolomic results were consistent with the validation measurement, indicating that the metabolomic results were reliable.
Combined analysis of transcriptome and metabolome of ergosterol biosynthesis in F. velutipes at different developmental stages Systematic analyses of metabolites and genes were intended to investigate the relationship between genetic control of metabolite levels and metabolic impact on gene expression [43]. This article makes a comparison between the profiles of metabolites and gene expression in F. velutipes in different developmental stages using Pearson᾽s correlation coefficient (Additional file 11: Table S6 and Additional file 12: Table S7).
Metabolite datasets of FrII vs. FrI group and FrIII vs. FrII group were extracted and clustered using correlativity. And the cluster of gene expression profiles was made using the same method. Then, we combined the metabolite and gene expression profiles by calculating correlation coefficients. The heatmap in Fig. 9A and 9B illustrated the combined results. The horizontal axis represents metabolites and the vertical axis represents genes with the distance between each metabolite or gene reflecting similarity between their profiles. Red indicates positive correlation between metabolites and genes, and blue negative correlation. The metabolite-gene heatmap helps us understand the correlations between them.

Discussion
Dynamic changes of metabolites in F. velutipes at different developmental stages F. velutipes is one of the most medicinal and edible fungi, which is secure and has been well received by consumers. Since the successful cultivation of F. velutipes in 1980s, its outputs and consumptions have rapidly increased cross the world [2,44]. As the yield of F. velutipes increases, its nutritional and medicinal values become the focus of interest. In this study, by using metabonomic analysis, we could find metabolites involved in the ergosterol biosynthesis underwent significant change in their content during the three developmental stages (Additional file 8: Table S4 and Additional file 9: Table S5). In Fig.   8A and 8B, the results indicated that the biological replicates of F. velutipes in three different developmental stages were good and significantly different. In the course of three different growth stages (Fig. 8C), the expression level of the only metabolite i.e. lanosterol increased in the first two stages and then decreased in the latter two stages, and the expression level of the metabolites, including ergosta-5,7,22,24-tetraenol, squalene-2-3-epoxide and 14-demethyl lanosterol, decreased in the first two stages and then decreased in the latter two stages. We found that some metabolites such as mevalonate-5-phosphate, 4-methylzymosterol, episterol, ergosta-5,7,24(28)-trienol and ergosterol were gradually increasing during the three growth stages , and the only metabolite mevalonate was gradually decreasing during the three growth stages. These indicated that the content of sterols in the fruiting bodies tended to increase in mature fruiting body stage. These results might be closely related to the expression of genes and the substrate content in their metabolic processes in F. velutipes [45]. Based on the above results, we found that the nutritional and medicinal values of F. velutipes were higher in the mature fruiting stage. Our results provide scientific grounds for health and food guidance, but more definite conclusion takes more effort and research.

Ergosterol biosynthesis pathway in F. velutipes
At present, effective approaches to efficient ergosterol production from the mycelium or fruiting body of fungus cannot be devised until the genetic engineering based on the metabolic pathway is well understood. Although the biosynthesis pathway of ergosterol in S. cerevisiae has been well characterized, few efforts have explored the ergosterol biosynthesis in F. velutipes [33]. In this paper, the extracts of F. velutipes samples in different developmental stages were analyzed by LC-MS, and the sterols were studied thoroughly. At the same time, the expression of key enzymes in the ergosterol biosynthesis in F. velutipes indifferent developmental stages was established by RNA-Seq.
The results showed that the four genes ( ERG10s, ERG1s, ERG11s and ERG26s) might be the most important genes involved in the biosynthesis of ergosterol. In previous studies, ERG1 and ERG11 are identified as the key regulators for the post-squalene biosynthesis and in the feedback regulation of ergosterol biosynthesis in S. cerevisiae and Trichoderma harzianum [46,47], which our results were consistent with and whose conclusion our results supported. In S. cerevisiae, the deletion of the ERG26 was lethal and disrupt the synthesis of ergosterol [48,49]. These results indicated that ERG26 was essential for cell growth, impacting the synthesis of ergosterol. The ERG10 gene encodes an acetoacetyl-CoA thiolase that catalyzed the formation of acetoacetyl-CoA by two acetyl-CoA molecules.
Studies have found that acetoacetyl-CoA reduced the activity of Erg10p when sterols were more than required [50], which our study supported. others᾽ work confirmed that this regulation occurred during the transcriptional process. When the levels of some sterols in the cell are low, ERG10 gene is expressed to a higher level and then regulates the activation pathway [51]. The efficiency of ergosterol biosynthesis was determined by ratelimiting enzymes and more crucially by optimal coordination of all enzymes [52].
Transcription factor UPC2 was reported to upregulate target genes involved in the biosynthesis of sterol through activating the sterol response elements in their promoter regions [53,54]. The core motifs of sterol-response elements had been identified in nine responsive Candida albicans ERGs ( ERG1, ERG2, ERG5, ERG6, ERG10, ERG11, ERG24, ERG26 and ERG27) . Our results revealed that the expression levels of most of these genes had differed, which was consistent with the regulatory effect of UP2C on these genes. And most of these genes were found to be related to the post-squalene pathway, which promised to improve the sterol biosynthesis in F. velutipes [26,[54][55][56]. However, overexpression of specific enzymes could result in imbalance of sterol intermediate accumulation, and significantly reduced total cell biomass yield and repressed endproducts. Therefore, it was necessary to avoid the over-expression of some specific enzymes and control the accumulation of some cytotoxic intermediates [57]. In summary, a balance between precursor supplies and catalytic activities of enzymes should be struck to achieve optimal production of ergosterol.
In Fig. 9, combining analyses of differential metabolites and genes were expected to find regulation relationships between them. We made a comparison between profiles of metabolites and gene expression between three different stages using Pearson᾽s correlation coefficient (Additional file 11: Table S6 and Additional file 12: Table S7). The horizontal axis and the vertical axis respectively represent metabolites and genes. With the distance between each metabolite or gene reflecting the similarity between their profiles [43]. Red indicates positive correlation between metabolites and genes, and blue negative correlation. This could be a useful method for comparing correlations of metabolites or genes between different groups. The heatmap of metabolites and genes helps us to understand the correlations between them. Thus, the systematic analysis in this study can be used for further research of metabolic and gene functions for bioproduction. The heatmap in Fig. 9A

Availability of data and materials
The datasets supporting the conclusions of this article are included with in the article and its additional files.

Competing interests
The authors declare that they have no competing interests.

Funding
This study was financially supported by Science and Technology Coordination Innovation  Table 1 Transcript quality standards of the RNA-Seq analysis of F. velutipes Table 2 The unigenes related to ergosterol biosynthesis of F. velutipes in different developmental stages   Each unigene has three biological replicates and three technical replicates