All strains, plasmids and primers used in this study are listed in Supplementary Table S1. P. acidipropionici CGMCC1.2232 was propagated anaerobically at 30°C in culture medium containing the following: yeast extract 5 g⋅L^–1; peptone 5 g⋅L^–1; dipotassium hydrogen phosphate 0.25 g⋅L^–1; manganese sulfate 0.05 g⋅L^–1; glucose 5 g⋅L^–1. Resazurin was added to a final concentration of 0.05 % as an oxygen indicator. Escherichia coli DH5α was used for cloning and plasmid propagation, and was grown aerobically in Luria-Bertani (LB) medium at 37°C. E. coli S17–1 (λ pir) was used as the donor strain to transform P. acidipropionici via conjugation. All media were sterilized through autoclaving (121°C for 20 min), after which 50 μg⋅mL^–1 ampicillin, kanamycin, and erythromycin or 25 μg⋅mL^–1 spectinomycin was added when necessary.

The chromosomal target genes were deleted seamlessly according to the method we reported previously (Jiang et al., 2015). The upstream and downstream regions flanking the lactate dehydrogenase (ldh) and pyruvate oxidase (poxB) genes were amplified using the genome of P. acidipropionici as the template. The resulting products with the 2.0 kb Smr/Spcr Ω cassette were amplified using the high copy vector pBluescriptII SK+ (Stratagene). The Kan^R-suicide vector pJP5603 was selected to construct the recombinant plasmids pJLDH and pJPOXB that were transferred into P. acidipropionici cells by conjugation using the donor strain E. coli S17–1 (λ pir). The overexpression of genes was conducted as described before (Jiang et al., 2015). The methylmalonyl-CoA carboxyltransferase (mmc) gene was amplified from the genome of P. acidipropionici using the primers mmc-for1 and mmc-rev1. The shuttle vector pBRESP36A was selected to construct the mmc overexpression plasmid pBMMC that was transferred into P. acidipropionici via conjugation.

Evolutionary engineering of P. acidipropionici mutants was carried out using multiple stress in the form of different pH values (7, 6, 5, 4), glucose concentrations (60, 90, 120 g⋅L^–1) and oxygen flux values (99.999% N₂, 79% N₂ + 11% CO₂ + 10% O₂, 79% N₂ + 21% O₂) during anaerobic cultivation in a 2 L NBS fermenter at 30 °C and 120 rpm (Siero et al., 2015; Wu et al., 2017). Samples were taken every 24 h to assess cell growth and stress tolerance. Fresh medium was added to maintain a normal 2 L continuous culture system. multiple stress was applied stepwise to maintain the cells at an OD₆₀₀ above 1.0. The final mutant was then used as the experimental group, and the wild-type was considered as the control group. The two groups were cultured under single stress (pH, glucose concentration, oxygen flux) to evaluate cell growth and compare their stress tolerance.

For PA fermentation, strains were pre-cultured in 100 mL anaerobic bottles with 50 mL medium with different pH, glucose concentrations and oxygen flux until the OD₆₀₀ reached 2.0. Then, 5% of seed liquid was used to inoculate a 2-L NBS fermenter with corresponding medium, followed by anaerobic culture at 30°C and 120 rpm for batch fermentations. Samples were taken regularly to measure the concentrations of LA, AA and PA. Additionally, the wild type and the evolved mutants were subjected to transcriptomic and proteomic analyses to further investigate the tolerance mechanism of P. acidipropionici.

Total RNA of the wild type and the mutant III obtained via evolutionary engineering from P. acidipropionici-Δldh-ΔpoxB+mmc were extracted using an RNeasy Mini Kit (Qiagen, Hilden, Germany). The concentration and purity of the RNA was measured using a NanoDrop instrument. Then, mRNA was isolated and purified using the Poly (A) Purist^TM MAG Kit (Ambion, United States). The cDNA of the transcriptomes was synthesized using the SMARTer^TM PCR cDNA Synthesis Kit (Clontech, United States). The NEBNext^® mRNA Library Prep Reagent Set for Illumina^® (NEB, UAS) was performed to construct a cDNA library (∼300 bp), which was subsequently sequenced on the Illumina Xten platform (PE150 mode).

The clean data collected after filtering was aligned to the reference transcript sequences using SOAP2. The alignment was conducted via RPKM conversion to obtain the expression level of the transcript.

Differential expression genes were analyzed using the R package edgeR, and the screening threshold was p-value < 0.05, log FC (fold change (condition 2/condition 1) for a gene) ≥ 1 or log FC ≤ −1.

