Comparative analysis of fatty acid metabolism based on transcriptome sequencing of wild and cultivated Ophiocordyceps sinensis

Specimen collection

Naturally growing O. sinensis were harvested at Aba Prefecture, Sichuan Province, China (4,500 m, N31°08′51.9″, E102°21′26.78″) in May. Cultivated samples were collected from the artificial cultivation workshop at Chengdu Eastern Sunshine Co. Ltd. The genetic background of Hepialidae larvae and H. sinensis used in the cultivation is consistent with the wild genetic background. Cultured conditions: one month after the artificial infestation has successfully formed sclerotium, keeping light for 10 h a day, air humidity above 80%, soil humidity above 20%, temperature18 °C, and adding nutrients to the soil according to the proportion of the large elements of MS medium. All samples were stored at −80 °C until further processing.

Morphological comparison

To identify morphological differences between the samples, we first compared the appearance of wild and cultivated O. sinensis and recorded photographs. Subsequently, crosscutting was performed at the top and tail of the worm body and the fruiting body at one cm, including the middle part, and photographed using Toup View imaging system. Finally, microscopic analysis was performed to further observe the microstructure characteristics, (Jie-Pin, 2006). Paraffin sections were prepared and observed under an electron microscope. Three replicates were conducted for each set of samples.

Transcriptome data acquisition

Total RNA was extracted using Eastep Super Total RNA Extraction Kit (Promega Shanghai) according to the manufacturer’s instructions. After assessing the RNA quality, 1 µg of total RNA from each group was used as input material for RNA sample preparation using the NEBNext^R Ultra™ Directional RNA Library Prep Kit for Illumina^R (NEB, USA) to construct cDNA libraries. Pair-end sequencing was performed using the Illumina Hiseq2000 platform. A total of 100-base paired-end clean reads were generated after removing low-quality reads, adaptor-polluted reads, and reads with high levels of an unknown base. Two biological replicates were performed for RNA-seq.

Gene expression analysis

