Lactococcus lactis subsp. cremoris NCDO712 is an industrial strain containing six plasmids of which one, pLP712, carries lactose and casein utilization genes (Tarazanova et al., 2016). L. lactis MG1363 is its plasmid-free derivative (Gasson, 1983). L. lactis NCDO712 was grown as a standing culture at 30°C in M17 broth (Difco, Becton Dickinson, Le Pont de Claix, France) containing 0.5% (w/v) glucose (GM17), and on GM17 agar plates.

All stress conditions were applied in two biological replicates by starting with two single colonies of L. lactis NCDO712 grown on GM17 agar plates. These were each used to inoculate 10 ml of fresh GM17 media and grown overnight. Each overnight culture was diluted 100-fold in a bottle with 500 ml GM17 and grown until an optical density at 600 nm (OD₆₀₀) of 0.9 was reached. The content of each bottle was divided over seven 50-ml Greiner tubes and centrifuged at 4,000 g for 1.5 min. Subsequently, the cultures from each bottle were subjected to all conditions tested. The cell pellets were re-suspended in fresh GM17 containing 0.25% glucose (G^*M17) and the following properties: control (G^*M17 at 30°C), cold (G^*M17 at 10°C), heat (G^*M17 at 42°C), acid (G^*M17 set at pH 4.5 with lactic acid), osmotic stress (G^*M17 containing 2.5% NaCl), and oxidative stress (G^*M17, shaking at 250 rpm). The cold stress was applied in an incubator while the heat stress was performed in a water bath to maintain the cold and heat conditions stable. For starvation stress, the cell pellets were re-suspended in filter-sterilized phosphate-buffered saline (PBS). The stress conditions were applied for 5 min, after which the cells were harvested by centrifugation at 10,000 rpm for 1 min, snap frozen in liquid nitrogen and stored at −80°C prior to RNA isolation.

RNA was isolated as described previously (van der Meulen et al., 2016). All procedures were executed at 4°C unless otherwise stated and all solutions were DEPC-treated and subsequently autoclaved. Frozen cell pellets were re-suspended in 400 μl TE-buffer (10 mM Tris-HCl, 1 mM EDTA, pH 7.4) and added to 50 μl 10% sodium dodecyl sulfate (SDS), 500 μl phenol/chloroform (1:1 v/v), and 0.5 g glass beads (75–150 μm, Thermo Fischer Scientific, Rockford, IL, United States). The cells were disrupted by shaking 2 times for 45 s in a Biospec Mini-BeadBeater (Biospec Products, Bartlesville, OK, United States) with cooling on ice for 1 min between the shaking steps. Subsequently, the cell suspension was centrifuged at 14,000 rpm for 10 min. The upper phase containing the nucleic acids was treated with 500 μl chloroform and centrifuged as above. Nucleic acids in the water phase were precipitated by sodium acetate and ethanol. The nucleic acid pellet was re-suspended in 100 μl buffer consisting of 82 μl MQ, 10 μl 10x DNase I buffer, 5 μl RNase-free DNase I (Roche Diagnostics GmbH, Mannheim, Germany), and 3 μl RiboLock RNase inhibitor (Fermentas/Thermo Scientific, Vilnius, Lithuania), and treated for 30 min at 37°C. The RNA was then purified using standard phenol/chloroform extraction and sodium acetate/ethanol precipitation. RNA pellets were re-suspended in 50 μl elution buffer from the High Pure RNA Isolation Kit (Roche Diagnostics, Almere, the Netherlands) and stored at −80°C.

RNA concentration was measured with a Nanodrop ND-1000 (Thermo Fischer Scientific). As a measure of RNA quality, the integrity of the 16S/23S rRNA and the presence of any DNA contamination were assessed by using an Agilent 2100 Bioanalyser (Agilent Technologies, Waldbronn, Germany). cDNA library was prepared by employing a ScriptSeq^TM Complete Kit for Bacteria (Epicentre, Madison, WI, United States) including Ribo-Zero™ for rRNA removal. The cDNA libraries were sequenced at Otogenetics Corporation (Norcross, GA, United States) on an Illumina HiSeq2000.

Raw sequence reads were analyzed for quality and trimmed with a PHRED score >28. Read alignment was performed on the genomic DNA of L. lactic NCDO712 (nucleotide sequences of the chromosome and all six plasmids of NCDO712; Tarazanova et al., 2016) using Bowtie 2 (Langmead and Salzberg, 2012). RKPM values were used as an input for the T-REx analysis pipeline (de Jong et al., 2015) together with a text file describing the factors, contrasts, and classes. In the class file, genes from the NCDO712 plasmids were colored green while sRNAs were colored red (see Supplementary Information S1). T-REx was used to perform all statistical analyses (de Jong et al., 2015).