Gene Ontology analysis of the DEGs was performed using a hyper-geometric distribution. Each GO term with false discovery rate (FDR) ≤ 0.05 was chosen as a significant enriched GO entry. Significant enrichment by pathway can determine the most important biochemical metabolic pathways and signal transduction pathways related to the DEGs. KEGG (Kyoto Encyclopedia of Genes and Genomes) classification was analyzed using R software by setting the parameter-fdr to BH to find pathways (FDR ≤ 0.05) where differential genes are significantly enriched relative to all annotated genes.

where N is the number of genes with pathway annotation; n is the number of differentially expressed genes in N; M is the number of genes annotated as belonging to a particular pathway; m is the number of differentially expressed genes in M.

A high strength ultrasonic processor (Scientz, China) was used to treat the samples on ice three times in the lysis buffer (8 M urea, 1% protease inhibitor mixture). The cell debris was removed by centrifugation at 12,000 g and 4 °C for 10 min, and the supernatant was collected. The protein concentration was determined using a BCA kit according to the manufacturer’s instructions. The peptide obtained by trypsin digestion was reconstituted in 0.5 M Triethylammonium Bicarbonate (TEAB) and processed by a TMT kit/iTRAQ kit (AB Sciex, Foster City, CA, United States). The tryptic peptides resuspended in 0.1% formic acid (buffer A) were directly loaded onto a self-made reversed-phase analytical column (15 cm length, 75 μm inner diameter). The gradient consisted of 6-23% buffer B (0.1% formic acid in 98% acetonitrile) over 26 min, 23–35% in 8 min, 35–80% in 3 min, then hold at 80% for the last 3 min before reversing to buffer A for re-equilibration. The flow rate was set at constant 400 nL⋅min^–1 on an EASY-nLC 1000 UPLC system.

The peptides were exposed to a nano-spray ionization (NSI) source followed by tandem mass spectrometry (MS/MS) in a Q Exactive^TM Plus instrument (Thermo Fisher Scientific, United States) coupled online to the UPLC. The applied electrospray voltage was 2 kV. The m/z scan range was from 350 to 1,800 for full scan, and whole peptides were detected at 70,000 resolution in the OrbiTrap. Peptides were then selected for MS/MS using an NCE setting of 28, and the fragments were detected at 70,000 resolution in the OrbiTrap. A data-dependent program was performed alternately after an MS scan, followed by 20 MS/MS scans, and eliminated dynamically after 15 s. Automatic gain control (AGC) was set at 5E4. The first fixed mass was set as 100 m/z.

The MaxQuant search engine (v.1.5.2.8) was used to analyze the obtained MS/MS data. Tandem mass spectra were searched against the UniProt Proteome database of P. acidipropionici concatenated with a reverse decoy database. Trypsin/P was designated as the enzyme, allowing up to two missed cleavage sites. The mass tolerance of precursor ions was set to 20 ppm in the first retrieval, 5 ppm in the main retrieval and 0.02 Da for fragment ions. Oxidation of methionines was designated as variable modification, while carbamidomethylation of cysteines was designated as a fixed modification. The FDR for modification sites, peptides and proteins was adjusted to < 1% and the lowest score for peptides was set at > 40. The DEPs were selected when the difference ratio of mutant III/wild type was greater than 1.5 (significantly up-regulated proteins) or less than 1.15 (significantly down-regulated proteins) with a t-test p-value < 0.05.

The proteins were classified into three categories: biological process, cellular component and molecular function by GO annotation¹. Two-tailed Fisher’s exact test was used to test the enrichment of the DEPs based on all identified proteins from the transcriptome of P. acidipropionici. The p-value from the test was subjected to negative logarithmic (−log10) conversion. Moreover, KEGG² analysis was performed to predict enriched pathways in which the DEPs were involved using the uniform two-tailed Fisher’s exact test. KEGG mapper was used to map the annotation results onto the KEGG pathway database. These pathways were classified into hierarchical categories according to the KEGG website. Items with a corrected p-value < 0.05 were considered significant.

