High Throughput Identification of the Potential Antioxidant Peptides in Ophiocordyceps sinensis

Ophiocordyceps sinensis, an ascomycete caterpillar fungus, has been used as a Traditional Chinese Medicine owing to its bioactive properties. However, until now the bio-active peptides have not been identified in this fungus. Here, the raw RNA sequences of three crucial growth stages of the artificially cultivated O. sinensis and the wild-grown mature fruit-body were aligned to the genome of O. sinensis. Both homology-based prediction and de novo-based prediction methods were used to identify 8541 putative antioxidant peptides (pAOPs). The expression profiles of the cultivated mature fruiting body were similar to those found in the wild specimens. The differential expression of 1008 pAOPs matched genes had the highest difference between ST and MF, suggesting that the pAOPs were primarily induced and play important roles in the process of the fruit-body maturation. Gene ontology analysis showed that most of pAOPs matched genes were enriched in terms of ‘cell redox homeostasis’, ‘response to oxidative stresses’, ‘catalase activity’, and ‘ integral component of cell membrane’. A total of 1655 pAOPs was identified in our protein-seqs, and some crucial pAOPs were selected, including catalase, peroxiredoxin, and SOD [Cu–Zn]. Our findings offer the first identification of the active peptide ingredients in O. sinensis, facilitating the discovery of anti-infectious bio-activity and the understanding of the roles of AOPs in fungal pathogenicity and the high-altitude adaptation in this medicinal fungus.