The RNA-seq data is publically available in GEO under accession number GSE98499.

The transcriptomic response of L. lactis subsp. cremoris NCDO712 (hereafter named L. lactis NCDO712) after 5 min of cold, heat, acid, osmotic, oxidative, or starvation stress was determined by high-throughput RNA sequencing. This resulted in a total of 246M of reads of which 209M reads (85%) were successfully mapped on the genome and plasmids of L. lactis NCDO712. The libraries varied between 11M and 19M reads per individual sample (Figure 1A). The data was normalized using the T-REx software (de Jong et al., 2015) and plotted in a box plot of normalized signals for all samples (Figure 1B). A Principle Component Analysis (PCA) shows that the stress conditions are statistically well-distributed from each other (Figure 1C). The transcriptome of cells exposed to osmotic stress or to starvation were most different from that of the control. The absolute numbers of differentially expressed genes underpin these observations; the largest transcriptome changes were observed after starvation (756 genes involved), while oxidative stress had the least impact (91 affected genes). Table 1 gives an overview of the counts of the differentially expressed genes and Figure 2 shows the distribution of affected genes for each stress condition. To gain insight in the distribution of the stress-responsive genes, those of which the expression changed highly (fold change ≥ 5 and p ≤ 0.01) under all conditions were visualized in a heatmap. T-REx was used to pinpoint nine clusters that vary strongly in size and in the functions of the constituting genes. One cluster is rich in stress-related genes, while another one contains predominantly sRNAs. See Figure 3 for the heatmap and the complete list of cluster descriptions. Genes that were affected by a fold change ≥ 10-fold and p ≤ 0.01) are listed in Table 2 and are discussed in more detail below.

We assessed the expression of the 186 sRNAs that have recently been identified in the genome of L. lactis MG1363 (van der Meulen et al., 2016). The majority of the sRNAs that were significantly changed after applying acid, oxidative, or starvation stress showed a decrease in expression. In contrast, after cold stress more sRNAs were upregulated instead of downregulated (see Figure 4). Of the 186 sRNAs, the expression of 110 was significantly changed after applying at least one stressor, while 42 sRNAs were differentially expressed under at least three stress conditions (Table 3). This list of sRNAs was restricted to those with a logCPM value >1. The expression of some of these 110 sRNA genes was highly affected under only one specific stress condition. For example, LLMGnc_152 (10.7-fold) and LLMGnc_153 (6.0-fold) were increased specifically after cold stress, while LLMGnc_176 (6.3-fold) showed a higher expression after salt stress. The sRNA LLMGnc_092 (−9.2-fold) was only decreased after starvation, while LLMGnc_025/064/065/073 were downregulated in all conditions.

Transfer RNAs (tRNAs) play a crucial role in the translation of mRNAs and it is important for cells to balance their tRNA levels, as well as to ensure optimal utilization of amino acids. Most of the L. lactis tRNAs are downregulated under all of the stresses applied. A number of tRNAs, such as NCDO_2402 (Val-CAG) and NCDO_2022 (Lys-TTT), are upregulated under some of the conditions. A distinct tRNA response is observed upon cold stress; seven tRNAs are upregulated by at least 2-fold. Exposure to acid or starvation stress induced the most severe changes in tRNA expression. See Figure 4 for a complete overview of tRNA expression under the various stress conditions.

During a 5-min cold stress at 10°C, 412 genes were differentially expressed (Table 1). The most highly upregulated transcripts include those from the zit operon, a gene cluster involved in Zn²⁺ uptake and regulation (Llull and Poquet, 2004). Also, expression of the gene sugE, encoding a presumed multidrug resistance protein, was increased (~18-fold), as well as genes involved in nucleotide synthesis (pur and pyr operons) and the fab and acc operons for saturated fatty acid biosynthesis. Transcripts encoding Cold Shock Proteins A, B, C, and D were upregulated at least 4-fold, as would be expected in cells under cold stress (Wouters et al., 1998). The gene of an uncharacterized protein (Llmg_1848) with high sequence similarity to a bacteriocin in other L. lactis species, is located downstream of cspA and was also upregulated by a factor of ~10. Downregulation was observed e.g., for the fructose utilization fru operon (~18-fold), the lysine-specific permease lysP gene (~8-fold) and the ribosomal RNA 5S (~9-fold). Thirty-nine sRNAs from intergenic regions were significantly affected, of which 12 were changed at least 5-fold (see Figure 4 and Table 3).