The qRT-PCR experiment was carried out in QuantStudio 3 and 5 Real-Time PCR Systems and Applied Biosystems^TM TaqMan^TM* (Thermo Scientific, Wilmington, DE, United States) as previously reported (Guan et al., 2014). Cells of wild type and the mutant III were cultured overnight and then transferred into fresh liquid medium to control OD₆₀₀ at ∼0.6. Total cellular RNA was extracted by using the RNeasy Mini Kit (Qiagen, Germany). cDNA was synthesized using a PrimeScript RT reagent Kit With gDNA Eraser (Takara, Dalian, China). The primers of verification genes were designed by using Primer Express 3.0 (Supplementary Table S2) and followed with PCR amplification (pre-incubation at 95°C for 30 s; 40 cycles at 95°C for 5 s, 60°C for 20 s, and cooling at 50°C for 30 s). The 20.0 μL reaction system was composed of 10.0 μL SYBR Premix Ex Taq (2×), 1.2 μL DNA template, 0.4 μL forward and reverse primers (10 μM), and 8 μL dH₂O. The rimM, encoding 16S rRNA processing protein, was selected as a reference gene for normalization.

The optical density at 600 nm (OD₆₀₀) was measured using a standard spectrophotometer (Ultrospec 3300 pro, Amersham Bioscience) to analyze cell growth. The concentrations of organic acids were determined by high-performance liquid chromatography (HPLC) on an Aminex HPX-87H ion-exclusion column (Bio-Rad, United States, 9 μm × 300 mm × 7.8 mm) and ERC 7515A refractive index detector (ERC, Saitama, Japan). 10 mM H₂SO₄ solution at a flow rate of 0.4 mL/min was used as a mobile phase. The column temperature was set at 55°C and refractometer temperature was set at 30°C. Samples were centrifuged at 12,000 g for 5 min followed by dilution in ultrapure water and then boiled for 15 min. The supernatant obtained by centrifugation at 13,000 g for 15 min was filtered through 0.22 μm membrane pore size (diameter, 25 mm) for analysis.

The ldh and poxB genes, which are main contributors of the generation of the by-products acetate and lactate (Liu et al., 2016), were deleted. Additionally, the mmc gene from P. acidipropionici, which controls the carbon flow into the propionate-producing Wood-Werkman cycle (Wang et al., 2015b), was overexpressed. These manipulations yielded the strains P. acidipropionici-Δldh, P. acidipropionici-Δldh-ΔpoxB, and P. acidipropionici-Δldh-ΔpoxB+mmc, respectively. In evolutionary engineering, specific stress is applied gradually on the basis of the growth rate of the adapted strain, which gives a buffer period so that the strain can adequately mobilize the in vivo response mechanisms to resist stress (Figure 1). During the evolution of P. acidipropionici-Δldh into the mutant I, the growth of the experimental strain displayed a peak at OD₆₀₀ 17.1 on the fourth day of culture, followed by a downward trend because of the decline of pH and the increase of glucose concentration. Subsequently, the growth rate of P. acidipropionici-Δldh tended to be stable at OD₆₀₀ 13.9 and began to descend when the pH was decreasing again (6 to 5), and was ultimately stable at OD₆₀₀ 15.9. From P. acidipropionici-Δldh-ΔpoxB to the mutant II (deletion of poxB followed by evolutionary engineering), the formation of the subsequent two peaks at OD₆₀₀ 16.5 and 15.4 on the 6th and 15th day corresponded to the progressive deepening of the levels of single stress factors. From P. acidipropionici-Δldh-ΔpoxB+mmc to the mutant III (overexpression of mmc followed by evolutionary engineering), the growth of the strains started to pick up a steady upward trend, and the strain showed a normal growth curve after supplementation with fresh medium. The highest OD₆₀₀ was 14.6, while the lowest OD₆₀₀ of 9.8 was observed with pH 4, 120 g⋅L^–1 glucose and 21% oxygen sparging. These results indicate that the experimental strain could still maintain a considerable growth rate under the superimposed multiple stress, even when exposed to simultaneously increasing stress levels. This indicated that the strain had developed a certain level of ‘multiple-stress’ protection and acquired a fitness advantage based on the gradual application of stress by artificial control so that it could deal with multiple stress.