Introduction
Ophiocordyceps sinensis (Berk.), syn. Cordyceps sinensis, belongs to the family Ascomycetes and is a highly valued Traditional Chinese Medicine (TCM) fungus used in Asian countries for over 2000 years [1,2]. Its common name is the Chinese caterpillar fungus, and it naturally inhabits the alpine Qinghai-Tibetan Plateau in South East Asia with an altitude of 3000-5000 m above sea level [3]. Over 20 bioactive ingredients have been reported in this species, including adenosine, cordycepic acid, ergosterol, and polysaccharides, which are thought to be responsible for a range of health benefits, including anti-inflammatory, anti-tumor, immunomodulating, and antioxidative activities [4]. Bioactive peptides have been discovered from a range of prokaryotes and eukaryotes [5]. While bioactive peptides from fungi source account for a small proportion, mainly from Trichoderma viride [5]. To date, bioactive peptide compounds have not been identified in O. sinensis.
O. sinensis is an entomopathogenic fungus and has a complex life cycle. This parasitic fungus forms a unique complex with the larva of ghost moth caterpillars (Thitarodes spp.), and the larva progressively becomes stiff and coated with mycelia to mummify the larvae [6]. A small stroma bud emerges from the head of the sclerotium and forms the stalked fruiting body [6]. The immune response against the fungi would produce plenty of reactive oxygen species (ROS) to prevent pathogens [7]; furthermore, the fungi and the host itself would promote its ROS antioxidant defense system [7]. Omics studies showed Reads were aligned to the genome from each sample and then assembled into transcripts by StringTie (version 1.3.4d) using default parameters [19,20]. To globally characterize the expression patterns of diverse RNA-Seq samples, paired-end reads were aligned back to the assembled transcripts using Bowtie 2.0 as the aligner [19]. Gene expression patterns were quantified using htseq-count (version: 0.11.3) based on the read numbers that were mapped to each gene. The mapped read numbers of each assembled transcript were estimated and were normalized by RESM-based algorithm to obtain fragments per kilobase of transcript per million mapped fragments (FPKM) values for each RNA-Seq sample using perl scripts in the Trinity package (v2.11.0) [20]. RSEM results of each replicate of the sample were merged as one matrix for downstream analyses.
Differentially expressed genes (DEGs) were identified by using the edgeR package (empirical analysis of digital gene expression in R), with a threshold of adjusted log2FC(log2fold-change) ≥ 1 and false discovered rate (FDR) < 0.05 as statistically significant. A Venn diagram was used to analyze the pDAPs distribution in different comparisons. The RNA expression patterns of the common pDAPs with p-value < 0.05 in the three comparisons (IL vs. ST, ST vs. MF, and MF vs. YF were further clustered by one-way hierarchical clustering implemented in hCluster R package (Euclidean distance, average linkage clustering). The Gene Ontology (GO) based annotation suite BLAST2GO v. 3.0 [21] was used to functionally annotate genes and assess the quality of our assembly. The transcripts were aligned against three public databases (NR, Swiss-Prot, and KEGG).

Protein Preparation
The samples were derived from the three growth stages (ST, PR, MF) of the cultivated O. sinensis and the wild specimens (YF). The protein preparation method was used according to a previous study [26]. Briefly, 0.1 g of each fresh sample was ground into a powder in liquid nitrogen, suspended with a solution with 100:1 (w/v) urea lysis buffer containing 8 M urea, 2 mM EDTA, 10 mM iodoacetamide (DTT), 25 mM iodoacetamide (IAA), and 1% protease inhibitor cocktail (Roche), and it was sonicated for 1 min. The protein in the supernatant was precipitated by adding 3 volumes of acetone at −20 • C for 2 h. After centrifugation at 12,000 rpm for 15 min at 4 • C, the protein pellet was collected and air dried. Protein was extracted by resuspending the dry pellet in UT buffer (8 M urea, 100 mM triethyl ammonium bicarbonate (TEAB)). Protein concentrations were determined using the Bradford protein assay (Bio-Rad, Hercules, CA, USA) according to the manufacturer's instructions, and then samples were kept in −80 • C prior to use.

Peptide Isobaric Labeling
For each sample, 100 mg protein was reduced with 10 mM DTT for 1 h at 37 • C and then alkylated with 55 mM iodoacetamide (IAM) for 30 min in darkness. The pool protein was diluted by adding 100 mM TEAB to a final urea concentrations < 2 M and then digested with sequencing grade modified trypsin (Promega, Madison, WI, USA) at a ratio of 1:50 trypsin/protein (w/w) t for the first digestion at 37 • C overnight, and a ratio of 1:100 trypsin/protein (w/w) for a second 4 h digestion [27]. Subsequently, protein was desalted using a Strata X SPE column (Phenomenex Inc; Torrance, CA, USA) and vacuum-dried using a SpeedVac concentrator (Thermo, San Jose, CA, USA). Peptide was reconstituted in 20 mL 500 mM TEAB, processed according to the manufacturer's protocol for an 8-plex iTRAQ kit (Sigma, Aldrich, St. Louis, MO, USA), and incubated for 2 h at room temperature. The labeled peptide was desalted and dried by vacuum centrifugation.

HPLC Fractionation and High-Resolution LC-MS/MS Analysis Based on Q Exactive
Briefly, the iTRAQ-labeled samples were reconstituted in 0.1% formic acid (FA), injected onto an Acclaim PepMap R 100 C18 reversed-phase pre-column (3 µm, 100 Å, 75 µm × 2 cm, Thermo Fisher Scientific, San Jose, CA, USA) at 5 µL/min in 100% solvent A (0.1% FA in water). Separation of peptides was performed using a reversed-phase column (Acclaim PepMap R RSLC C18, 2 µm, 100 Å, 50 µm × 15 cm) with an increase gradient of 0-8% solvent B (0.1% FA in 98% ACN) over 5 min, 8-25% B over 35 min, and 25-98% B for 10 min, and then kept in 98% for 8 min. The flow rate was kept constant at 300 nL/min on an ultimate 3000 system. The eluent was sprayed via an NSI source at an electrospray voltage of 2.5kV and then analyzed by MS/MS in Q Exactive HF(Thermo Fisher Scientific, San Jose, CA, USA) [28]. The mass spectrometer was operated in data-dependent mode by tandem mass spectrometry (MS/MS). Full-scan MS spectra (from m/z 350 to 1800) were acquired in the Orbitrap with a resolution of 70,000. MS data were obtained by selecting up to the 15-most-abundant precursors ions present in the survey scan (300-1800 m/z) for decision-tree-based ion trap higher-energy collisional dissociation (HCD) fragmentation [27]. The HCD collision energy was set at 32% and a dynamic exclusion duration of 10.0 s.

Data Processing
The MS/MS raw data were analyzed using the Sequest software integration in Proteome Discoverer (version 1.3, Thermo Scientific) and searched against O. sinensis (CGMCC 3.14243). Search parameters were as follows: trypsin as the cleavage enzyme, maximum of 2 missed cleavages, carbamidomethyl used as a fixed modification and oxidation (M), N-Term Acetylation and iTRAQ labeling used as fixed modifications [27]. The databank searches were performed using a peptide mass tolerance of 20 ppm and a product ion mass tolerance of 0.05 Da and an FDR cutoff 0.05. The pAOPs were searched against the peptides from our proteome data.

Statistical Analysis
Statistical analysis of relative gene expression levels were calculated via the 2 −∆∆Ct method, using R (V3.2) to estimate the correlations between the mRNA expression levels of 9 the AOPs determined through qRT-PCR [29]. The histograms were plotted using GraphPad Prism software 8.0 (GraphPad, La Jolla, CA, USA).
Proteome and transcriptome analyses were carried out on at least three independent biological repeats for each sample. Other experiments were performed using at least three independent biological repeats and for each biological repeat, at least three technical repeats were performed.

Overview of Transcriptome and Differential Expression Analysis
Clean reads were obtained for each replicate of IL, ST, and MF, and aligned to the O. sinensis genome using HISAT2 [18], and assembled by StringTie, resulting in 65,534 exons and transcripts. For each replicate, over 20,694,375 reads were successfully mapped. The assembled results are shown in Table S1. The efficiency of comparison with the reference genome was between 91.69% and 97.68% (Table 1). The correlation analysis of differential expression patterns revealed that the three biological samples in each group showed similar high performance, and the samples of the three stages could be clearly assigned to four clusters ( Figure 1A). The IL and ST stages were found to group together, while the stages of YF and MF grouped together, indicating that the asexual hyphae stage (IF vs. ST) and the sexual stages (YF vs. MF) have more similar expression patterns than that in the comparison to ST and MF. DEGs were calculated based on fragments per kilobase per million (FPKM), with an FDR of <0.05 and |log2(fold change, FC)| ≥ 1. This threshold resulted in a total of 1176 genes as significant DEGs in IF vs. ST (657 genes up-regulated, 519 genes downregulated), 263 in MF vs. YF (503 genes up-regulated, 260 genes down-regulated), and 1916 in ST vs. MF (932 genes up-regulated, 984 genes down-regulated) ( Figure 1B, Table S2). Differential expression in our dataset is represented as a volcano plot ( Figure 2). These results revealed that the expression profiles of the asexual and sexual growth stages (YF and MF) and (ST and IL) were more similar to each other than when compared across growth stages (ST and MF). Previous studies showed that the DEGs in the adjacent growth stages were enriched primarily with the response to oxidative stress [6]. In this study, we mainly focused on the AOPs that might be generated during the infection and sexual development process.

AOP Identification and Functional Annotation
A total of 607,300 ORFs was predicted by ORF finder (Table S3). A total of 154,507 smORFs was identified by searching against databases. The length distribution of the smORF is shown in Figure 3A and Table S4, showing that most of ORFs were distributed within 40-60 bp. The 607,300 ORFs were used to identify putative AOPs (pAOPs) using

AOP Identification and Functional Annotation
A total of 607,300 ORFs was predicted by ORF finder (Table S3). A total of 154,507 smORFs was identified by searching against databases. The length distribution of the smORF is shown in Figure 3A and Table S4, showing that most of ORFs were distributed within 40-60 bp. The 607,300 ORFs were used to identify putative AOPs (pAOPs) using

AOP Identification and Functional Annotation
A total of 607,300 ORFs was predicted by ORF finder (Table S3). A total of 154,507 smORFs was identified by searching against databases. The length distribution of the smORF is shown in Figure 3A and Table S4, showing that most of ORFs were distributed within 40-60 bp. The 607,300 ORFs were used to identify putative AOPs (pAOPs) using the homology-based prediction method against the UniProtKB/Swiss-Prot databases. The results showed that totals of 2276 pAOPs were identified by the homology-based predic-tion method and 6269 pAOPs were identified by the denovo-based prediction method. Four AOPs were shared by the two methods listed in Table 2. A total of 8541 AOPs was identified and listed in Table S5. About 50% of the pAOPs contained more than 100 aa with an average length of 138 aa (Table S6). The longest antioxidant peptide had 1870 aa and the shortest one had merely 9 aa (Table S6). The length distribution of the pAOPs was shown in Figure 3B.

The Differential Expression Analysis of AOPs
The identified AOPs were searched against the DEGs. A total of 1008 pAOPs was identified between stages ( Figure 5, listed in Table S8 The number of DAPs in ST vs. MF was the highest, while that in MF vs. YF was the smallest, indicating that the differences in DAP expression primarily occurred between the stages of ST and MF. The Venn map showed that there were 33 common DAPs in the three comparisons ( Figure 5). Furthermore, RNAs profiles of the 33 DAPs were analyzed ( Figure 6). The results showed that there were five clusters with visible difference expression patterns ( Figure 6A-E) in the three comparisons. In cluster one and five, with six and two DAPs, respectively, there were lower expression levels in IL compared to in ST, and higher expression levels in ST compared to MF, indicating that they had a gradual increase during the serial growth stages from IL to MF and might be associated with the whole growth process. In cluster 2 and 3, with 15 and 5 pDAPs, respectively, there were higher expression levels in IL compared to ST and lower expression levels in ST compared to MF, a decrease expression profile and then an increase expression profile upon the transition from the asexual-to sexual-stage, illustrating that these pDAPs might play important roles in fungus-larva interaction and fruit-body maturation. In clusters three and five, with five and two pDAPs, respectively, there were lower expression levels in MF compared to YF, indicating that these DAPs might be associated with habitat-adaptation in this fungus. Combined with the results of the Heatmap (Figure 6F)

Validation by RT-PCR
To confirm the reliability of transcriptome analysis using Illumina sequencing, nine transcripts encoding pAOPs were selected for qPCR validation. Specific primers were designed as listed in Table S9. Except for two genes (MSTRG.5870, MSTRG.10795), the other results of qPCR conformed to that of RNA-seqs (Figure 7), supporting the accuracy of the RNA-Seq and differential expression AOPs (DAPs) analysis.

Validation by RT-PCR
To confirm the reliability of transcriptome analysis using Illumina sequencing, nine transcripts encoding pAOPs were selected for qPCR validation. Specific primers were designed as listed in Table S9. Except for two genes (MSTRG.5870, MSTRG.10795), the other results of qPCR conformed to that of RNA-seqs (Figure 7), supporting the accuracy of the RNA-Seq and differential expression AOPs (DAPs) analysis.  Figure 7. qRT−PCR validation of the pAOPs. The first four sets of bars represent the relative expression levels in each candidate AOPs identified in ST vs. IL (relative to IL), the middle three sets of bars represent that in ST vs. MF (relative to MF, and the latter three sets of bars represent that in MF vs. YF (relative to YF). Yellow bars represent qRT-PCR results (2 −∆∆Ct ). Error bars indicate the standard error. Purple bars represent the protein-seq results. 18sRNA was the internal reference. The X-axis represents the relative mRNA expression levels, and the Y-axis represents the AOPs gene ID. IL represents the mycoparasite complex of the cultivated O. sinensis; ST represents the mummified larvae coated with mycelia of the cultivated O. sinensis; MF represents the mature fruiting body of the cultivated O. sinensis; YF represents the wild-grown O. sinensis, and YF represents the wild-grown O. sinensis with mature fruit-body.

iTRAQ Analysis and pAOPs Validation
The proteome of three growth stages of the cultivated O. sinensis were investigated. For visualization, PCA was performed. A total of 1414 proteins was identified using iTRAQ in these samples of different stages, and 1129 proteins were quantified. We detected 2343, 2511, 2972, and 2972 unique peptides; 3747, 3384, 2808, and 4080 unique peptides; 1725, 1830, 2210, and 3266 unique peptides; and 2921, 2960, 2892, and 3400 unique peptides in ST, PR, MF, and YF, respectively. The peptides identified are listed in Table S10. In this study, the most spectra results had a mass error within ± 4 ppm, indicating that the mass spectrometer's mass accuracy was normal ( Figure 8A). The length of the peptides identified were distributed between 7 aa and 41 aa (Figure 8B), and 90% of the peptide length was distributed within 28 aa. MS data were deposited in iProX with the primary accession code PXD030687.
detected 2343, 2511, 2972, and 2972 unique peptides; 3747, 3384, 2808, and 4080 unique peptides; 1725, 1830, 2210, and 3266 unique peptides; and 2921, 2960, 2892, and 3400 unique peptides in ST, PR, MF, and YF, respectively. The peptides identified are listed in Table S10. In this study, the most spectra results had a mass error within ± 4 ppm, indicating that the mass spectrometer's mass accuracy was normal ( Figure 8A). The length of the peptides identified were distributed between 7 aa and 41 aa (Figure 8B), and 90% of the peptide length was distributed within 28 aa. MS data were deposited in iProX with the primary accession code PXD030687.
A total of 1655 pAOPs was found in our protein-seqs, shown in Table S11. Combined with previous studies, some crucial pAOPs are listed in Table 3,    Peroxide stress-activated histidine kinase mak1 OS A total of 1655 pAOPs was found in our protein-seqs, shown in Table S11. Combined with previous studies, some crucial pAOPs are listed in Table 3

Discussion
O. sinensis has been reported to have various biological activities that are of relevance for development in pharmaceutical products. Recent studies have mainly focused on its sexual development and pharmacological studies [8]. Omics studies suggested that antioxidants would be largely generated during interactions between fungal pathogens and host insects and the fruit-body formation, as well as the high-altitude adaptation [16,30,31]. Here, using transcriptome and proteome data of three serial growth stages in cultivated O. sinensis and the wild-grown mature fruit-body, the putative AOPs were firstly predicted and analyzed. The correlation analysis of differential expression genes revealed that the three biological samples in each group showed a higher similarity, indicating that samples within each group had good repeatability. The expression patterns in IL vs. ST and YF vs. MF had a higher similarity compared to that in ST vs. MF, indicating that a higher difference occurred in the stages of asexual growth and sexual development. Here a combination of homology-and denovo-based prediction methods were used to identify AOPs. Homologous-based approaches used BLAST searching against sequences deposited in the present antioxidant peptides databases. As the present databases with antioxidant peptides are sparse, a denovo-based prediction method was also used to identify AOPs. A total of 8541 pAOPs was identified. Only four pAOPs were shared by the two analysis methods, indicating that most of the pAOPs have not yet been identified and needed to be verified, and the four pAOPs needed to be further studied.
GO analysis showed that 27 pAOPs enriched in biological processes were involved in 'cell redox homeostasis' and 'response to oxidative stresses'. A total of seven pAOPs was enriched in the mature fruiting body for 'catalase activity'. It is well known that AOPs are induced in response to oxidative stresses, the maintenance of cell homeostasis, and the involvement of large molecular antioxidants such as superoxide dismutase (SOD), catalase (CAT), glutathione peroxidase (GPx), and glutathione reductase (GR). There were five pAOPs that were found to be involved in 'actin cytoskeleton reorganization' and 'actin filament severing', respectively. Many cytoskeletal proteins are sensitive to ROS and are critical for the function of vascular cells, serving mechanical, organizational, and signaling roles [32], suggesting these AOPs might play important roles in redox regulation of actin cytoskeleton. In addition, most of the AOPs enriched in cellular components were involved in the term 'integral component of membrane'. Four AOPs were located in 'endoplasmic reticulum membrane (ERM)'. Furthermore, four AOPs were classified in the molecular function term of 'metal ion transmembrane transporter activity'. A previous study showed that antioxidants reacting with ROS, RNS, and radicals were produced in association with electron transport in the endoplasmic reticulum [33], consistent with our results. Moreover, several AOPs were found to be involved in the biological processes of 'transporter' and 'FAD/FMN binding'. With the response to oxidative stresses, mitochondria would produce a limited amount of ROS to initiate a molecular stress response by inducing defense enzymes such as SOD and catalase, as well as other stress defense pathways [34]. Thus, these AOPs might reduce oxidative stress via mitochondrial regulation. Moreover, AOPs would exert effective metal ion (Fe 2+ /Cu 2+ ) chelating activity to defend againstoxidative stress.
The differential expression profiles of the 1008 pAOPs matched genes were analyzed. Similarly, the numbers of pDAPs in ST vs. MF was the highest, secondly in IL vs. ST, and the least in MF vs. YF, indicating that the differences primarily existed in different growth stages. Furthermore, the highly differentiated expressed pDAPs matched genes between stages were selected. The result showed that the highest difference existed in the comparisons of growth stages, suggesting that these pDAPs genes might be play crucial roles in the mummification of the infected larvae and fruit-body maturation. For example, the AOP matched catalase (MSTRG.8461) was found to be up-regulated in the highest level in ST compared to that in IL and MF. Catalase is one of the crucial antioxidant enzymes that mitigates oxidative stress to a considerable extent by destroying cellular H 2 O 2 to produce H 2 O and O 2 . To avoid the infection of fungal pathogens, insect hosts often rapidly produce plenty of ROS to directly kill pathogen; the increase of catalase would detoxify ROS for host infection [8,35]. Most of oxygen-dependent oxidoreductases use FAD as a co-factor. A previous study showed that a large FAD-linked oxidase encoded by RSc0454 is required for pathogenicity [36]. In the present study, four AOPs were identified to match FAD linked oxidase (MSTRG.3794) and were remarkably up-regulated in ST compared to that in the other two stages, suggesting that the AOPs might be related to this fungal pathogenicity in O. sinensis. Moreover, three AOPs from MSTRG.10756, a zinc finger protein GIS2 OS, were remarkably up-regulated in MF compared to ST. C 2 H 2 -type zinc finger proteins (ZFPs) are thought to play important roles in modulating the responses of plants to drought, salinity, and oxidative stress, and different members have different roles during stresses. Moreover, the mutant of ste A, a transcription factor homeodomain C 2 H 2 zinc finger protein was found to not form cleistothecia in A. nidulans, suggesting that these AOPs might be related to the adaptation to the high latitude and sexual development [37]. Serine/threonine-protein kinase SKY1(MSTRG.7638) has a much higher expression level in IL and MF compared to that in ST. MAP kinases belong to the family of serine/threonine protein kinases and are activated by a MAPKKK-MAPKK-MAP kinase cascade, which plays critical roles in pathogenicity and in fungal development [38]. It was suggested that the ROS signal and MAPK cascade might be cross-linked and crucial for the fungal pathogenicity and sexual development in O. sinensis. Sexual differentiation process protein isp7 (MSTRG.6754) has a much higher expression level in MF than in ST (Table 3) and could be related to sexual reproduction in this fungus. Generally, these differentially expressed antioxidant peptides/proteins would be induced by oxidative stresses during the periods of fungal-insect interaction and sexual development and participate in regulating fruit-body development, the mechanism of which need to be further studied.
A total of 1655 AOPs was identified in our proteomic data. For examples, alnexin (CNX)(MSTRG.2357) is an integral membrane protein of the endoplasmic reticulum, and it is also associated with mitochondria oxidative stresses [39,40]. Here, some AOPs were matched to a putative peroxiredoxin (MSTRG.2560, IOZ07G2061; Table 3). Peroxiredoxins are central elements of the antioxidant defense system and the dithiol-disulfide redox regulatory network of plant and cyanobacterial cells [41]. They serve in the context of photosynthesis and respiration, but also in metabolism and development of all tissues, in nodules as well as during seed and fruit development [42]. Peroxiredoxins were firstly identified in yeast [43]. Here, we firstly identified the putative AOPs matched to peroxiredoxin, which might serve in the fruit-body development and defense against oxidative stresses in O. sinensis. The AOPs (IOZ07G1895) matched to SOD [Cu-Zn], SOD, and CAT play a critical role in scavenging free radicals and therefore are used as indicators to evaluate the effect on the antioxidant defense system [44]. Previous studies revealed that SOD was involved in cold response and fruiting body development [10].
To conclude, the integration of transcriptome-and proteome-seqs to analyze the pAOPs in O. sinensis from this study pave the way for developing antioxidant peptides in this medicinal fungus and understanding the biology underlying fungal pathogenicity and the fruit-body development at high latitudes. : Table S1. The assembled result of three biological replicates in each sample from three growth stages in the cultivated O. sinensis and the wild specimens by using String Tie; Table S2. Differentially expressed genes between stages in O. sinensis; Table S3. The information of all identified ORF and smORF;  Table S10. A number of spectra, peptide, and protein identified by using ITraq approach from three adjacent growth stages (ST, PR, and MF) in the cultivated O. sinensis and the wild specimens; Table S11. The pAOPs identified in proteome data.

Supplementary Materials
Author Contributions: X.T. designed and performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. J.G. participated in designing the study, provided the samples, and reviewed drafts of the paper. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: The following information was supplied regarding data availability: Raw data of RNA-seqs are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive under the accessions GSE160504. Raw proteomics data was deposited in integrated proteome resources (iProX) with the primary accession code PXD030687. The O. sinensis strain is available in the China General Microbiological Culture Collection Center, accession number CGMCC 3.14243.