At high temperatures, proteins may be at risk of denaturation and cells may encounter difficulties in the synthesis of new proteins (Nguyen et al., 1989; Parsell and Lindquist, 1993). The genes of several protein chaperones, such as GroEL, GroES, DnaJ, DnaK, and GrpE are usually upregulated after heat stress. The chaperones aid in protein folding and maturation. Their genes were upregulated ~14–18 times (except dnaJ, which was increased ~5-fold) after L. lactis was placed for 5 min at 42°C (see Table 2). The Clp protease genes clpE, clpP, and clpB were also upregulated, as was the gene with an unknown function upstream of clpE, llmg_0527. Surprisingly, expression of genes belonging to the arginine deaminase pathway (arc operon) were increased ~11-fold, while this operon was downregulated after 30 min of incubation of L. lactis IL1403 at 42°C in previous work (Xie et al., 2004). An operon predicted to be involved in the utilization of maltose (llmg_0485-llmg_0490), the ribose operon (rbsABCDK) and llmg_1210, predicted to encode the drug resistance transporter EmrB, were all upregulated. We also observed that the lacR-lacABCDEF gene cluster, located on the largest plasmid pLP712 of L. lactis NCDO712 (Wegmann et al., 2012) and involved in lactose utilization, was upregulated ~3-fold. The heat treatment caused a decrease by a factor of 2–3 of several rRNA species, as well as of most tRNA-transcripts (Figure 5). An operon (llmg_2513-llmg_2515) containing a gene for the universal stress protein A (UspA) and two uncharacterized genes was decreased most severely after heat shock, followed in severity of fold change by the multidrug resistance protein B gene (lmrB, llmg_1104). In total, 37 sRNAs were differentially expressed (Table 3). Expression of two of them was decreased over 5-fold: LLMGnc_065 (−6.0) and LLMGnc_079 (−6.2).

Genes for de novo synthesis of pyrimidines (pyr) and purines (pur) were highly upregulated after 5 min of exposure to pH 4.5. Expression of these operons was also seen during starvation and, to a lesser extent, upon cold shock. Among the top upregulated genes were metCcysK encoding a cystathionine gamma-lyase and cysteine synthase, respectively, involved in cysteine and methionine metabolism. The genes llmg_0333-0340 were all affected by the acid stress applied, although they are not transcribed from the same operon. Among these genes are those of a putative methionine ABC transporter system (llmg_0335-0340 upregulated), and thiT (llmg_0334) encoding the thiamine transporter (Erkens and Slotboom, 2010) and a putative transcriptional regulator gene (llmg_0333) (both downregulated). The fab operon, responsible for the biosynthesis of saturated fatty acids (SFAs) for membrane phospholipids, was upregulated. Expression of the operon llmg_2513-2515 was increased, while these genes were strongly downregulated after heat (see above) and salt stress (see further). The putative Mn^2+/Fe²⁺ transporter gene mntH, which is located on plasmid pNZ712, was upregulated 2.5-fold under pH 4.5 stress. Interestingly, the chromosomal gene mgtA, putatively specifying Mg²⁺ transport, was downregulated 3.7-fold. Upregulation was furthermore observed for llmg_1915-1917 (ykgEFG) and for the following genes with predicted functions; llmg_1066 (unknown), llmg_1133 (exonuclease), and llmg_1702 (glutathione reductase). The fruAKR operon for fructose transport, conversion of fructose 1-phosphate to fructose 1,6-bisphosphate and their regulation by the FruR repressor (Barriere et al., 2005), was downregulated on average ~28-fold. The fhuABCD operon for putative ferric siderophore transport was downregulated by a factor 3–7. Of the 51 sRNA genes of which the expression was significantly changed, only 8 were upregulated; among these was LLMGnc_121, which was exclusively upregulated after acid stress (Table 3).