To further investigate the growth rate of the strain obtained by metabolic and evolutionary engineering under multiple stress, the wild type and mutant III were compared. Compared with the 0.09 h^–1 specific growth rate of the wild type under normal fermentation conditions (pH 7), mutant III had a significant improvement at 0.14 h^–1, which was apparently higher than the 0.131 ± 0.007 h^–1 previously reported for growth on soy molasses and soy molasses hydrolysate (Yang et al., 2018). However, when the strains were exposed to low pH, the growth of both the wild type and the mutant III was suppressed accordingly. Nevertheless, the growth of the mutant III was significantly superior to that of the wild type. With the decrease of pH, the growth rate difference between the wild type and mutant III gradually increased, and mutant III could maintain normal growth, while the growth of the wild type was clearly suppressed (Supplementary Figure S1A). As the concentration of glucose in the medium increased, the difference of the growth rate between the wild type and mutant III gradually increased. At normal glucose concentrations (30 g⋅L^–1 glucose), the specific growth rate of the wild type was about 0.09 h^–1, while that of mutant III reached 0.2 h^–1, which was 1.4-fold higher than the previously reported value for the same glucose concentrations (Wang et al., 2015a). While the wild type almost stopped growing at a high glucose concentration of 120 g⋅L^–1, mutant III could maintain growth at OD₆₀₀ 0.4 (Supplementary Figure S1B). Furthermore, the strain exhibited similar growth characteristics when exposed to oxygen stress compared with acid stress and osmotic stress. A 2-fold growth rate difference has been observed under anaerobic fermentation conditions (oxygen content 0.001%). Under these conditions, the specific growth rate of mutant III reached 0.18 h^–1, while that of the wild type was 0.09 h^–1. By contrast, with 21% oxygen content in the sparged gas mixture, the specific growth rate of the wild type was only 0.03 h^–1, while that of mutant III was 0.05 h^–1 (Supplementary Figure S1C). As a result, the growth performance of mutant III under multiple stress was clearly superior to that of the wild type. Importantly, the generation of the mutant III has proved that evolutionary engineering is an effective method for improving the resistance of P. acidipropionici to multiple stress.

Batch fermentations were carried out using wild-type P. acidipropionici as well as the mutants I, II and III to investigate their PA production capacity. As shown in Table 1, under the optimal culture conditions (pH 7, 99% N₂, 5 g⋅L^–1 glucose), the PA titer and productivity of mutant I were increased by 12.4 and 13.5% compared to the wild type. Its performance was therefore also significantly better than that of P. jensenii-Δldh with 1.6 and 8.1% enhancement (Liu et al., 2016), respectively. At the same time, a 68.8% of decrease of LA yield (3.06 ± 0.82 g⋅L^–1 vs. 1.04 ± 0.17 g⋅L^–1) was obtained. However, there was a 3% increase of AA titer (5.24 ± 0.83 g⋅L^–1 vs. 5.63 ± 0.62 g⋅L^–1), which may be related to the deletion of ldh, which increased carbon flow from pyruvate into the synthesis of LA (Luna-Flores et al., 2018). The PA production of the mutant II only increased by 3.2%, while the titers of AA and LA were both decreased by more than 70%, which was consistent with the variation tendency of the PA, LA and AA titers of P. acidipropionici-ΔpoxB-Δldh (Guan et al., 2018). A further 37.1% increase of PA titer (28.1 ± 0.96 g⋅L^–1 vs. 38.7 ± 1.14 g⋅L^–1) and 37.8% increase of PA productivity (0.216 ± 0.006 g⋅L^–1 vs. 0.298 ± 0.008 g⋅L^–1⋅h^–1) was achieved in the mutant III. This increase was much more significant than the reported value achieved via the overexpression of mmc in P. freudenreichii in batch fermentation (Wang et al., 2015b). At the same time, there was a 78.6 and 87.8% decrease of the byproduct yields because of the deletion of ldh and poxB, and the overexpression of mmc drove the carbon flow into the synthesis of PA to the greatest extent.

To validate the transcriptome and proteome data, 16 genes involved in pyruvate metabolism with relatively large differences in expression between the wild type and the mutant III were analyzed by qRT-PCR. As can be seen in Figure 2, the expression levels of the tested genes were only slightly different from the results of the omics analyses. However, the 2.52-fold up-regulation expression of sucA obtained from qRT-PCR was a considerably higher than that in the proteomics data, while the expression level of porA with 2.41-fold up-regulation through qRT-PCR analysis was lower than in the proteomics data. Moreover, the 2.51- and 3.02-fold respective down-regulations of fabI and fabD, was consistent with the expression at the transcription and protein level. Overall, the results of qRT-PCR for the 16 selected genes verified the reliability of the transcriptomic and proteomic data.

