ABSTRACT
Psychrotrophic lactic acid bacteria (LAB) are the prevailing spoilage organisms in packaged cold-stored meat products. Species composition and metabolic activities of such LAB spoilage communities are determined by the nature of the meat product, storage conditions, and interspecies interactions. Our knowledge of system level responses of LAB during such interactions is very limited. To expand it, we studied interactions between three common psychrotrophic spoilage LAB (Leuconostoc gelidum, Lactococcus piscium, and Lactobacillus oligofermentans) by comparing their time course transcriptome profiles obtained during their growth in individual, pairwise, and triple cultures. The study revealed how these LAB employed different strategies to cope with the consequences of interspecies competition. The fastest-growing bacterium, Le. gelidum, attempted to enhance its nutrient-scavenging and growth capabilities in the presence of other LAB through upregulation of carbohydrate catabolic pathways, pyruvate fermentation enzymes, and ribosomal proteins, whereas the slower-growing Lc. piscium and Lb. oligofermentans downregulated these functions. These findings may explain the competitive success and predominance of Le. gelidum in a variety of spoiled foods. Peculiarly, interspecies interactions induced overexpression of prophage genes and restriction modification systems (mechanisms of DNA exchange and protection against it) in Lc. piscium and Lb. oligofermentans but not in Le. gelidum. Cocultivation induced also overexpression of the numerous putative adhesins in Lb. oligofermentans. These adhesins might contribute to the survival of this slowly growing bacterium in actively growing meat spoilage communities.
IMPORTANCE Despite the apparent relevance of LAB for biotechnology and human health, interactions between members of LAB communities are not well known. Knowledge of such interactions is crucial for understanding how these communities function and, consequently, whether there is any possibility to develop new strategies to interfere with their growth and to postpone spoilage of packaged and refrigerated foods. With the help of controlled experiments, detailed regulation events can be observed. This study gives an insight into the system level interactions and the different competition-induced survival strategies related to enhanced uptake and catabolism of carbon sources, overexpression of adhesins and putative bacteriocins, and the induction of exchange of genetic material. Even though this experiment dealt with only three LAB strains in vitro, these findings agreed well with the relative abundance patterns typically reported for these species in natural food microbial communities.
INTRODUCTION
Psychrotrophic lactic acid bacteria (LAB) are prevailing spoilage organisms in a variety of refrigerated modified-atmosphere-packaged (MAP) meat products (1–3). Different spoilage LAB species form dense communities on meat surfaces and are assumed to interact with each other (4). These interactions eventually affect the population structure and global activity of a microbial community. Generally, interactions between LAB have been shown to result in various outcomes, from growth promotion and increased survival due to symbiotic relationships (5–7) to growth inhibition due to competition (7–10). The individual mechanisms behind bacterial interspecies interactions include increasing nutrient availability and creating more-favorable growth conditions during symbiotic relationships (11, 12) or developing attack (release of harmful molecules, such as bacteriocins or reactive oxygen species [ROS], and contact-induced killing [13–15]) defense mechanisms and increasing the efficiency of nutrient uptake during competitive relationships (16). Despite the apparent relevance of LAB for biotechnology and human health, knowledge of their system level interactions is limited. Only few studies have investigated transcriptional responses of food fermentation, probiotic, or pathogenic LAB (6, 17–20) during their cocultivation with other LAB and non-LAB species. To our knowledge, no information is available on transcriptome responses during interactions between food spoilage-associated LAB constituting communities in MAP foods.
The most abundant LAB genera associated with meat spoilage are Leuconostoc, Carnobacterium, Lactobacillus, Lactococcus, and Weissella (2). The relative abundance of different species depends on the type of meat and preservation technology (21). However, in the majority of packaged and refrigerated meat products, as well as in many other foods, Leuconostoc gelidum has been found to be the most predominant species, followed by Lactococcus piscium, which also belongs to predominant microbiota at the late shelf life stage (2, 3, 22). The documented spoilage activities for these two bacteria include excessive formation of buttery and sour off-odors (23, 24). In addition, Le. gelidum spoilage has been associated with discoloration and the production of gas and slime (23). For Le. gelidum, rapid proliferation and outgrowth from a low level of contamination have been detected during storage of MAP foods (25). This, together with the ability of some strains to adhere strongly to food contact surfaces (26) and the wide metabolic capabilities of this species (23), might explain Le. gelidum's competitive success and prevalence.
Alongside prevailing bacteria, microbial meat spoilage communities contain less abundant species that are not completely replaced by the predominant species and are able to coexist with them. One such species isolated mainly from meat sources (27) is Lactobacillus oligofermentans. Unlike the majority of LAB, Lb. oligofermentans strains did not originally grow well on glucose; instead, they efficiently utilized pentoses (27), the concentration of which in their free form is very low in comparison to glucose in meat (28, 29). Nevertheless, Lb. oligofermentans is able to maintain quite high levels (10 to 18% of spoilage LAB) in some actively growing cocultures in meat (27).
To obtain insights into system level interactions between food spoilage LAB, we chose to study the transcriptome responses of three spoilage LAB species during their cocultivation pairwise and in triple culture in vitro in a medium containing glucose as the main carbon source. Le. gelidum subsp. gasicomitatum LMG 18811T (23, 30), Lc. piscium MKFS47 (31, 32), and Lb. oligofermentans LMG 22743T (27, 33) originating from MAP broiler meat were selected for the study. These species have been detected together in various food spoilage communities, and we aimed to compare their responses in the presence of other species and under monoculture conditions.
RESULTS
Growth of the three LAB in individual and mixed cultures.The relative growth dynamics of the three LAB grown individually differed markedly depending on the type of growth measurement used, i.e., optical density values at 600 nm (OD600) or CFU per milliliter. According to the OD600-based growth curves (assumingly correlated with dry biomass growth [34]) (Fig. 1A, B, C, and D), Lc. piscium culture had the highest accumulated biomass concentration until 9 h, after which the biomass concentrations of Lc. piscium and Le. gelidum became similar. Based on the same type of measurement, the culture of Lb. oligofermentans had notably lower biomass than the other two species. In terms of viable-cell concentration (CFU/milliliter) (Fig. 1E, F, G, and H), the growth rate of Le. gelidum was higher than that of Lc. piscium and Lb. oligofermentans, the last two species showing similar growth rates. The observed difference in relative growth rates of Lc. piscium obtained using different measurements might indicate that part of its cell population was nonculturable (dead or dormant) or formed chains/clumps. Alternatively, the correlation between the sample optical density and overall cell concentration for Lc. piscium might differ significantly from that in the other two LAB. For all the mixed cultures, their total OD600 values and viable cell counts did not exceed the values obtained for the individual cultures of the coculture members (Fig. 1). Concentrations of viable cells of each member in the cocultures were obtained using selective media (Fig. 2). Generally, for pairwise cocultures (Fig. 2A, B, and C), abundance relationships between the species constituting the cocultures were similar to that between their individual cultures (Fig. 1E, F, G, and H). Le. gelidum seemed to dominate in both of its pairwise cocultures (based on the retention of the differences in means at the late time points, though not to the level of statistical significance). Nevertheless, coculturing resulted in a plateauing of Le. gelidum concentration sooner (Fig. 2A, C, and D) than in its individual culture (Fig. 1E and H), which suggests that the two other species have a restrictive effect on the growth of Le. gelidum.
Growth curves based on the optical density (OD600) values (A to D) and viable cell counts (CFU/milliliter) (E to H) of the individual and mixed cultures of the three LAB. G, O, P, individual cultures of Le. gelidum, Lb. oligofermentans, and Lc. piscium, respectively; GP, OP, GO, GOP, cocultures (e.g., GP is a coculture of Le. gelidum and Lc. piscium). Panels A and E represent cultures obtained in the first batch, while other panels represent cultures obtained in the second batch (see Materials and Methods). The values are the averages for three biological replicates ± standard errors of the means (SEM). Statistically significant differences between individual cultures of each bacterium are marked: *, P value ≤ 0.05, two-tailed t test.
Growth curves of Le. gelidum (G), Lb. oligofermentans (O), and Lc. piscium (P) in their cocultures (GP, OP, GO, and GOP) based on viable cell counts (CFU/milliliter). G_in_GP indicates cell counts of Le. gelidum in the coculture of Le. gelidum with Lc. piscium. The values are the averages for three biological replicates ± SEM. Statistically significant differences (P value ≤ 0.05, two-tailed t test) between bacterial species in their pairwise cocultures are marked with black asterisks. In the triple culture, red and purple asterisks indicate statistically significant differences of cell counts of Le. gelidum and Lc. piscium, respectively, from the two other bacterial species.
General characteristics of the cocultivation-induced transcriptome responses.Transcriptome profiles of the LAB grown in individual cultures were compared with their corresponding profiles obtained during growth in their mixed cultures (two pairwise and one triple for each species) independently for each LAB species and time point (3, 5, and 11 h). The numbers of differentially expressed (DE) genes identified for each bacterium (Tables 1, 2, and 3) differed greatly depending on the time point and type of the coculture and constituted maximally 12%, 16%, and 16% of the genomes of Le. gelidum, Lc. piscium, and Lb. oligofermentans, respectively. For each species, we calculated the intersections between DE genes in different cocultures at each time point independently (see the Venn diagrams in Fig. S1 in the supplemental material). The number and proportion of DE genes that were common for the interactions with two different cohabiters in the pairwise interactions and had the same direction of changes (see Table S1 in the supplemental material) were much higher for Lb. oligofermentans (19 to 34%) than for Le. gelidum and Lc. piscium at all time points. This might be caused either by the higher consistency of Lb. oligofermentans data or by its more universal response to the coculturing with two other LAB than such responses from Le. gelidum and Lc. piscium, which might be more coculture specific. However, the number and proportion of such “overlapped” DE genes in Le. gelidum and Lc. piscium increased with time to reach a maximum at 11 h (5 to 7%).
Numbers and enriched COG categories of cocultivation-induced DE genes of Le. geliduma
Numbers and enriched COG categories of cocultivation-induced DE genes of Lc. pisciuma
Numbers and enriched COG categories of cocultivation-induced DE genes of Lb. oligofermentansa
To validate the differential expression results obtained using RNA sequencing (RNA-seq) data, the expression levels of 13 genes that spanned a range of expression fold changes (FC) induced in different cocultures were measured using droplet digital reverse transcription-PCR (dd-RT-PCR). Expression levels of the housekeeping genes infB and ftsQ as well as 16S rRNA were used to normalize the gene expression levels between samples. Overall, the dd-RT-PCR results were very consistent with the RNA-seq results (see Fig. S2 in the supplemental material). Moreover, the relative gene expression changes obtained using RNA-seq and dd-RT-PCR correlated linearly for the tested genes (Pearson's correlation coefficient, 0.984) (Fig. 3).
Comparison of relative expression changes (log2 FC) for 13 genes obtained using droplet digital reverse transcription-PCR (dd-RT-PCR) versus RNA sequencing (RNA-seq). A letter in parentheses after a gene name represents the source organism (G, P, or O represents Le. gelidum, Lc. piscium, or Lb. oligofermentans, respectively). Expression fold changes (FC) represent the ratio of gene expression (averaged across three replicates) between mixed and individual cultures, taken at the same time point. To test gene expression with dd-RT-PCR, the following mixed cultures were used: GP-11 h for nagB(G), glmS(G), pyrP(G), pfl(P), fruK(P), ilvA(P); GOP-11 h for ccpA(G), pdhD(G), mapA(G); OP-3 h for LACPI_0417(P), OP-11 h for ccpA(O), infC(O); GO-11 h for LACOL_1683(O), where, e.g., GP-11 h represents the mixed culture of Le. gelidum and Lc. piscium at 11 h. For dd-RT-PCR, the concentration of the housekeeping gene infB was used to normalize the gene concentrations between samples. For data points with a square marker, the changes in gene expression were statistically significant based on both RNA-seq (DESeq2 adjusted P value < 0.05) and dd-RT-PCR (two-tailed t test P value < 0.05). A full functional description of the genes can be found in Data Set S4.
To obtain a global overview of functions affected by the cocultivations, the lists of DE genes were tested for the enrichment of the clusters of orthologous groups (COG) functional categories (Tables 1 to 3 and S1), gene ontology (GO) terms, and metabolic pathways (see Data Sets S1, S2, and S3 in the supplemental material). The wide array of functions was affected by the interspecies interactions, and the enriched categories were more coculture specific for Le. gelidum and Lc. piscium than for Lb. oligofermentans. Cocultivation of LAB affected mostly the genes in the “Carbohydrate transport and metabolism” (“G”), “Posttranslational modification, protein turnover, chaperones” (“O”), “Translation, ribosomal structure and biogenesis” (“J”), and “Energy production and conversion” (“C”) COG categories. A selection of the enriched functional groups and pathways (Fig. 4) are described in more detail in the following sections.
Differential expression of genes, pathways, and functional groups in each of the three LAB during their cocultivation in the pairwise and triple cultures. For Le. gelidum, the columns “+ P,” “+ O,” and “+ P + O” represent differential expression of its genes/pathways during growth with Lc. piscium, Lb. oligofermentans, and both (Lc. piscium and Lb. oligofermentans), respectively. Upward red arrows indicate upregulation, and downward blue arrows indicate downregulation. The time points (3 h, 5 h, and 11 h) for the differential expression events are marked next to the arrows. Superscript a indicates catabolism and transport of hexoses and di-/oligosaccharides; superscript b indicates that cationic peptides might play the role of bacteriocins; superscript c indicates complete RM systems or their restriction enzymes.
Carbohydrate catabolism/transport, fermentation, and energy metabolism.In Le. gelidum, upregulated genes in the triple culture at 11 h were enriched in GO terms associated with carbohydrate catabolism/transport, pyruvate metabolism, and energy generation (Fig. 5A and Data Set S1). This was also reflected in the enrichment of the “G” COG category (Table 1). Particularly, many genes involved in carbohydrate-specific catabolism/transport (those encoding mannose, galactose, glucosamine, sucrose, trehalose, lactose, maltose, oligo-1,6-glucosides, and beta-glucosides), the phosphoketolase pathway (gnd2, xfp, gap, pgk, gpmA2, eno, and pyk), and carbon catabolite control (CCC) (ccpA, ptsH, and hprK) and two genes encoding pyruvate-dissipating enzymes (PDEs), pdhABCD and poxB, were upregulated (some of the genes are illustrated in Fig. 3). Notably, this effect was not as pronounced in the pairwise cocultures. Nevertheless, some genes involved in the phosphoketolase pathway were also overexpressed in the pairwise cultures with Lc. piscium (gap, pgk, gpmA2, eno) or Lb. oligofermentans (gap, gpmA2, eno, pyk), and acetate kinase genes ackA1 and ackA2 were upregulated only in the coculture with Lb. oligofermentans at 11 h. Peculiarly, during growth with Lc. piscium, two genes, encoding glutamine-fructose-6-phosphate aminotransferase (glmS) and glucosamine-6-phosphate deaminase (nagB), that catalyze the biosynthesis of glucosamine-6-phosphate for cell wall formation and catabolic degradation of glucosamine-6-phosphate, respectively, were strongly downregulated and upregulated, respectively, during the entire experiment (Fig. 3). The interaction with Lc. piscium clearly shifted the balance in the glucosamine metabolism of Le. gelidum from cell wall formation toward glucosamine catabolic degradation.
TreeMaps of the GO terms enriched in the DE genes during interspecies interactions of Le. gelidum, Lc. piscium, and Lb. oligofermentans. Only the most representative interactions are shown. The size of the box for each GO term is proportional to the absolute value of its log10 P value. Semantically similar GO terms form clusters distinguished by color. TreeMaps were built using the REVIGO Web server (57).
In Lc. piscium, DE genes during growth with Le. gelidum and in the triple culture at 11 h were enriched in the “G” COG category (Table 2). Genes that belonged to this category consisted mainly of downregulated carbohydrate catabolism/transport genes (for maltodextrin/maltose, mannitol, fructose, beta-glucosides, and unknown sugars). This was also reflected in the results of GO term enrichment analysis (Data Set S2 and Fig. 5D). The same phenomenon was not observed in the coculture with Lb. oligofermentans. Three PDE-encoding genes (pdhABCD at 3 to 11 h, alsS at 3 to 5 h, and pfl at 11 h [Fig. 3]) were downregulated in all the cocultures, but at different time points depending on the coculture and the gene. In addition, growth with Lb. oligofermentans and in the triple culture led to downregulation of several ATP synthase subunits at 3 and 11 h, respectively, and strong upregulation (log2 FC, ∼3) of aldehyde-alcohol dehydrogenase (adhE) at 11 h.
For Lb. oligofermentans, the functional enrichment analyses (Table 3, Data Set S3, and Fig. 5E) and closer investigation of DE genes during its growth in all the cocultures revealed downregulation of carbohydrate catabolism/transport (3 to 11 h) and pyruvate metabolism and energy generation (3 to 5 h), but at different time points depending on the coculture and the gene. In particular, downregulated genes included those involved in carbohydrate-specific catabolism/transport (glucose, fructose, mannose, gluconate, maltose, N-acetylglucosamine, and xylose; at 3 to 11 h), the phosphoketolase pathway (zwf, gnd, yqeC, rpe, xfp, gap, eno; at 3 to 5 h at least in one coculture), and CCC (ccpA [Fig. 3] and ptsH; at 5 to 11 h). Growth in the triple culture led to downregulation of three PDEs, pdhABCD (also with Le. gelidum), ldhD (also with Lc. piscium), and ldh, at 3 to 5 h, and three ATP synthase subunits (5 h). Genes involved in NAD(P)H reoxidation (adhE, radH, adhA, butA2, dhaT, and nox) were downregulated at different times in all the cocultures, but the effect seemed to be stronger at 5 h.
Translation.Strong enrichment of GO terms related to protein translation and peptide metabolism was observed in the Lb. oligofermentans genes upregulated at 3 h and downregulated at 11 h during growth in all the cocultures (Data Set S3 and Fig. 5F). The enrichment in all the cases was due mostly to the differential expression of the ribosomal proteins (RPs) (mostly from the cluster of LACOL_0481 to _0518 [LACOL_0481-0518]) that were upregulated at 3 h (50% of RPs in the genome) and downregulated at 11 h in all the cocultures (34 to 50%).
Even though the enrichment of the protein translation/peptide metabolism-related GO terms and the “J” category was less prominent for the interactions of Le. gelidum and Lc. piscium (Data Sets S1 and S2, Fig. 5B and D, and Tables 1 and 2), their RPs showed differential expression as well. In Le. gelidum, 16% and 35% of the RPs were mildly upregulated in the pairwise cocultures with Lc. piscium and Lb. oligofermentans at 11 h, respectively. In Lc. piscium, 23% of the RPs (mostly from the cluster LACPI_0247-0267) in the triple culture were first mildly upregulated at 5 h and then downregulated at 11 h.
Fatty acid biosynthesis.In Lc. piscium, the fatty acid biosynthesis cluster (LACPI_0829-0839) was overexpressed in all its cocultures at 11 h, and this was also reflected in the enrichment of the “Lipid transport and metabolism” (“I”) COG category (Tables 2 and S1).
Stress response and protein turnover.For Le. gelidum, downregulation of stress response and protein turnover genes during growth with Lc. piscium at 3 h was observed. This is consistent with enrichment of the DE genes in the “O” COG category (Table 1) in the same coculture and at the same time point.
In Lc. piscium, downregulation of stress response-related genes was much stronger and started earlier (3 h) for the coculture with Lb. oligofermentans and the triple culture than for growth with Le. gelidum (5 h). This observation was consistent with the observed enrichment of GO terms associated with response to (oxidative) stress/protein folding/cellular homeostasis (Data Set S2 and Fig. 5C), and the “O” and “Defense mechanisms” (“V”) COG categories (Table 2).
For Lb. oligofermentans, several stress response- and protein turnover-related genes had lower expression at 3 to 5 h (enriched GO terms included protein folding/cell redox homeostasis) and higher expression at 11 h (enriched GO terms included proteolysis/protein folding/response to oxidative stress) in all the cocultures (Data Set S3). The latter was also reflected in the enrichment of the “O” category for the DE genes at 11 h (Table 3).
More-detailed information about the downregulated stress response and protein turnover genes can be found in Note S1 in the supplemental material.
Bacteriocins, putative antibacterial peptides, and corresponding immunity proteins.In Le. gelidum and Lc. piscium, genes previously predicted to be involved in the production of bacteriocins and other antimicrobial substances (23, 32, 35) did not show differential expression during their interactions with other species. Nevertheless, two candidates for antimicrobial cationic peptides, LEGAS_1650 (50 amino acids [aa]; predicted pI, 10.29; total network charge, +9) in Le. gelidum and LACOL_0036 (46 aa; predicted pI, 9.46; total network charge, +3) in Lb. oligofermentans, which were predicted to form alpha helices and interact with membranes, were upregulated in Le. gelidum in response to cocultivation with Lb. oligofermentans at 11 h and in Lb. oligofermentans in response to cocultivation in all its cocultures at 11 h (log2 FC, ∼1.7 to 2.1), respectively. LEGAS_1650 had only one uncharacterized homolog in the UniProtKB database (as of February 2017), while LACOL_0036 had more than 80 homologs with identity of 49 to 72% with uncharacterized function found only in Lactobacillus species. Notably, LEGAS_1650 was a part of a cluster of uncharacterized proteins, LEGAS_1647-1652 (putatively secreted and involved in cell wall degradation) in Le. gelidum, which had the same DE profile. As for bacteriocin immunity, a protein (LACOL_1493) containing the OmdA domain, which is found in proteins mediating bacteriocin protection function (36), was overexpressed in all three cocultures of Lb. oligofermentans at 11 h. In addition, antiholin-like protein-encoding genes (lrgAB) that were shown to reduce penicillin-induced killing (37) were overexpressed in the coculture with Le. gelidum at all time points.
Horizontal gene transfer (HGT) and protection against it.In Lc. piscium, structural and DNA regulation genes in two of the four predicted prophages (32) were upregulated in the cocultures at 3 to 5 h. Eleven and 23 genes in the prophage LACPI_1919-1981 (63 genes) were upregulated in the pairwise culture with Lb. oligofermentans and in the triple culture, respectively, while 14 genes in the prophage LACPI_2131-2184 (54 genes) were upregulated in the triple culture. In the first prophage, DE genes included also lysis-related genes. All five restriction modification (RM) systems or their restriction enzymes were upregulated in response to Lb. oligofermentans (3 h), and four of them were so in the triple culture (3 to 5 h). Two RM systems (type I [LACPI_0877, LACPI_0881, LACPI_0884] and type II [LACPI_1774, LACPI_1776-1777]) were upregulated as a whole, while in the three other RM systems (types III, II, and I), only restriction enzymes (LACPI_0417 [Fig. 3], LACPI_0984, LACPI_2391) were overexpressed.
In Lb. oligofermentans, several phage structural genes (17 had log2 FC of ≥0.6 and 7 of them had an adjusted P of ≤0.05) of the only prophage in the genome (LACOL_0790-0844, 55 genes overall) (33) were upregulated in the triple culture at 5 h. In addition, three competence-related proteins (ComFA, ComEA, and LACOL_1261) were overexpressed in all the cocultures at 5 h. Two of four clustered regularly interspaced short palindromic repeat(s) (CRISPR) proteins (LACOL_0188 and LACOL_0191) were downregulated at 11 h in the pairwise culture with Lc. piscium. In addition, the restriction enzyme (LACOL_1094) of a type II RM system (the only one in the genome) was upregulated at 11 h in the triple culture.
None of the three predicted prophages in Le. gelidum (35) and one RM system of type I were found to be DE. Only one competence-related protein ComEA was downregulated at 3 and 11 h during growth with Lc. piscium.
Adhesins.Eleven of 22 predicted putative adhesin proteins in Lb. oligofermentans (33) were upregulated at least in the two cocultures. Nine such proteins, overexpressed in the pairwise cocultures with Le. gelidum and Lc. piscium and/or in triple culture, included eight serine-rich proteins (LACOL_0450 [3 and 11 h] and LACOL_1668-1174 [3 and 5 h]) and one mucus-binding domain-containing protein (LACOL_1109 [3 h]). The expression of the two other putative adhesins was more culture specific: fibronectin-binding protein A, encoded by fbpA, was upregulated in the pairwise culture with Lc. piscium and in the triple culture at 11 h, while LACOL_1683, containing a mucus-binding domain and an LPxTG-like motif, had higher expression during growth with Le. gelidum (Fig. 3) and in the triple culture at 11 h. In addition, the uncharacterized operon LACOL_0664-0667, putatively coding for transmembrane/cell surface proteins, was upregulated at 11 h in all the cocultures. One of these proteins (LACOL_0664) contained the DUF4097 domain, annotated with a putative adhesin function in the PFAM database (pfam13349). Even though genomes of Le. gelidum and Lc. piscium also contained putative adhesins with fibronectin-binding, mucus-binding, serine-rich, and LPxTG domains/motifs (23, 32), they were not upregulated in response to cocultivation.
A full functional description of the genes mentioned in Results can be found in Data Set S4 in the supplemental material using their gene names.
DE genes during growth of the individual cultures.To see if the cocultivation-induced changes in gene expression could be a result of cocultivation-induced growth phase transitions, we performed DE analysis between different time points (5 and 3 h, 11 and 5 h) during growth of the individual cultures of each LAB species followed by functional enrichment analyses of the DE genes (Data Sets S1 to S4 and Fig. S3 in the supplemental material). A summary of the growth-related DE changes can be found in Note S2 in the supplemental material, and they are discussed below.
DISCUSSION
Nearly all ecological niches are inhabited by numerous microbial species that interact with each other. Such interactions, along with abiotic factors, determine the population structure and global metabolic activity of a microbial community. In the case of LAB, the knowledge of interspecies interactions would be of importance for developing new technologies/strategies enabling the shelf life of food products to be extended, and also in other LAB applications, such as developing stable starter cultures and designing the composition of probiotics.
Generally, mixing different bacterial species is predicted to generate competition, while mutualistic relationships emerge quite rarely (16). Competition creates different kinds of stress conditions, such as nutrient limitation, acidification, alteration in the atmosphere and medium composition, and bacteriocin- and ROS-induced damage of membranes and DNA (16). This explains the wide array of LAB functions affected during their cocultivation in this study.
The interspecies interactions strongly affected the expression of the genes involved in carbohydrate catabolism/transport, fermentation, energy generation, and translation in all three LAB studied. However, these pathways were regulated oppositely in these LAB. The observed upregulation of carbohydrate catabolism/transport and fermentation genes in Le. gelidum in the triple culture at 11 h, when glucose concentration is expected to be much lower than at 3 to 5 h, suggests that under carbohydrate limitation conditions this bacterium responds to the presence of other LAB by enhancing its nutrient-scavenging capabilities. Expression of RPs and concentration of ribosomes were shown to be positively correlated with the growth rate of microorganisms (38–40). Even though upregulation of RPs was detected only in its pairwise cocultures (11 h) and not in the triple culture, this finding might indicate that during growth with other LAB, Le. gelidum attempts to enhance its growth capabilities under glucose limitation conditions (at 11 h). Conversely, cocultivation with Le. gelidum in pairwise cocultures and triple culture caused downregulation of the carbohydrate catabolism/transport and fermentation/energy generation pathways in the two other LAB species. The observed early (3-h) coculture-induced upregulation of RPs in Lb. oligofermentans and to a lesser extent in Lc. piscium suggests that these two bacteria attempt to increase their growth rates to potentially gain a competitive advantage already in the beginning of their interactions with other LAB. Nevertheless, the late (11-h) coculture-induced downregulation of RPs in these two LAB indicates a slowdown in their growth. This is consistent with the overall inhibition of catabolic pathways in these LAB. Such opposite adjustments of metabolic and growth capabilities in Le. gelidum and the two other LAB species may give a competitive advantage to Le. gelidum at the late growth stage and explain its ability to dominate in packaged food products at the end of shelf life, even though the initial cell counts are not that high (25, 41). The absence of such a response in the pairwise cultures of Le. gelidum suggests that the factors elicited by both Lc. piscium and Lb. oligofermentans are necessary in order to induce such a response. The coordinated up- and downregulations of glucosamine-6-phosphate degradation and biosynthesis enzymes, respectively, in Le. gelidum in the presence of Lc. piscium were similar to those obtained in Escherichia coli in response to growth on amino sugars (42). This implies that Lc. piscium might provide a source of amino sugars to Le. gelidum. We consider it likely that such a source is a peptidoglycan of the nondividing fraction of Lc. piscium cells that presumably includes also dead cells. Strong overexpression of adhE in Lc. piscium at 11 h in the coculture with Lb. oligofermentans and in the triple culture without upregulation of acetyl coenzyme A (acetyl-CoA)-producing enzymes (pdhABCD and pfl) may indicate that adhE mediates additional functions other than ethanol production, such as protection against oxidative stress (including H2O2-induced stress) (43). This possibility seems interesting, given that Lb. oligofermentans is able to produce a large amount of H2O2 both aerobically and anaerobically (33).
During competition, bacteria can employ attack strategies, which involve production and secretion of antimicrobial substances, such as bacteriocins (44). Bacteriocins are very diverse antimicrobial peptides that share little sequence homology but often contain a positive network charge allowing them to permeabilize bacterial membranes and eventually disrupt them (45, 46). They are difficult to predict both based on homology and ab initio. Since they are released toward competitor species, investigation of gene expression changes induced by microbial interactions provides an opportunity to detect putative candidates for antimicrobial peptides. In this study, two such bacteriocin candidates, LEGAS_1650 in Le. gelidum and LACOL_0036 in Lb. oligofermentans, were detected based on their overexpression in the cocultures and predicted physicochemical properties. The presence of LACOL_0036 in many Lactobacillus species and its overexpression in all the cocultures of Lb. oligofermentans (11 h) suggest that this peptide can be a competitor nonspecific bacteriocin of lactobacilli. In turn, the upregulation of LEGAS_1650 was interaction specific (only in pairwise culture with Lb. oligofermentans). Self-immunity of bacteriocin-producing bacteria is achieved through expression of the corresponding bacteriocin immunity proteins (47). The putative bacteriocin immunity protein LACOL_1493 in Lb. oligofermentans upregulated in all the cocultures at the same time (11 h) as the cationic peptide LACOL_0036 could possibly mediate self-protection from this peptide. In addition, strong and long-lasting upregulation of antiholin-like proteins (shown to be involved in the modulation of murein hydrolase activity and penicillin tolerance [37]) in Lb. oligofermentans in response to growth with Le. gelidum is an indication that Le. gelidum might produce some cell wall-damaging substance(s) against Lb. oligofermentans. One such possibility is the upregulated cell wall degradation cluster in Le. gelidum (LEGAS_1647-1652) with the putative cationic antimicrobial peptide LEGAS_1650.
In response to the hostile environment created by the competitors, bacteria may set up defense responses, such as stress response and repair mechanisms. Therefore, the early (3- to 5-h) downregulation of stress response and protein turnover genes in all LAB in response to cocultivation is interesting. Further studies are required to determine whether the presence of other LAB species alleviates the level of certain stresses in the beginning of coculture growth.
The presence of several bacterial species in the same environment creates a possibility for genetic information exchange. Whereas this may allow the community to increase its ecological fitness and can thus be beneficial for the community as a whole, it might be harmful for the individual members (48). Therefore, mechanisms for both DNA exchange and DNA protection exist. Among the LAB studied, Lc. piscium and Lb. oligofermentans showed induction of prophages (one of the mechanisms of HGT [49]) and RM systems (innate immune system against foreign DNA entry) in their mixed cultures, whereas Le. gelidum did not. The early (3- to 5-h) upregulation of the prophage genes in Lc. piscium and Lb. oligofermentans occurred solely or more strongly in the triple culture, suggesting that induction of prophage genes in these two species can be a cumulative effect invoked by the presence of more than one species. Prophages are known to be induced by different stress conditions (DNA damage, heat and oxidative stress, starvation), and upregulation of prophage genes was also shown during interactions of soil bacteria (50). The induction of only restriction enzymes and not the whole RM systems in several cases in Lc. piscium and Lb. oligofermentans makes sense, since only restriction enzymes are directly involved in protection function, while methylation of a host DNA occurs beforehand.
Good adhesive properties allow bacteria to quickly colonize surfaces and may result in a competitive advantage. The observed interaction-induced upregulation of many predicted adhesive proteins, such as serine-rich, fibronectin-binding and mucus-binding domain-containing proteins, in Lb. oligofermentans suggests that adhesion can be an important survival mechanism for Lb. oligofermentans during its interspecies interaction, and this might explain its persistence despite its low growth rate (27).
A comparison of cocultivation-induced and growth-related DE changes for each LAB indicates that some of the cocultivation-induced DE changes might be linked to the growth phase transitions. Particularly in Le. gelidum, upregulation of the carbohydrate catabolism and fermentation in the triple culture may result from cocultivation-induced promotion of its transition to the later growth phase. This is suggested by the plateauing of its growth curve as a result of cocultivation with other LAB and upregulation of similar pathways during its growth in the late period (5 to 11 h). Even though no clear cocultivation-induced changes in growth curve shapes were observed for Lc. piscium and Lb. oligofermentans, some metabolic changes (e.g., downregulation of carbohydrate catabolism and fermentation and upregulation of the fatty acid biosynthesis in Lc. piscium; downregulation of carbohydrate catabolism and fermentation and upregulation of some adhesins in Lb. oligofermentans) during growth in their cocultures may be explained by the cocultivation-induced delay of their growth transition to the later stage.
To conclude, the study revealed different growth and survival strategies employed by three meat spoilage-associated LAB species to cope with the consequences of interspecies competition. In the presence of other LAB, the fastest growing LAB, Le. gelidum, attempted to access carbon sources more efficiently and enhanced its growth potential under glucose limitation conditions by the upregulation of carbohydrate catabolic and fermentation pathways as well as RPs. The opposite responses were observed for the slower-growing LAB Lc. piscium and Lb. oligofermentans. Such behavior of Le. gelidum may partly explain its documented rapid outgrowth from small initial cell counts (25) and predominance in spoilage microbial communities in a variety of food products. Le. gelidum may also be able to use amino sugars, presumably released from the peptidoglycan of dead cells of other bacterial species, as an alternative source of carbon and nitrogen, and this may give a competitive advantage to Le. gelidum in food microbial communities. Finally, cocultivation-induced overexpression of the numerous putative adhesins in Lb. oligofermentans could contribute to its survival in actively growing communities in meat despite its low growth rate. Even though the study cannot be directly generalized to the “real-life” situations occurring in meat products, which have a more complex nutrient composition and more-diverse microbial communities, these findings provide hypotheses for future studies conducted with more strains and in a real food environment.
MATERIALS AND METHODS
LAB cocultivation.Three LAB, Le. gelidum subsp. gasicomitatum LMG 18811T, Lc. piscium MKFS47, and Lb. oligofermentans LMG 22743T, were grown separately in individual cultures and in mixed cultures pairwise (three combinations) and all together (here referred to as triple culture). The experiments were performed in two batches. First, cocultivation of Le. gelidum and Lc. piscium accompanied with their individual cultures was performed and processed. Later, two other pairwise cultures (Le. gelidum with Lb. oligofermentans and Lc. piscium with Lb. oligofermentans) and the triple culture were grown together with the individual cultures of all three LAB. For all the growth experiments, each LAB was precultured separately as follows: one colony from a plate was suspended in 9 ml of de Man-Rogosa-Sharpe (MRS) liquid medium without acetate and grown for 48 h, followed by reinoculation of 0.1 ml of the obtained culture into 9 ml of the same fresh medium, where it was grown for 12 h (Le. gelidum and Lc. piscium) or for 24 h (Lb. oligofermentans) to reach approximately equal OD600 values for all three LAB. The obtained precultures were inoculated into 250 ml of fresh MRS broth without acetate individually or as mixtures, with the total amount of inoculum being 0.1 ml and equal ratios of each species in mixed cultures, therefore, 0.1 ml of a preculture for individual cultures, 0.05 ml of any two precultures for pairwise cultures, and 0.033 ml of all precultures for the triple culture. The obtained cultures were grown (50 rpm) until they reached an OD600 of ∼0.2 to 0.4 and were then diluted to an OD600 of ∼0.02 to 0.04, the total amount being 500 ml. The cultures obtained were grown for 11 h under microaerobic conditions (50 rpm) in three replicates. The OD600 values were measured between 0 and 11 h of growth. The CFU per milliliter was determined by plating 10-fold serial dilutions on an MRS medium without acetate under anaerobic conditions after 5 days. While measuring OD600 accounts for all types of cells in a sample (dividing, dormant, and dead) and is generally well correlated with the dry weight of cells regardless of their size (34), the CFU/milliliter-based method measures the concentration of actively dividing cells only. The viable cell counts of each species in the cocultures were determined by plating on species-selective media: MRS medium with fusidic acid (30 μg/ml) for Le. gelidum, glucose-supplemented M17 (GM17) medium with polymyxin (60 μg/ml) for Lc. piscium, and MRS medium with kanamycin (85 μg/ml) and xylose instead of glucose for Lb. oligofermentans. All cultures were grown at 25°C.
Reverse genetics assays with these species have been proven to be difficult. Knockouts have not been obtained by crossover recombination even after excessive trials.
RNA extraction and sequencing.Samples for RNA extraction were taken at three time points, 3, 5, and 11 h, for all the cultures. RNA extraction and library construction were done in a manner analogous to that of our previous study (32). Briefly, cell metabolism and RNase activity were inhibited by treatment with a cold 10:1 mixture of ethanol-phenol. The samples were centrifuged and immediately frozen. The cells were disrupted by a mixer mill, and RNA was extracted by the RNeasy plant minikit (Qiagen) with DNase treatment. rRNAs were omitted from the total RNA, and RNA-seq libraries were produced using Ovation RNA-seq system V2 (NuGEN). The obtained double-stranded DNA (dsDNA) was sheared using sonication, followed by the polishing of the obtained fragments with T4 DNA polymerase and the ligation of adapters for SOLiD sequencing. Libraries, numbering 81 overall [(31st batch + 62nd batch) cultures × 3 time points × 3 replicates], were size selected and thereafter sequenced in five lanes using SOLiD 5500XL (Life Technologies, Foster City, CA, USA) to produce 75-bp single-end reads.
RNA-seq data preprocessing and differential expression analysis.RNA-seq reads obtained from all samples were mapped using Lifescope software (Life Technologies), specifically designed for SOLiD RNA-seq data, against all three genomes of the LAB at the same time to allow proper map and count reads that align to more than one genome. During the counting process of coding DNA sequence (CDS) mapped reads performed by Lifescope, a read with alignments of different qualities to more than one genome/genomic locations is assigned to the location where it has the highest alignment quality or is discarded. Mapping quality takes into account the length of alignments, the number of possible alignments, and the number of mismatches (51). The overall read-mapping percentages varied from 69% to 95% for different samples. Hierarchical clustering of samples based on Euclidean distances and principal-component analysis (PCA) of samples were done for each bacterium as described in the manual for the DESeq2 package (52) on “regularized log” (rlog)-transformed read count data to visually explore sample relationships. Based on visual investigation of the obtained clusters and PCA plots, five samples were discarded as potential outliers, as they were clearly separated and most distant from all other samples in the same batch and at the same time point (for more details, see Table S2 and Fig. S4 in the supplemental material). In addition, two samples contained a very low number of Lb. oligofermentans CDS-mapped reads (∼0.04 million and 0.07 million). To avoid introducing higher variability to the data, especially for low-expressed genes, these reads were not used for further analysis. The number of CDS-mapped reads was in the range of 0.29 to 9.9 million for Le. gelidum, 0.18 to 8.5 million for Lc. piscium, and 0.17 to 8.0 million for Lb. oligofermentans per sample for each genome (see Table S3 in the supplemental material). Hierarchical clustering revealed a strong batch effect, and so differential expression (DE) analysis was performed independently for each batch.
DE analysis for each LAB species was performed between RNA-seq reads mapped to its CDSs in the individual culture (serving as a control) and in its mixed cultures (two pairwise and one triple for each species) using R package DESeq2 (52) (see Fig. S5 in the supplemental material). This was done independently for each of the two batches and for each time point (3, 5, and 11 h). In addition, differential expression was tested between different time points (5 h over 3 h and 11 h over 5 h) for the individual cultures of each LAB species controlling for the batch effect. Default parameters were used. Genes (CDSs) that had very low raw read counts (average per sample across all samples for each bacterium, <1 read) were considered not expressed and were filtered out before DE analysis. Genes were considered to be differentially expressed if their adjusted P value was less than or equal to 0.05 and their log2 FC, obtained from the DESeq2 output, was higher than or equal to 0.6 (therefore, FC was higher than or equal to ∼1.5). Information on the processed RNA-seq data (raw read counts) and DE analysis (adjusted P values and log2 FC) can be found in Data Sets S5 and S4 in the supplemental material, respectively.
Droplet digital reverse transcription-PCR.For absolute gene expression quantification using dd-RT-PCR, primers (see Table S4 in the supplemental material) were designed using Primer3Plus (53) and manufactured by Integrated DNA Technologies. Prior to reverse transcription (RT), residual genomic DNA (gDNA) was removed from the RNA samples (600 to 800 ng) by the Heat&Run gDNA removal kit (ArcticZymes). One half of the RNA samples were used to synthesize cDNA using 200 ng random hexamer primer (Thermo Fisher Scientific), 40 units RiboLock RNase inhibitor (Thermo Fisher Scientific), and 200 units SuperScript III reverse transcriptase (Invitrogen). The RT reaction mixtures were incubated at 50°C for 60 min. The other half of the RNA samples were used for no-RT control assays, in which water was added instead of the RT enzyme. dd-PCRs included 1× QX200 ddPCR EvaGreen Supermix (Bio-Rad), 200 nM both primers, and 5 μl of cDNA or no-RT controls in a total volume of 22 μl. Template cDNA/no-RT reaction mixtures were diluted either 1:200 (single-copy gene assays) or 1:1,000,000 (16S rRNA assays). Droplets were generated using the Automated droplet generator (Bio-Rad). Thermal cycling conditions were 95°C for 5 min, 40 cycles of 95°C for 30 s and 58°C for 1 min, and then 4°C for 5 min and 90°C for 5 min. Amplification was performed in a Veriti thermal cycler (Applied Biosystems) using a 50% ramp rate and a 40-μl sample volume. The fluorescence of each droplet was measured by a QX droplet reader (Bio-Rad), and the results were analyzed using QuantaSoft software (version 1.7.4.0917; Bio-Rad). Thresholds to separate PCR-positive and PCR-negative droplets were set manually. To be able to compare expression levels in different samples, concentrations of the target genes (cDNA copies/microliter) were normalized by the concentrations of species-specific genes: infB, ftsQ, and 16S rRNA. Genes infB and ftsQ are housekeeping genes, and their expression was found to be mostly unchanged across different cultures and time points for all three LAB species.
Functional enrichment analyses of the DE genes.Genes (based on protein translations) in the genomes were classified into COG functional categories by an RPS-BLAST search against the COG database (March 2016, with an E value threshold of 0.01 and one best hit taken into account for each gene) (Data Set S4). GO terms of the genes were predicted using the Blast2GO platform (54) with default parameters based on both BLAST and InterPro searches, where the BLAST search was performed against the Swiss-Prot database (Data Set S4) (November 2016). The databases of the metabolic pathways were constructed using Pathway Tools software (“pathological” component) (55) for Lc. piscium and Lb. oligofermentans or downloaded from the BioCyc PGDB (pathway/genome database) registry for Le. gelidum (LGAS762550cyc, version 16.1). The lists of DE genes obtained were tested for the enrichment of the COG functional categories, GO terms (belonging to the Biological Process ontology), and metabolic pathways using one-tailed Fisher's exact test implemented in an R software environment (using the function fisher.test), topGO package (56), and Pathway Tools software (55), respectively. Enrichment for GO terms and metabolic pathways was tested separately for up- and downregulated genes. The reference set included all the annotated (in terms of COG categories, GO terms, and pathways, respectively) genes in the genome. In the cases of the COG categories and GO term enrichment analyses, only expressed genes (see “RNA-seq data preprocessing and differential expression analysis” above for the criteria for nonexpressed genes) were included in the reference set. The functional categories/terms/pathways were considered to be enriched if their P value of enrichment was less than or equal to 0.05. For the COG categories enrichment analyses, the P values obtained were also adjusted for multiple comparisons using the Benjamini and Hochberg correction method (1995) implemented in the R function p.adjust. TreeMaps of the enriched GO terms were built using the REVIGO Web server (57) (the database with GO term sizes for Bacillus subtilis was used).
Functional gene annotation.Isoelectric points for the peptides were predicted using the Expasy Compute pI/Mw tool (http://web.expasy.org/compute_pi/) (58). Prediction of antimicrobial peptides was done using APD3, antimicrobial peptide calculator and predictor (http://aps.unmc.edu/AP/prediction/prediction_main.php) (59).
Accession numbers of the RNA-seq data.Raw RNA-seq read sequences are available in the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-5499, except for the reads from the individual culture of Lc. piscium in the batch 1, which are available under accession number E-MTAB-3245.
ACKNOWLEDGMENTS
Authors' contributions are as follows: M.A. performed RNA-seq data preprocessing and DE analyses, functional enrichment analysis, and interpretation of the results and wrote the manuscript. E.J. and P.J. designed and performed the growth experiments. A.Y. performed dd-RT-PCR, including the design of primers. L.P. planned and performed RNA-seq. J.B. and P.A. designed the experiments and participated in the interpretation of the results and the writing of the manuscript.
We appreciate the excellent technical assistance of Erja Merivirta, Henna Niinivirta, Eeva-Marja Turkki, Kirsi Lipponen, Harri Kangas, and Päivi Laamanen. Timo Nieminen is thanked for useful comments concerning the manuscript. Last but not least, we express our gratitude to Stephen Skate for the language editing of the manuscript.
This project was funded by the Academy of Finland (grant numbers 267755, 311717, 307856, 307855, and 267623), Integrative Life Science Doctoral Program (ILS), and Institute of Biotechnology. The funders had no role in study design, data collection and interpretation, or in the decision to submit the work for publication.
FOOTNOTES
- Received 15 March 2018.
- Accepted 12 April 2018.
- Accepted manuscript posted online 20 April 2018.
Supplemental material for this article may be found at https://doi.org/10.1128/AEM.00554-18.
- Copyright © 2018 American Society for Microbiology.