Osmotic stress was induced by adding 2.5% NaCl to the cell culture for 5 min. Among the highest upregulated genes are those of an operon (llmg_2163 - llmg_2164) specifying a putative stress-responsive transcriptional regulator with a PspC domain (Llmg_2163). Both genes are ~10-fold upregulated; it was also induced upon overproduction in L. lactis of the membrane protein BcaP (Pinto et al., 2011) and after exposure to the bacteriocin Lcn972 (Martínez et al., 2007). A deletion mutant of llmg_2164 was shown to be very sensitive to NaCl (Roces et al., 2009). In E. coli, the psp operon is induced after application of various types of stresses including salt stress (Brissette et al., 1991). As expected, induction was seen of the genes specifying the glycine betaine ABC transport system BusAA-BusAB (~7- to 8-fold). Some of the responses observed during the exposure to the high concentration of salt were similar to those seen after heat stress. In particular, transcripts encoding the chaperones GroEL, GroES, DnaJ, DnaK, and GrpE were upregulated in the same fold change range. Induction of these proteins has been reported previously for both heat and salt stress (Kilstrup et al., 1997). Both operons for oligopeptide transport were also upregulated under the salt stress applied here.

Genes encoding transporters for various substrates were strongly downregulated, among which the PTS IIA component genes fruA and ptcA, that of the IIB PTS component, bglP, lmrB specifying multidrug resistance protein B (both llmg_0967 and llmg_1104), msmK, encoding a multiple sugar ABC transporter and amtB involved in ammonium transport. Also, the gene encoding the glycerol uptake facilitator protein GlpF2 was downregulated by ~10-fold, as was the lac operon on pLP712. As shown before (Xie et al., 2004), the potABCD operon involved in spermidine/putrescine transport and the fatty acid biosynthesis operons fab and acc were downregulated under osmotic stress. As mentioned above, the llmg_2513-2515 operon was downregulated. From the 44 significantly affected sRNA genes, the expression of LLMGnc_036 (214-fold down), LLMGnc_118 (75-fold down), and LLMGnc_176 (6.3-fold up) was specifically only changed after salt stress using a threshold of 5-fold change.

From all stress conditions tested, oxidative stress applied by shaking of the culture resulted in the least number of differentially expressed genes (Table 1). The ones that did change did so with relatively minor fold changes. Among the few upregulated genes were those of the pathway for saturated fatty acid biosynthesis, including fabT, the transcriptional repressor of this route (Eckhardt et al., 2013). Downregulation was seen of arcABD1C1C2, llmg_1915-1917 (ykgEFG), and of genes involved in the uptake and/or utilization of maltose, trehalose, and lactose. The heat shock chaperone genes groEL and groES were downregulated ~2-fold. Unexpectedly, the gene for the manganese superoxide dismutase SodA, an enzyme well-known for its role in oxidative stress, was decreased 2-fold relative to the unstressed control. Twenty sRNAs were affected by at least 2-fold. LLMGnc_019 was downregulated 9.4-fold, while it was upregulated in all the other stress conditions (see Table 3).

Incubation of the cells for 5 min in PBS greatly affected the transcriptome of L. lactis, as witnessed by the large number of 150 genes of which the expression had changed significantly, and at least by 5-fold. Different biological functions were switched on or off in response to the sudden absence of basically all nutrients. For example, the operons for the de novo biosynthesis of nucleotides, pyr, pur, and car, were all highly upregulated. Moreover, the genes of several transport systems were upregulated in an apparent attempt to import a (any) carbon source (lactose: lac operon, multiple sugar ABC transporter: msmK, ribose: rbs operon, cellobiose/glucose: ptcC), ions (manganese: mtsA, zinc: zit operon), amino acids (arginine: arc, methionine: metQ, branched-chain amino acid transporter: ctrA/bcaP), and vitamins (riboflavin: rib operon). On the other hand, some uptake system genes were downregulated, such as thiT, specifying the thiamin transporter, and the fru operon for fructose uptake and utilization. Strong downregulation was also observed for the genes of all six cold shock proteins. Dozens of (predicted) transcriptional regulators were affected upon the starvation stress applied here, among which all of the spxA genes. Interestingly, the expression profiles were very different. For instance, spxA with locus tag llmg_0640 was decreased ~19-fold, while expression of llmg_1703 and llmg_1130 was increased 12.5- and 53-fold, respectively. Notably, the gene for the putative transcriptional repressor CadC, which is located on plasmid pSH73 (Tarazanova et al., 2016), was 4.9-fold upregulated. Besides a strong decrease in the expression of different tRNA genes, also transcripts for ribosomal proteins were downregulated. Starvation changed the expression of 75 sRNAs, of which the majority was downregulated. The nine sRNA genes that were upregulated after starvation were also increased under at least one of the other conditions tested here, with the exception of LLMGnc_060, which was only affected after starvation.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.