Pyruvate from glucose or glycerol is a vital precursor in the metabolic pathway of propionic acid synthesis in Propionibacterium spp. In this study, we found that pyruvate metabolism in the mutant III was both significantly enriched from the analyses of combination of transcriptomics and proteomics compared to the wild type. A variety of metabolic pathways are involved in pyruvate metabolism, including the glycolysis pathway, fatty acid biosynthesis, TCA cycle, mevalonate pathway and so on (Fernandez et al., 2008). The differentially expressed genes and proteins related to pyruvate metabolism were both identified. As shown in Figure 3, the pfkA gene encoding 6-phosphofructokinase, a key gene associated with glycolysis pathway was 3.62-fold down-regulated in the mutant III relative to the wild type (Okada et al., 2015). However, at the protein level, the expression of pfkA gene has not been identified but the enzyme aldo that converts 1.6-diphosphate-D-fructose to 3-phosphate-D-glyceraldehyde increased 2.24-fold. Another two key enzymes enolase and pyruvate phosphate dikinase catalyzing the formation of phosphoenolpyruvate from glycerate-2P in the glycolysis pathway was also observed to be significantly up-regulated (Virmani et al., 2019). Generally, the absence of mRNA-protein correlation manifests that the relation between mRNA and protein is not strictly linear, but has a more inherent and complex dependence (Wang et al., 2013). The inequivalence of the identified genes in the transcriptome and proteome may be related to the dynamic regulation of cells to maintain metabolic balance. As a result, the significant up-regulation of expression of the three key genes involved in the glycolysis pathway at the protein level facilitated the considerable accumulation of ATP and pyruvate under accelerated metabolic flow from phosphoenolpyruvate, which provided abundant precursors and energy for pyruvate metabolism.

The DEGs and DEPs in the mutant III involved in TCA cycle were also investigated (Figure 3). At the transcript level, the genes sdhA and frdA were up-regulated 2.03- and 1.61-fold, compared to 2.94- and 2.93-fold up-regulation at the protein level in the mutant III. The obviously higher expression level in the proteome indicated that high protein abundance was a prerequisite for the mutant III to confront environment and DNA damage stress. Additionally, the confirmed DEGs fumA and fumB, DEPs such as aceE, porA, gltA, sdhA, frdA and so on all showed obvious up-regulation relative to the wild type, which implied that the TCA cycle in the mutant III was highly active. Interestingly, the quantitative superiority of DEPs compared with DEGs suggested that gene expression profiles were mainly highlighted in the form of proteome in P. acidipropionici after multiple stressors and genetic manipulation. The results showed a poor correlation between the quantities of these genes, which indicated a complex regulatory mechanism controlling expression both at the RNA and the protein levels.

Fatty acids are excellent sources of energy and precursors for phospholipids, sterols, sphingolipids, as secondary metabolites and signaling molecules (Janßen and Steinbüchel, 2014). In this study, we paid particular attention to the differentially expressed genes and proteins related to the most enriched fatty acid biosynthesis pathway. As can be seen in Figure 3, a total of five genes participating in fatty acid biosynthesis were down-regulated both at the mRNA and protein level to different degrees. At the RNA-Seq and proteome level, the essential gene fabD encoding malonyl-CoA:ACP transacylase for fatty acid biosynthesis was found to have 3.79-fold and 1.60-fold down-regulation, respectively. Malonyl-CoA:ACP transacylase can catalyze the formation of malonyl-ACP, which promotes fatty acid neogenesis and fatty acid chain elongation (Janßen and Steinbüchel, 2014). The down-regulated expression of malonyl-CoA:ACP transacylase indicated that the metabolic flux toward the synthesis of malonyl-ACP was repressed. The gene fabH encoding β-ketoacyl-ACP synthase III that catalyzed malonyl-ACP and acetyl-CoA units in the initial elongation step was also down-regulated 3.50- and 1.55-fold, which further restricted the metabolic flux from malonyl-ACP and led to the low accumulation of β-Ketoacyl-ACP (Handke et al., 2011). Remarkably, β-ketoacyl-ACP reductase encoded by fabG, which promotes the production of β-hydroxyacyl-ACP at the expense of NADPH was found to be as much as 4.17- and 1.65-fold down-regulated after the combination of metabolic and evolutionary engineering. Compared to the significant down-regulation at the RNA-Seq level (3-∼4-fold), the proteome analysis revealed lower expression profiles (1-∼2-fold), which indicated the down-regulation of the fatty acid biosynthesis process may be slowed down. Since the generation of malonyl-CoA requires ATP, the viability of microbes adopted to withstand multiple stresses and genetic manipulation is highly dependent on ATP. Based on this, the ATP in the mutant III seemed to be economically consumed for fatty acid biosynthesis to reduce the pernicious effects of environmental and DNA damage stress. The low consumption of the important precursors pyruvate and ATP, indicated that the ATP accumulation and metabolic flux to PA synthesis may be increased during fatty acid biosynthesis. Therefore, we hypothesized that the repressed metabolic flux toward fatty acid synthesis may be associated with the high yield of target product in this study.