Raw data were preprocessed using the fastp tool (http://github.com/openGene/fastp) to obtain clean reads. The Tophat2 aligner tool (http://ccb.jhu.edu/software/tophat/index.shtml) was then used to align sequencing reads to the reference sequences (Daehwan et al., 2013). To eliminate the effect of different sequencing discrepancies and gene lengths on the gene expression calculation, fragments per kilobases of exon per million fragments mapped (FPKM) values were used to normalize the expression level (Robinson & Oshlack, 2010). To infer the transcriptional changes in both groups, differential expression analysis was performed using the DESeq2 R package according to the express count matrix (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) (Love, Huber & Anders, 2014). A false discovery rate (FDR) ≤0.01 and —log2(fold change, FC)—≥2 were set as thresholds for the selection of differentially expressed genes (DEGs).

Functional annotation and enrichment of DEGs

All DEGs were annotated using five databases, including the NCBI nonredundant protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG), Universal Protein, and the Cluster of Orthologous Groups of proteins (COG). Gene Ontology (GO) annotations were classified into three ontology categories: molecular function, cellular component, and biological process (Ashburner et al., 2000). KEGG enrichment analysis of the DEGs was performed using the KOBAS software (http://kobas.cbi.pku.edu.cn/kobas3/) (Minoru & Susumu, 2000) and showed the top 20 pathways with the smallest significant Q value.

qRT-PCR

qRT-PCR was performed on the CFX ConnectTM Optics Module using 2X Ultra SYBR Mixture (TransGen, Beijing, China) according to the manufacturer’s instructions. The β-actin gene was used as a reference to normalize the expression data, and the 2^−ΔΔCT values were shown as relative expression levels (Thomas, Kenneth & Livak, 2008). All primers used for the assessment of transcript levels are listed in Table S1.

Fatty acid extraction and identification using GC-MS

We collected three batches of wild and cultivated O. sinensis specimens for fatty acid testing (Table S2). The samples were prepared dried at 60 °C and grounded to obtain 80∼100 mg of powder. two mL 5% hydrochloric acid solution-methanol (1 :1, V/V) and three mL chloroform-methanol were added, then water bathed at 85 °C for 1 h. Next, when cooled to room temperature, one mL of n-hexane was mixed and extracted for 2 min. 100 µL of the supernatant was taken into a one mL volumetric flask and n-hexane was used to make the volume constant. Finally, the fatty acid preparation sample was obtained by filtration through a 0.45 µm microporous membrane (Zhihui et al., 2017). All the extracts were analyzed using an Agilent 7890B/5977B gas chromatograph. A TG-5MS capillary column (30 m × 0.25 mm × 0.25 µm) was used with helium (99.999%) as the carrier gas at a constant flow rate of 1.2 mL min⁻¹. The sample (0.1 µL) was injected in the split-less mode, and the inlet temperature was 290 °C. The oven temperature program was as follows: the initial temperature was maintained at 80 °C for 1 min, increased to 200 °C at the rate of 10 °C/min, increased from 200 °C to 250 °C at the rate of 10 °C/min, and maintained at 250 °C for 3 min. Finally, the temperature was increased to 270 °C at the rate of 2 °C/min. The qualitative and quantitative analysis of fatty acids was performed according to the GC-MS NIST Mass Spectral Library and the chromatogram of mixed acid fatty standard (NU-CHEK, USA, M20-D, 25mg).

Data availability

The O. sinensis reference genome was derived from the Ensembl Fungi database, ASM44836v1 (http://fungi.ensembl.org/Ophiocordyceps_sinensis_co18_gca_000448365/Info/Index). The raw data of wild O. sinensis were deposited in the Genome Sequence Archive under the accession number PRJCA000970, and the data of cultivated O. sinensis are available at the NCBI Sequence Read Archive under accession numbers SRR5282577 and SRR5282578.

Morphological comparison between wild and cultivated O. sinensis

To perform a preliminary comparison between wild and cultivated O. sinensis, the physical characteristics, cross-sectional characteristics of each part, and microscopic observations of paraffin sections were evaluated to determine the morphological differences associated with the developmental phenotype of O. sinensis. The structural characteristics of wild and cultivated O. sinensis are consistent with the description in the Chinese Pharmacopoeia: “The insect body is connected to the fungus seat, the surface is dark yellow to yellowish brown, and eight pairs of feet” (Fig. S1), and the main difference is reflected in the morphology. The cross-section and microscopic results of each part are shown in Fig. 1, Figs. S2 and S3. Compared with cultivated O. sinensis, the outer wall structure is thinner, and the digestive tract and other organs in the body are less complete. In addition, the inner hyphae are arranged in parallel in the fruiting body, the hyphae in the wild specimens appear to be arranged more closely, and the density is significantly higher than that in cultivated O. sinensis specimens. Such differences are mainly caused by their physiological functions and growth environment.

Overview of transcriptome and differential analysis between wild and cultivated O. sinensis

After performing quality control of the sequencing data, a total of 28.85 Gb of clean data was obtained, and the percentage of Q30 bases was ≥93.68% (Table S3), indicating that the higher the accuracy of base identification during the sequencing process, the sequence quality can be used for subsequent analysis. The efficiency of comparison with the reference genome was between 81.00% and 84.40% (Table S4), and 9360 transcripts were obtained after mapping, including 997 newly predicted genes.

DEGs were calculated based on FPKM, with an FDR of <0.01 and —log2(fold change, FC)— ≥2. A total of 563 DEGs in the wild and cultivated O. sinensis were identified. Compared with the wild group, the expression of 415 genes was downregulated and the expression of only 148 genes was upregulated in the cultivated group. The number of DEGs in the wild group was significantly higher (Fig. 2).

Functional enrichment analysis of DEGs

Overall, 563 DEGs were annotated by Blast2GO to a total of 274 annotation entries. The GO term cellular component was mainly enriched into 14 terms, including membrane, organelle, macromolecular complex, etc. In the GO term molecular function, a total of 225 DEGs were enriched into 8 terms, mainly catalytic activity, binding, transcription regulator activity and transporter activity, etc. In biological process, 193 DEGs were significantly enriched in 14 terms, mainly metabolic process, cellular process, and single-organism process (Fig. 3A). KEGG annotation and enrichment results are shown in Table S5 and Fig. 3B. Fatty acid degradation, fatty acid metabolism, butanoate metabolism, and β-alanine metabolism showed the most significant differences in the metabolic pathways, of which the most abundant was biosynthesis of antibiotics.

Based on these results, we found that fatty acid metabolism was the most significant difference between wild and cultivated O. sinensis. On comparing these data with major functional databases, the DEGs in the fatty acid pathway were identified and were found to code for acyl-CoA dehydrogenase (ACAD) (OCS_02728), enoyl-CoA hydratase (ECH) (OCS_06591), 3-ketoacyl-CoA thiolase (KAT) (OCS_04932), and acetyl-CoA acetyltransferase (ACAT) (OCS_01242/OCS_05099).

Fatty acid metabolism-related DEGs involved in wild and cultivated O. sinensis

Based on the results of functional enrichment analysis of DEGs, we selected fatty acids as the focus of research on the differences between wild and cultivated O. sinensis. Acetyl-CoA (Ac-CoA) is the precursor of the Malongl-CoA (MalCoA) in fungi and is the final product of fatty acid degradation. Under the catalysis of fatty acid synthase (FANS) and other enzymes, one molecule of MalACP is added in each step, increasing the length of two carbon atoms each time; after seven additions, hexadecanoyl-[acp] is formed and, finally, palmitic acid is formed under the action of FASN (Fig. 4A) (Nowinski et al., 2018). In addition, hexadecanoyl-CoA can be used as a substrate for the synthesis of longer-chain fatty acids, generating different types of fatty acids under the catalysis of different enzymes.

In this study, the enzymes encoded by the DEGs of the significantly enriched metabolic pathways participated in the O. sinensis fatty acid β-oxidation pathway, including four steps: dehydrogenation, hydration, re-dehydrogenation, and thiolysis (Shen & Burger, 2009). ACAD, ECH, and KAT were significantly downregulated in cultivated O. sinensis. ACAD participates in the first reaction of fatty acid and amino acid β-oxidation under the action of flavin adenine dinucleotide cofactor (Kim & Miura, 2004). ECH catalyzes the second step of fatty acid β-oxidation, promoting the cis-addition of water molecules on the trans-2-enoyl-CoA thioester double bond to form β-hydroxyacyl-CoA thioester (Willadsen & Eggerer, 2008). KAT catalyzes the final step of fatty acid β-oxidation to produce one molecule of Ac-CoA and one molecule of acyl-CoA (Mano et al., 1996). The β-oxidation that occurs in mitochondria and peroxisomes is the most important form of oxidation and provides a large amount of energy for organisms by coupling with the tricarboxylic acid cycle and electronic respiratory chain (Neely & Morgan, 1974). The results showed that the genes encoding various enzymes involved in the degradation and synthesis of fatty acids were downregulated in cultivated O. sinensis, which indicates that fatty acid metabolism is more active in wild O. sinensis (Fig. 4B). Moreover, the product of palmitic acid decomposition is Ac-CoA. When fatty acid metabolism is weakened, the downstream citric acid cycle/tricarboxylic acid cycle, alanine aspartate metabolism, and butyrate metabolism also weakens (Fig. 4A). Because of the long growth cycle of wild O. sinensis and the lack of nutrients in the living environment, we speculate that O. sinensis relies mainly on its own fatty acid oxidative decomposition to provide energy for itself, and on β-oxidation, especially via palmitic acid decomposition.

Validation of transcriptome data by qRT-RCR

To confirm the reliability of the sequencing data, 12 DEGs were randomly selected to validate the RNA-Seq expression profiles. As expected, qRT-PCR results showed that most of the mRNAs shared similar expression with those from the sequencing data. The expression levels of ten genes were downregulated in the cultivated O. sinensis (Fig. 5) and those of two genes were upregulated, indicating that our data were accurate and reliable in the subsequent analyses.

Quantitative analysis of fatty acid content

GC-MS was performed to analyze the fatty acids of O. sinensis. The results showed that the types of fatty acids in the wild and cultivated specimens were basically the same (Fig. 6) but that the content was quite different. The content of fatty acids in the wild group was much higher than that in the cultivated group, which is consistent with the previously reported (Shen & Burger, 2009). Furthermore, among all types of fatty acids, oleic acid (C18:1n9c), linoleic acid (C18:2n6c), and palmitic acid (C16:0) are the most important components and account for up to 89.38%–94.88%. From the perspective of fatty acid types, the content of unsaturated fatty acids in O. sinensis was higher than the content of saturated fatty acids. In addition, the content of fatty acids in O. sinensis has a certain aging effect. With the extension of storage life, the content of fatty acids in medicinal materials is continuously reduced (Table S6).