Comparative Transcriptome Analysis Reveals Different Molecular Mechanisms of Bacillus coagulans 2-6 Response to Sodium Lactate and Calcium Lactate during Lactic Acid Production

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 [1][2][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 [2,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 [5]. Previous studies indicated that the maximum yield of lactic acid can exceed 200 g/L when calcium lactate accumulates in the fermentation broth [6,7]; however, the lactic acid yield is generally less than 125 g/L when sodium lactate is formed [5,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 [9,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 [5,11]. Under the same concentration of lactates, sodium lactate results in stronger stress on B. coagulans 2-6 compared with calcium lactate [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. [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.

Bacterial strains and cultivation conditions
The B. coagulans wild-type strain 2-6 [11,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 [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: NC_015634) using the short sequence alignment software Bowtie 2 [14]. To identify differentially expressed genes (DEGs) between the two different treatment conditions, the expression level for each transcript was calculated by RSEM (http:// deweylab.biostat.wisc.edu/rsem/) using the fragments per kilobase of exon per million mapped reads method [15]. The R statistical package software EdgeR (http://www.bioconductor.org/ packages/2.12/bioc/html/edgeR.html) was employed to quantify differential gene expression [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 [17,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 [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 S1 Table. The relative gene expression data were analyzed using the 2 -ΔΔCt method as described by Tan et al. [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.

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).
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 Fig 2. 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 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 Table 1, S1 Data and 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.  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 [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 [21] and Bacillus subtilis [22].
KEGG pathway enrichment analysis revealed that 'ABC transporters' was significantly affected by calcium lactate stress. A total of 10 significantly upregulated genes ( Table 2) and 27 significantly down-regulated genes (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 (Table 3) and 35 significantly down-regulated genes (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 (S2 Table and 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).
Under the two types of lactate stress conditions, two significantly upregulated genes involved in 'glycolysis/gluconeogenesis' were found, namely, alcohol dehydrogenase GroES domaincontaining 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 [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) ( 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), 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) ( 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 Llactic acid. Under sodium lactate stress condition (NA group), 28 genes were upregulated greater than fivefold (S7 Table).
Only three genes were upregulated by greater than fivefold under calcium lactate stress condition (CA group) (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 (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) (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 (S9 Table).

Discussion
B. coagulans strains have attracted significant attention for L-lactic acid production in the past few years [5,11,[25][26][27][28]. Most of the research on this strain has focused on fermentation process. Qin et al. [11] and Ye et al. [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 Llactate 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 [5,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 [29] and major facilitator superfamily protein [30], while others were previously unknown proteins, whose functions could be characterized in future studies.
Wang et al. [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 (ACR02673), cysteine synthase A (AEH52120) and ribosomal protein L7/L12 (AEH52167) were upregulated under calcium lactate stress but down-regulated under sodium lactate stress. Gene expression of aldo/keto reductase (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.