Abstract Lactate production is enhanced by adding calcium carbonate or sodium hydroxide during fermentation. However, Bacillus coagulans 2-6 can produce more than 180 g/L L-lactic acid when calcium lactate is accumulated, but less than 120 g/L L-lactic acid when sodium lactate is formed. The molecular mechanisms by which B. coagulans responds to calcium lactate and sodium lactate remain unclear. In this study, comparative transcriptomic methods based on high-throughput RNA sequencing were applied to study gene expression changes in B. coagulans 2-6 cultured in non-stress, sodium lactate stress and calcium lactate stress conditions. Gene expression profiling identified 712 and 1213 significantly regulated genes in response to calcium lactate stress and sodium lactate stress, respectively. Gene ontology assignments of the differentially expressed genes were performed. KEGG pathway enrichment analysis revealed that ‘ATP-binding cassette transporters’ were significantly affected by calcium lactate stress, and ‘amino sugar and nucleotide sugar metabolism’ was significantly affected by sodium lactate stress. It was also found that lactate fermentation was less affected by calcium lactate stress than by sodium lactate stress. Sodium lactate stress had negative effect on the expression of ‘glycolysis/gluconeogenesis’ genes but positive effect on the expression of ‘citrate cycle (TCA cycle)’ genes. However, calcium lactate stress had positive influence on the expression of ‘glycolysis/gluconeogenesis’ genes and had minor influence on ‘citrate cycle (TCA cycle)’ genes. Thus, our findings offer new insights into the responses of B. coagulans to different lactate stresses. Notably, our RNA-seq dataset constitute a robust database for investigating the functions of genes induced by lactate stress in the future and identify potential targets for genetic engineering to further improve L-lactic acid production by B. coagulans. Introduction Lactic acid has wide applications in the industries of food, medicine, brewing, leather, textiles, environmental protection and agriculture [[35]1–[36]3]. The major industrial application of lactic acid is the production of polylactic acid, a biodegradable polymer. Given the biocompatibility, biodegradability and excellent mechanical properties of polylactic acid, it is considered the most promising novel packaging material of this century. It can replace non-degradable petrochemical materials such as polyethylene, polypropylene and polystyrene. Importantly, widespread use of polylactic acid has the potential to solve ‘white pollution’ and mitigate the current energy crisis [[37]2, [38]4]. During fermentative production of lactic acid, a neutraliser such as calcium carbonate or sodium hydroxide is required to maintain neutral or mildly acidic conditions for the fermentation liquid. Sodium lactate, calcium lactate or another lactate can accumulate in the fermentation broth depending on the type of neutraliser selected. High concentrations of lactate result in stress effects on lactic acid producers, and are among the key factors limiting the enhancement of lactic acid yields [[39]5]. Previous studies indicated that the maximum yield of lactic acid can exceed 200 g/L when calcium lactate accumulates in the fermentation broth [[40]6, [41]7]; however, the lactic acid yield is generally less than 125 g/L when sodium lactate is formed [[42]5, [43]8]. Therefore, the influences of lactate on lactic acid production may be different depending on the type of lactic acid formed. Bacillus coagulans strains are L-lactic acid producers with industrial potential because they can tolerate high fermentation temperatures, their culture media do not need sterilisation and their products exhibit high optical purity [[44]9, [45]10]. B. coagulans 2–6 can produce more than 180 g/L L-lactic acid with high optical purity when calcium lactate is accumulated, but less than 120 g/L L-lactic acid when sodium lactate is formed [[46]5, [47]11]. Under the same concentration of lactates, sodium lactate results in stronger stress on B. coagulans 2–6 compared with calcium lactate [[48]12]. However, the mechanisms by which B. coagulans strains respond to different lactates remain unclear. Omics tools are effective for studying the response mechanisms of bacteria to stress conditions. Wang et al. [[49]12] reported a comparative proteomic analysis of B. coagulans 2–6 in response to lactate stress and identified four highly expressed proteins under calcium lactate stress. In the present study, we conducted comparative transcriptomics to study differential expression of genes among B. coagulans 2–6 cultured in normal, sodium lactate stress and calcium lactate stress conditions. Similarities and differences in the response of B. coagulans 2–6 to sodium lactate and calcium lactate were also investigated. Materials and Methods Bacterial strains and cultivation conditions The B. coagulans wild-type strain 2–6 [[50]11, [51]13] was used in this study. It was cultivated at 50°C with shaking at 150 rpm in non-stress medium (GY medium: 40 g/L glucose, 10 g/L yeast extract, 1 g/L CaCl[2], pH 6.7) until the cell dry weight (CDW) reached 2.0. Then 5 ml bacterial liquid was inoculated into 100 ml fresh non-stress medium and lactate stress medium (GY medium with 0.5 M sodium lactate or 0.25 M calcium lactate), and incubated at 50°C with shaking at 150 rpm until CDW reached 0.4 g/L at late exponential stage of the cell growth as reported before [[52]12]. Cells were harvested for preparation of total RNA. Library preparation and RNA sequencing (RNA-seq) Total RNA of B. coagulans 2–6 was isolated using TRIzol reagent according to the manufacturer’s protocol (Invitrogen, Carlsbad, CA, USA) and then treated with DNase I. The quality and quantity of purified RNA were determined by measuring the absorbance at 260 nm/280 nm (A260/A280) using a Nanodrop ND-1000 spectrophotometer (LabTech, Wilmington, MA, USA). RNA integrity was further verified by electrophoresis through 1.5% (w/v) agarose gel. mRNA was enriched from total RNA using a Ribo-Zero Magnetic Kit (Gram-Positive Bacteria) (EpiCentre, Madison, WY, USA). Synthesis of double-stranded cDNA from the extracted mRNA was done using the Truseq RNA sample prep kit (Illumina, San Diego, CA, USA), followed by 15 cycles of enrichment PCR. Sequencing was carried out on the Illumina HiSeq 2000 platform (Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China). All sequencing data have been deposited in the NCBI Sequence Read Archive under accession numbers SRX700697 (GY), SRX700698 (CA) and SRX 700710 (NA). Transcriptome analysis The raw data were initially processed to obtain clean reads by removing the adapter sequences and low quality bases. The clean reads were then aligned to the reference genome sequence of B. coagulans 2–6 (accession number: [53]NC_015634) using the short sequence alignment software Bowtie 2 [[54]14]. To identify differentially expressed genes (DEGs) between the two different treatment conditions, the expression level for each transcript was calculated by RSEM ([55]http://deweylab.biostat.wisc.edu/rsem/) using the fragments per kilobase of exon per million mapped reads method [[56]15]. The R statistical package software EdgeR ([57]http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html) was employed to quantify differential gene expression [[58]16]. Genes were defined as differentially expressed if they exhibited two-fold or greater change between the samples under a false discovery rate (FDR) of 5% or less. In addition, functional enrichment analyses, including gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG), were performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways compared with the whole-transcriptome background using the formula described in previous studies [[59]17, [60]18]. A Bonferroni-corrected P-value ≤ 0.05 was considered statistically significant. Quantitative real-time PCR (qRT-PCR) verification Five genes belonging to 5 different operons were chosen for confirmation of RNA-seq data by qRT-PCR. Total RNA was extracted using the EasyPure RNA Purification Kit (TRANSGEN BIOTECH, Beijing, China). Then, the total RNA sample was treated with DNase and qRT-PCR was performed on RNA samples to confirm the absence of DNA contamination using 16S primers [[61]19]. qRT-PCR was performed using the TransScript II Green One-Step qRT-PCR SuperMix (TRANSGEN BIOTECH, Beijing, China) and the Corbett Rotor Gene 3000 (Corbett Robotics, Sydney, Australia) according to the manufacturer’s instructions. The genes and primers used for qRT-PCR are shown in [62]S1 Table. The relative gene expression data were analyzed using the 2^-ΔΔCt method as described by Tan et al. [[63]19]. Rotor-Gene Real-Time Analysis Software 6.0 was used for data analysis (Corbett Robotics, Sydney, Australia). All qRT-PCR runs were conducted with three biological and three technical replicates. Results Illumina draft reads Large volumes of data were generated with Illumina HiSeq 2000 sequencing of the three B. coagulans cDNA samples, namely, GY (untreated control), CA, (calcium lactate-stressed sample) and NA (sodium lactate-stressed sample). Each cDNA fragment yielded two 101 base pair long paired-end (PE) reads. After removing low quality reads and trimming off adapter sequences, 13,747,126 (GY), 13,464,898 (CA) and 14,974,752 (NA) high-quality, clean PE sequencing reads with a total of 1,334,784,442 (GY), 1,306,406,094 (CA) and 1,456,148,628 (NA) nucleotides were obtained for the three samples. Clean reads were aligned to the reference genome. The mapping rates were 97.9% (GY), 95.9% (CA) and 95.6% (NA). Gene ontology (GO) functional annotation of transcripts GO assignments were used to classify the functions of the B. coagulans genes under lactate stress. Based on sequence homology, 2053 transcript sequences were assigned at least one GO term, including 43 functional groups at the second level ([64]Fig 1). Among these groups, ‘catalytic activity’ (1359 sequences, 66.2%), ‘binding’ (1091 sequences, 53.1%) and ‘single-organism process’ (587 sequences, 28.6%) were the most abundant subcategories for ‘molecular function’. ‘Cell’ (739 sequences, 36.0%), ‘cell part’ (739 sequences, 36.0%) and ‘membrane’ (431 sequences, 21.0%) terms were dominant in ‘cellular component’. ‘Metabolic process’ (1472 sequences, 71.7%), ‘cellular process’ (1396 sequences, 68.0%) and ‘single-organism process’ (587 sequences, 28.6%) were represented the highest under ‘biological process’. Fig 1. Gene Ontology (GO) functional annotation of transcripts. [65]Fig 1 [66]Open in a new tab All 2,053 transcripts were assigned to at least one GO term and were grouped into three main GO categories and 43 groups, 19 groups in biological process, 12 in cellular component, and 12 in molecular function. The right-hand Y-axis represents the number of genes in a sub-category. The left-hand Y-axis indicates the percentage of a specific sub-category of genes in each main category. Differentially expressed genes (DEGs) involved in response to lactate stress in B. coagulans Comparison of gene expression showed that a total of 1577 unigenes were differentially expressed between any two-way comparison of GY, CA and NA (fold changes ≥2 or ≤−2; FDR<0.05), including 712 isogenes between GY and CA (GY vs CA), 702 isogenes between NA and CA (NA vs CA) and 1213 isogenes between NA and GY (NA vs GY). The amounts of significantly up- and down-regulated genes in NA were 543 and 670, respectively, compared with those in GY. The amounts of significantly up- and down-regulated genes in CA were 398 and 314, respectively, compared with those in GY. Moreover, 364 unigenes were significantly differentially expressed in both CA vs GY and NA vs GY; 121 unigenes were significantly differentially expressed in both CA vs GY and NA vs CA; and 381 unigenes were significantly differentially expressed in both NA vs GY and NA vs CA. Only 92 unigenes were significantly differentially expressed in CA vs GY, NA vs GY and NA vs CA. The numbers of all differentially expressed genes among the three groups of GY, CA and NA are illustrated in [67]Fig 2. Fig 2. Differentially expressed genes that are unique or shared among three groups of GY, NA and CA. [68]Fig 2 [69]Open in a new tab NA vs GY refers to the comparison between sodium lactate stressed (NA) and untreated (GY) groups. CA vs GY refers to the comparison between calcium lactate stressed (CA) and GY groups. NA vs CA refers to the comparison between NA and CA groups. Numbers in each section of the figure indicate the numbers of differently expressed genes in the indicated comparison. GO assignments of DEGs GO assignments of DEGs in response to calcium lactate stress and sodium lactate stress were performed. The results are shown in [70]Table 1. The 398 significantly upregulated genes in CA compared with those in GY were assigned 34 GO terms, namely, 14 subcategories of biological process, 9 subcategories of cellular component and 11 subcategories of molecular function. The 314 significantly down-regulated genes in CA compared with those in GY were assigned 34 GO terms, namely, 15 subcategories of biological process, 10 subcategories of cellular component and 9 subcategories of molecular function. The 543 significantly upregulated genes in NA compared with those in GY were assigned 39 GO terms, namely, 17 subcategories of biological process, 12 subcategories of cellular component and 10 subcategories of molecular function. The 670 significantly down-regulated genes in NA compared with those in GY were assigned 40 GO terms, namely, 18 subcategories of biological process, 10 subcategories of cellular component and 12 subcategories of molecular function. The significantly up- or down-regulated genes in the corresponding GO terms are listed in [71]Table 1, [72]S1 Data and [73]S2 Data. Taken together, the results indicate differences in gene expression changes in response to calcium lactate stress and sodium lactate stress in B. coagulans 2–6. Table 1. Gene ontology (GO) assignments of differentially expressed genes. GO type GO term GO id CA vs GY NA vs GY up down up down Biological process biological regulation 0065007 39[74]^* 24 63 61 cellular component organization or biogenesis 0071840 8 8 25 40 cellular process 0009987 161 145 237 359 developmental process 0032502 3 2 14 18 establishment of localization 0051234 30 69 52 105 immune system process 0002376 0 1 0 1 localization 0051179 30 71 52 117 locomotion 0040011 1 4 2 19 metabolic process 0008152 195 149 258 368 multi-organism process 0051704 0 2 1 3 negative regulation of biological process 0048519 2 2 4 3 positive regulation of biological process 0048518 1 1 regulation of biological process 0050789 37 23 62 60 response to stimulus 0050896 34 23 35 50 Signaling 0023052 3 6 7 16 single-organism process 0044699 153 133 201 311 biological adhesion 0022610 1 0 1 1 reproduction 0000003 3 1 Cellular component Cell 0005623 79 83 129 201 cell part 0044464 79 83 129 201 extracellular region 0005576 1 3 1 2 extracellular region part 0044421 0 1 macromolecular complex 0032991 34 29 58 89 Membrane 0016020 36 78 65 131 membrane-enclosed lumen 0031974 1 2 2 3 membrane part 0044425 22 56 46 89 Nucleoid 0009295 1 1 Organelle 0043226 6 10 13 51 organelle part 0044422 2 5 7 26 Virion 0019012 2 0 virion part 0044423 2 0 Molecular function Binding 0005488 145 123 213 266 catalytic activity 0003824 185 143 236 340 electron carrier activity 0009055 8 6 7 2 metallochaperone activity 0016530 0 1 0 1 molecular transducer activity 0060089 3 5 5 14 nucleic acid binding transcription factor activity 0001071 12 7 26 15 receptor activity 0004872 1 2 1 6 structural molecule activity 0005198 4 5 3 31 transporter activity 0005215 20 55 35 83 antioxidant activity 0016209 5 0 2 1 enzyme regulator activity 0030234 2 0 0 1 protein binding transcription factor activity 0000988 3 0 4 3 [75]Open in a new tab *These are numbers of significantly up- or down-regulated genes in the corresponding GO term GO and pathway functional enrichment analysis of DEGs GO function and pathway enrichment analyses were performed on all DEGs. According to GO functional enrichment analysis between CA and GY, the GO term ‘response to stress’ was represented the most and included 96 unigenes and 41 significantly differentially expressed unigenes. Analysis between NA and GY revealed one significantly enriched GO term. This GO term was ‘biological process’ (GO: 0008150), which involved 1813 unigenes and 793 significantly differentially expressed unigenes. No significantly enriched GO term was found between NA and CA. Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis showed that the DEGs between the CA and GY groups were significantly involved in the ‘ATP-binding cassette (ABC) transporters’ (KEGG: ko02010) pathway, which included 96 unigenes and 37 significantly differentially expressed unigenes. The DEGs between the NA and GY groups were significantly involved in the ‘amino sugar and nucleotide sugar metabolism’ pathway, which included 35 unigenes and 26 significantly differentially expressed unigenes. The DEGs between the NA and CA groups were significantly involved in the ‘amino sugar and nucleotide sugar metabolism’ (KEGG: ko00520) pathway, which included 35 unigenes and 22 significantly differentially expressed unigenes. Identification and comparison of DEGs involved in ‘ABC transporters’ ‘ABC transporters’, which are also called traffic ATPases, use the binding and hydrolysis of ATP to power the translocation of a diverse assortment of substrates across biological membranes, ranging from ions to macromolecules [[76]20]. ‘ABC transporters’ either import or export substrates with relative specificity for a given substrate. Dozens of ABC transporters were identified by genome analysis in Escherichia coli [[77]21] and Bacillus subtilis [[78]22]. KEGG pathway enrichment analysis revealed that ‘ABC transporters’ was significantly affected by calcium lactate stress. A total of 10 significantly upregulated genes ([79]Table 2) and 27 significantly down-regulated genes ([80]S2 Table) were involved in ‘ABC transporters’ under calcium lactate stress condition (CA group), compared with those under non-stress condition (GY group). Fourteen significantly upregulated genes ([81]Table 3) and 35 significantly down-regulated genes ([82]S3 Table) were involved in ‘ABC transporters’ under sodium lactate stress condition (NA group) compared with those under non-stress condition (GY group). Six genes, namely, glycine betaine/L-proline ABC transporter ATPase (BCO26_0425), binding-protein-dependent transport system inner membrane protein (BCO26_0553), NMT1/THI5 like domain-containing protein (BCO26_0554), ABC transporter-like protein (BCO26_0555), periplasmic solute binding protein (BCO26_2084) and ABC transporter (BCO26_2523), were significantly upregulated in both the CA and NA groups compared with those in the GY group. Four genes, namely, phosphate binding protein (BCO26_0705), phosphate ABC transporter permease (BCO26_0706, BCO26_0707) and ABC transporter-like protein (BCO26_0743), were significantly up-regulated in CA but not in NA compared with GY. Eight genes, namely binding-protein-dependent transport system inner membrane protein (BCO26_0422), glycine betaine ABC transporter substrate-binding protein (BCO26_0423), binding-protein-dependent transport system inner membrane protein (BCO26_0424), transport system permease (BCO26_0685, BCO26_0686), ABC transporter-like protein (BCO26_0687), periplasmic binding protein (BCO26_0688) and ABC transporter-like protein (BCO26_2588), were significantly upregulated in NA but not in CA compared with GY. The significantly down-regulated genes involved in ‘ABC transporters’ under the two kinds of lactate stress were mostly the same ([83]S2 Table and [84]S3 Table). However, sodium lactate induced more significantly down-regulated genes than calcium lactate (35 to 27) and higher average down-regulation (fold change ranging from 3.72 to 3.06). Table 2. Significantly upregulated genes involved in ‘ABC transporters’ under calcium lactate stress. Gene ID Description FDR Fold change BCO26_0425 glycine betaine/L-proline ABC transporter ATPase 6.01E-03 2.30 BCO26_0553 binding-protein-dependent transport system inner membrane protein 3.04E-03 2.55 BCO26_0554 NMT1/THI5 like domain-containing protein 1.80E-02 2.03 BCO26_0555 ABC transporter-like protein 2.38E-02 1.94 BCO26_0705 phosphate binding protein 2.62E-06 4.16 BCO26_0706 phosphate ABC transporter permease 8.04E-04 2.98 BCO26_0707 phosphate ABC transporter permease 6.70E-03 2.40 BCO26_0743 ABC transporter-like protein 5.19E-03 2.36 BCO26_2084 periplasmic solute binding protein 1.52E-02 2.05 BCO26_2523 ABC transporter 9.44E-03 2.25 [85]Open in a new tab Table 3. Significantly upregulated genes involved in ‘ABC transporters’ under sodium lactate stress. Gene ID Description FDR Fold change BCO26_0422 binding-protein-dependent transport system inner membrane protein 1.00E-02 2.02 BCO26_0423 glycine betaine ABC transporter substrate-binding protein 4.04E-02 1.60 BCO26_0424 binding-protein-dependent transport system inner membrane protein 4.16E-02 1.60 BCO26_0425 glycine betaine/L-proline ABC transporter ATPase 1.12E-02 1.98 BCO26_0553 binding-protein-dependent transport system inner membrane protein 1.15E-02 2.04 BCO26_0554 NMT1/THI5 like domain-containing protein 5.03E-03 2.22 BCO26_0555 ABC transporter-like protein 3.86E-04 2.82 BCO26_0685 transport system permease 5.90E-05 3.20 BCO26_0686 transport system permease 3.95E-04 2.81 BCO26_0687 ABC transporter-like protein 2.20E-03 2.42 BCO26_0688 periplasmic binding protein 1.03E-03 2.59 BCO26_2084 periplasmic solute binding protein 4.48E-02 1.57 BCO26_2523 ABC transporter 6.40E-03 2.21 BCO26_2588 ABC transporter-like protein 4.57E-02 1.94 [86]Open in a new tab Identification and comparison of DEGs involved in ‘amino sugar and nucleotide sugar metabolism’ In nucleotide sugar metabolism, nucleotide sugars act as donors for sugar residues in the glycosylation reactions that produce polysaccharides, which are important constituents of the cell wall [[87]23]. Nucleotide sugars are also intermediates in nucleotide sugar interconversions, which produce some of the activated sugars necessary for glycosylation reactions [[88]24]. KEGG pathway enrichment analysis revealed that ‘amino sugar and nucleotide sugar metabolism’ was significantly affected by sodium lactate stress but not by calcium lactate stress. Fifteen significantly upregulated genes ([89]Table 4) and 11 significantly down-regulated genes ([90]S4 Table) were involved in ‘amino sugar and nucleotide sugar metabolism’ under sodium lactate stress condition (NA group) compared with those under non-stress condition (GY group). Eight significantly upregulated genes and one significantly down-regulated gene were involved in ‘amino sugar and nucleotide sugar metabolism’ under calcium lactate stress condition (CA group) compared with those under non-stress condition (GY group) ([91]Table 5). Seven significantly regulated genes were found to exhibit the same trend under the two types of lactate stress conditions, namely, NAD(P)(+)-protein-arginine ADP-ribosyltransferase (BCO26_0186), glucosamine/fructose-6-phosphate aminotransferase (BCO26_0255), galactokinase (BCO26_0622), UDP-glucose 4-epimerase (BCO26_0623), N-acetylglucosamine-6-phosphate deacetylase (BCO26_0635), glucosamine-6-phosphate isomerise (BCO26_0636) and N-acylglucosamine-6-phosphate 2-epimerase (BCO26_2463). Eight genes were significantly upregulated specifically under sodium lactate stress, namely, UDP-N-acetylglucosamine pyrophosphorylase (BCO26_0036), phosphoglucosamine mutase (BCO26_0177), PTS system glucose subfamily transporter subunit IIA (BCO26_0204), phosphoglucomutase/phosphomannomutase alpha/beta/alpha domain I (BCO26_0507), galactose-1-phosphate uridylyltransferase (BCO26_0624), UDP-N-acetylenolpyruvoylglucosamine reductase (BCO26_1444), ROK family glucokinase (BCO26_1699) and UDP-N-acetylglucosamine 1-carboxyvinyltransferase (BCO26_2818). More genes were significantly down-regulated under sodium lactate stress than under calcium lactate stress (11 to 1). The four subunits of phosphoenolpyruvate-dependent sugar phosphotransferase system (PTS) mannose/fructose/sorbose family transporter (i.e., BCO26_0540, BCO26_0541, BCO26_0542 and BCO26_0543) were significantly down-regulated under sodium lactate stress ([92]S4 Table), whereas the PTS system glucose subfamily transporter subunit IIA (BCO26_0204) was significantly upregulated ([93]Table 4). Table 4. Significantly upregulated genes involved in ‘amino sugar and nucleotide sugar metabolism’ under sodium lactate stress. Gene ID Description FDR Fold change BCO26_0036 UDP-N-acetylglucosamine pyrophosphorylase 3.31E-03 2.29 BCO26_0177 phosphoglucosamine mutase 8.52E-03 2.04 BCO26_0186 NAD(P)(+)-protein-arginine ADP-ribosyltransferase 9.44E-07 4.03 BCO26_0204 PTS system glucose subfamily transporter subunit IIA 9.03E-03 2.16 BCO26_0255 glucosamine/fructose-6-phosphate aminotransferase 2.00E-09 5.17 BCO26_0507 phosphoglucomutase/phosphomannomutase alpha/beta/alpha domain I 3.58E-02 1.63 BCO26_0622 Galactokinase 1.44E-06 3.94 BCO26_0623 UDP-glucose 4-epimerase 3.77E-05 3.28 BCO26_0624 galactose-1-phosphate uridylyltransferase 1.80E-04 2.96 BCO26_0635 N-acetylglucosamine-6-phosphate deacetylase 1.84E-07 4.35 BCO26_0636 glucosamine-6-phosphate isomerase 5.23E-08 4.61 BCO26_1444 UDP-N-acetylenolpyruvoylglucosamine reductase 4.85E-02 1.54 BCO26_1699 ROK family glucokinase 4.25E-04 2.77 BCO26_2463 N-acylglucosamine-6-phosphate 2-epimerase 3.29E-09 5.08 BCO26_2818 UDP-N-acetylglucosamine 1-carboxyvinyltransferase 1.28E-03 2.51 [94]Open in a new tab Table 5. Significantly regulated genes involved in ‘amino sugar and nucleotide sugar metabolism’ under calcium lactate stress. Gene ID Description FDR Fold change BCO26_0186 NAD(P)(+)-protein-arginine ADP-ribosyltransferase 3.19E-02 1.82 BCO26_0255 glucosamine/fructose-6-phosphate aminotransferase 2.02E-02 1.96 BCO26_0622 Galactokinase 1.68E-03 2.64 BCO26_0623 UDP-glucose 4-epimerase 4.24E-02 1.72 BCO26_0635 N-acetylglucosamine-6-phosphate deacetylase 5.19E-03 2.36 BCO26_0636 glucosamine-6-phosphate isomerase 4.01E-02 1.78 BCO26_2126 PTS system glucose subfamily transporter subunit IIA 7.25E-03 2.27 BCO26_2463 N-acylglucosamine-6-phosphate 2-epimerase 1.51E-02 2.06 BCO26_0354 nucleotide sugar dehydrogenase 2.71E-02 -1.87 [95]Open in a new tab Identification and comparison of DEGs involved in ‘glycolysis/gluconeogenesis’ ‘Glycolysis/gluconeogenesis’ is the key pathway for lactic acid production in B. coagulans 2–6. Therefore, we compared DEGs involved in ‘glycolysis/gluconeogenesis’ under lactate stress and non-stress conditions. Thirteen significantly upregulated genes ([96]Table 6) and eight significantly down-regulated genes ([97]S5 Table) were involved in ‘glycolysis/gluconeogenesis’ under sodium lactate stress condition (NA group) compared with those under non-stress condition (GY group). The pyruvate dehydrogenase complex, which mediates the conversion of pyruvic acid to acetyl-CoA, was significantly upregulated under sodium lactate stress compared with non-stress condition. The complex consists of pyruvate dehydrogenase components (BCO26_0962 and BCO26_0963), dihydrolipoamide transacetylase (BCO26_0964) and dihydrolipoamide dehydrogenase (BCO26_0965 and BCO26_1409). Table 6. Significantly upregulated genes involved in ‘glycolysis/gluconeogenesis’ under sodium lactate stress. Gene ID Description FDR Fold change BCO26_0204 PTS system glucose subfamily transporter subunit IIA 1.73E-03 2.90 BCO26_0492 glycoside hydrolase 1.22E-02 2.09 BCO26_0507 phosphoglucomutase/phosphomannomutase alpha/beta/alpha domain I 3.58E-02 1.63 BCO26_0625 aldose 1-epimerase 8.69E-03 2.19 BCO26_0962 pyruvate dehydrogenase E1 component subunit alpha 9.62E-04 2.76 BCO26_0963 transketolase central region 7.31E-04 2.85 BCO26_0964 hypothetical protein BCO26_0964 1.21E-03 2.70 BCO26_0965 dihydrolipoamide dehydrogenase 9.46E-04 2.76 BCO26_1409 dihydrolipoamide dehydrogenase 1.52E-02 2.30 BCO26_1699 ROK family glucokinase 4.72E-02 1.68 BCO26_1901 alcohol dehydrogenase GroES domain-containing protein 5.33E-04 2.95 BCO26_2235 enolase 1.68E-02 1.86 BCO26_2817 fructose-1,6-bisphosphatase 3.96E-03 2.24 [98]Open in a new tab Eleven significantly upregulated genes ([99]Table 7) and five significantly down-regulated genes ([100]S6 Table) were involved in ‘glycolysis/gluconeogenesis’ under calcium lactate stress condition (CA group) compared with those under non-stress condition (GY group) ([101]Table 7). The key genes of the glycolysis pathway, including L-lactate dehydrogenase (BCO26_0531), pyruvate kinase (BCO26_1967), 6-phosphofructokinase (BCO26_1968), enolase (BCO26_2235), phosphoglycerate mutase (BCO26_2237), phosphoglycerate kinase (BCO26_2239) and fructose-bisphosphate aldolase (BCO26_2530), were all significantly upregulated ([102]Table 7). Table 7. Significantly upregulated genes involved in ‘glycolysis/gluconeogenesis’ under calcium lactate stress. Gene ID Description FDR Fold change BCO26_0531 L-lactate dehydrogenase 3.35E-02 1.79 BCO26_1901 alcohol dehydrogenase GroES domain-containing protein 5.52E-03 2.41 BCO26_1967 pyruvate kinase 5.20E-03 2.33 BCO26_1968 6-phosphofructokinase 4.09E-04 2.96 BCO26_2126 PTS system glucose subfamily transporter subunit IIA 7.25E-03 2.27 BCO26_2235 enolase 1.80E-02 1.99 BCO26_2237 phosphoglycerate mutase 1.33E-02 2.08 BCO26_2238 triosephosphate isomerase 2.45E-02 1.89 BCO26_2239 phosphoglycerate kinase 1.75E-02 2.00 BCO26_2240 glyceraldehyde-3-phosphate dehydrogenase, type I 1.37E-02 2.07 BCO26_2530 fructose-bisphosphate aldolase 7.53E-04 2.83 [103]Open in a new tab Under the two types of lactate stress conditions, two significantly upregulated genes involved in ‘glycolysis/gluconeogenesis’ were found, namely, alcohol dehydrogenase GroES domain-containing protein (BCO26_1901) and enolase (BCO26_2235). Four key genes of glycolysis (i.e., iron-containing alcohol dehydrogenase, lactate/malate dehydrogenase, AMP-dependent synthetase and ligase and hypothetical protein BCO26_2297) were significantly down-regulated both in CA and NA compared with those in GY. Pyruvate kinase (BCO26_1967) was significantly down-regulated in NA but upregulated in CA compared with that in GY. Identification and comparison of DEGs involved in ‘citrate cycle (TCA cycle)’ Analysis of the genome of B. coagulans 2–6 revealed that the strain possessed all components of the TCA cycle. The TCA cycle is an important aerobic pathway in the final steps of the oxidation of carbohydrates and fatty acids. As a critical pathway for energy production in cells, it may be beneficial for B. coagulans strains to resist lactate stress conditions. Our previous study indicated that trace amounts of aeration enhance the resistance of B. coagulans to sodium lactate [[104]5]. Therefore, we compared DEGs involved in ‘citrate cycle (TCA cycle)’ under lactate stress and non-stress conditions. Seven significantly upregulated genes and seven significantly down-regulated genes were involved in ‘citrate cycle (TCA cycle)’ under sodium lactate stress condition (NA group) compared with those under non-stress condition (GY group) ([105]Table 8). The pyruvate dehydrogenase complex (BCO26_0962, BCO26_0963, BCO26_0964, BCO26_0965 and BCO26_1409) was significantly upregulated as described above. Aconitate hydratase 1 (BCO26_1352), which responds to the conversion of citrate to isocitrate, and class II fumarate hydratase (BCO26_0741), which responds to the conversion of fumarate to L-malate, were also significantly upregulated. Succinyl-CoA synthetase (BCO26_1117 and BCO26_1118), 2-oxoglutarate dehydrogenase (BCO26_1275 and BCO26_1276), dihydrolipoamide dehydrogenase (BCO26_1640), succinate dehydrogenase (or fumarate reductase) cytochrome b subunit, b558 family (BCO26_1913) and isocitrate dehydrogenase (BCO26_1963) were all significantly down-regulated under sodium lactate stress. Table 8. Significantly regulated genes involved in ‘citrate cycle (TCA cycle)’ under lactate stress. Gene ID Description FDR Fold change BCO26_0741 class II fumarate hydratase 4.90E-06 3.69 BCO26_0962 pyruvate dehydrogenase E1 component subunit alpha 7.36E-05 3.15 BCO26_0963 transketolase central region 1.52E-03 2.48 BCO26_0964 hypothetical protein BCO26_0964 6.32E-03 2.13 BCO26_0965 dihydrolipoamide dehydrogenase 1.36E-02 1.92 BCO26_1352 aconitate hydratase 1 8.97E-05 3.10 BCO26_1409 dihydrolipoamide dehydrogenase 9.60E-03 2.17 BCO26_1117 succinyl-CoA synthetase subunit beta 1.58E-05 -3.51 BCO26_1118 succinyl-CoA synthetase subunit beta 1.58E-05 -3.51 BCO26_1275 succinyl-CoA synthetase subunit alpha 3.80E-04 -2.83 BCO26_1276 2-oxoglutarate dehydrogenase, E1 subunit 4.90E-02 -1.54 BCO26_1640 2-oxoglutarate dehydrogenase, E2 subunit, dihydrolipoamide succinyltransferase 4.19E-02 -1.59 BCO26_1913 dihydrolipoamide dehydrogenase 3.82E-03 -2.27 BCO26_1963 succinate dehydrogenase (or fumarate reductase) cytochrome b subunit, b558 family 1.29E-02 -2.09 BCO26_0741[106]^* class II fumarate hydratase 2.00E-02 1.96 [107]Open in a new tab *Significantly regulated gene under calcium lactate stress, others are significantly regulated genes under sodium lactate stress Only one ‘citrate cycle (TCA cycle)’ gene was found to be significantly upregulated under calcium lactate stress condition (CA group) compared with non-stress condition (GY group) ([108]Table 8). This gene was class II fumarate hydratase (BCO26_0741), which was also expressed at higher level under sodium lactate stress compared to non-stress condition. Therefore, the TCA cycle of B. coagulans 2–6 was less affected by calcium lactate stress compared to sodium lactate stress. Identification of genes with at least five-fold up-regulation under lactate stress The genes with at least five-fold up-regulation under lactate stress were identified because they are potential targets for genetic engineering to further improve B. coagulans as a producer of L-lactic acid. Under sodium lactate stress condition (NA group), 28 genes were upregulated greater than fivefold ([109]S7 Table). Only three genes were upregulated by greater than fivefold under calcium lactate stress condition (CA group) ([110]S8 Table). Lambda repressor-like DNA-binding domain-containing protein (BCO26_0364) and heat shock protein Hsp20 (BCO26_1778) were both upregulated under the two types of lactate stress. The function of small heat shock proteins is to mediate correct protein folding and coping with stress-induced cell damage (Guzzo 2012). Two of the three small heat shock proteins Hsp20 in B. coagulans 2–6 were significantly upregualted under the two types of lactate stress. Major facilitator superfamily permease (BCO26_0524) was specifically upregulated greater than fivefold under calcium lactate stress ([111]S8 Table). Five genes were specifically upregulated greater than six-fold under sodium lactate stress, namely, hypothetical protein (BCO26_0739), ATPase AAA-2 domain-containing protein (BCO26_0934), hypothetical protein (BCO26_0740), copper ion binding protein (BCO26_2102) and hypothetical protein (BCO26_2103) ([112]S7 Table). Verification of RNA-seq data Reliability of the RNA-seq data was assessed by qRT-PCR analysis. In total, 5 genes were selected for validation. Although the magnitude of fold-change in gene expression varied between the two data sets, similar trends were observed in qRT-PCR and RNA-seq results, which supported the validity of the RNA-seq data ([113]S9 Table). Discussion B. coagulans strains have attracted significant attention for L-lactic acid production in the past few years [[114]5, [115]11, [116]25–[117]28]. Most of the research on this strain has focused on fermentation process. Qin et al. [[118]11] and Ye et al. [[119]28] reported that Ca(OH)[2] and CaCO[3] are more suitable neutralizing agents than NaOH for getting higher lactic acid titer. They presumed that it was attributable to lower product inhibition by calcium lactate compared to sodium lactate. However, little is known about the molecular mechanisms by which B. coagulans responds to different lactates. In this study, the gene expression patterns of B. coagulans 2–6 responding to sodium lactate and calcium lactate were investigated by RNA-seq and comparative transcriptome analysis. Gene ontology assignments and KEGG pathway enrichment analysis of the DEGs were performed. It was found that the influence of sodium lactate on gene expression in B. coagulans 2–6 was greater than that of calcium lactate. Sodium lactate induced significant upregulation of 543 genes and down-regulation of 670 genes, while calcium lactate induced significant upregulation of 398 genes and down-regulation of 314 genes. It was also found that lactate fermentation was less affected by calcium lactate stress than by sodium lactate stress. Key genes of the glycolysis pathway were most significantly upregulated under calcium lactate stress, while L-lactate dehydrogenase (–1.52-fold), pyruvate kinase (–3.98-fold) and 6-phosphofructokinase (–1.41-fold) were down-regulated under sodium lactate stress. Pyruvate dehydrogenase complex, which converts pyruvic acid to acetyl-CoA, was significantly upregulated (2.8-fold) under sodium lactate stress. These changes are expected to have significant influence on the formation of lactic acid from pyruvic acid. Taken together, our findings provide new mechanistic insights into higher lactic acid titer obtained using CaCO[3] as neutralizing agent compared to that obtained using NaOH [[120]5, [121]11]. Our results also provide potential targets for genetic engineering to further improve B. coagulans as an L-lactic acid producer. ABC transporters play important roles in response to lactate stress. A large number of the ABC transporters in B. coagulans 2–6 were found to be significantly down-regulated under the two types of lactate stress. We also identified ABC transporters that were significantly upregulated under the two types of lactates stress or specifically upregulated by one lactate. High expression of ABC transporter genes may be of benefit to B. coagulans 2–6 to maintain intercellular homeostasis; as such, these genes are good potential targets for genetic engineering. Furthermore, we identified genes with at least five-fold upregulation under lactate stress. Some of them were well-known stress-induced proteins, such as heat shock protein Hsp20 [[122]29] and major facilitator superfamily protein [[123]30], while others were previously unknown proteins, whose functions could be characterized in future studies. Wang et al. [[124]12] carried out a comparative proteomic analysis of B. coagulans in response to lactate stress. Using two-dimensional electrophoresis coupled with mass spectrometric identification, they identified four highly expressed proteins under calcium lactate stress. The same four proteins were also identified in this study, which indicates that our results are consistent with previous findings on this topic. In addition, we found that expressions of lactate dehydrogenase ([125]ACR02673), cysteine synthase A ([126]AEH52120) and ribosomal protein L7/L12 ([127]AEH52167) were upregulated under calcium lactate stress but down-regulated under sodium lactate stress. Gene expression of aldo/keto reductase ([128]AEH52329) was upregulated under calcium lactate stress and sodium lactate stress. Thus, our RNA-seq dataset provides a rich account of gene expression changes in B. coagulans in response to lactate stress, which could be highly useful for future efforts to improve lactic acid production by genetic engineering. Conclusions In summary, RNA-seq technology was employed in this study to characterise the transcriptome of B. coagulans under non-stress, sodium lactate stress and calcium lactate stress conditions. The responses of B. coagulans 2–6 to sodium lactate stress include: more genes changed expression than that in response to calcium lactate stress, significant regulation of ‘amino sugar and nucleotide sugar metabolism’ genes and negative influence on ‘glycolysis/gluconeogenesis’ and ‘citrate cycle (TCA cycle)’ genes. The responses of B. coagulans 2–6 to calcium lactate stress include: less genes changed expression than that in response to sodium lactate stress, significant regulation of ‘ABC transporters’ genes and positive influence on ‘glycolysis/gluconeogenesis’ genes. These results provide new insights into the responses of B. coagulans to different lactate stresses. Notably, our findings not only provide a robust database for investigating the functions of genes induced by lactate stress but also identify potential targets for genetic engineering to further improve L-lactic acid production by B. coagulans. Supporting Information S1 Table. The primers for qRT-PCR analysis with target gene information. (DOC) [129]Click here for additional data file.^ (27.5KB, doc) S2 Table. Significantly down-regulated genes involved in ‘ABC transporters’ under calcium lactate stress. (DOC) [130]Click here for additional data file.^ (44KB, doc) S3 Table. Significantly down-regulated genes involved in ‘ABC transporters’ under sodium lactate stress. (DOC) [131]Click here for additional data file.^ (53KB, doc) S4 Table. Significantly down-regulated genes involved in ‘amino sugar and nucleotide sugar metabolism’ under sodium lactate stress. (DOC) [132]Click here for additional data file.^ (31.5KB, doc) S5 Table. Significantly down-regulated genes involved in ‘glycolysis/gluconeogenesis’ under sodium lactate stress. (DOC) [133]Click here for additional data file.^ (28.5KB, doc) S6 Table. Significantly down-regulated genes involved in ‘glycolysis/gluconeogenesis’ under calcium lactate stress. (DOC) [134]Click here for additional data file.^ (26KB, doc) S7 Table. Genes with at least fivefold upregulation under sodium lactate stress. (DOC) [135]Click here for additional data file.^ (45.5KB, doc) S8 Table. Genes with at least fivefold upregulation under calcium lactate stress. (DOC) [136]Click here for additional data file.^ (25KB, doc) S9 Table. The expression information of selected genes for qRT-PCR under lactate stress condition. (DOC) [137]Click here for additional data file.^ (29.5KB, doc) S1 Data. GO assignments of differentially expressed genes (CA vs GY). (XLS) [138]Click here for additional data file.^ (120.5KB, xls) S2 Data. GO assignments of differentially expressed genes (NA vs GY). (XLS) [139]Click here for additional data file.^ (205.5KB, xls) Data Availability Relevant data are within the paper and its Supporting Information files. All sequencing data have been deposited in the NCBI Sequence Read Archive under accession numbers SRX700697 (GY), SRX700698 (CA) and SRX 700710 (NA). Funding Statement This work was supported by the National Natural Science Foundation of China (31200078), the Shandong Provincial Natural Science Foundation, China (ZR2001CQ025), and the CAS Key Laboratory of Microbial Physiological and Metabolic Engineering, Institute of Microbiology, Chinese Academy of Sciences (KLIM-201302). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. References