Introduction

Several studies describe the lactic acid bacteria’s (LAB) ability to survive and respond to various environmental stresses (Tsuda et al., 2019; Bucka-Kolendo and Sokołowska, 2017) as they acquire preserving, probiotic, or spoilage properties (Bangar et al., 2022; Zapaśnik et al., 2022; Han et al., 2015). The most frequent phenotype described in the literature is the adapted cell (Papadimitriou et al., 2016), while adaptation refers to the effort of the cell to resist and persist under stress. However, the stress response to environmental factors can differ between species and depend on the applied stress (Mahmmodi et al., 2021; Van de Guchte et al., 2002). Extensive studies revealed the mechanisms involved in heat shock (De Angelis et al., 2004), bile (Bron et al., 2006), oxidative (Serrano et al., 2007), pH, and ethanol (Parente et al., 2010), where bacteria, through activation of the mechanisms involved in the stress response, adapt to the new conditions (Bucka-Kolendo and Sokołowska, 2017). Although the activated mechanisms may partially overlap, they are not identical (Papadimitriou et al., 2016), which can cause opposed results in LAB species, at species or even sub-species level.

Furthermore, the combination of stresses can trigger a cross-protection response (Yang et al., 2021b; Bucka-Kolendo and Sokołowska, 2017; Papadimitriou et al., 2016; Van de Guchte et al., 2002). The molecular mechanisms underlying the adaptation potential and response are based on the coordinated gene expression that can affect processes in cells, like cell division, transport, membrane composition, and DNA metabolism. Therefore, gene transcription, expression levels, and mechanisms engaged in bacteria growth under diverse stress conditions are greatly valued. Lactobacillus is a significant group of widespread organisms in different environments.

High hydrostatic pressure (HHP) is not a common stress factor for the LAB, as they are not generally exposed (Bucka-Kolendo and Sokołowska, 2017). However, HHP is a popular nonthermal preservation technique used in the food and beverage industry that reduces the number of microorganisms (Yaman et al., 2020; Chen et al., 2016) while preserving organoleptic molecules and providing “fresh” food. There is still limited knowledge about how LAB responds to the HHP, and its molecular mechanisms are not fully understood. HHP can negatively affect all molecular mechanisms in bacteria where DNA is involved, such as replication, transcription, and recombination (Salvador-Castell et al., 2020). The interaction between DNA and proteins may be disturbed due to the changes in the electrostatic and hydrophobic interactions. Pressure can dissociate ribosomal subunits and disturb the cytoskeletal proteins, resulting in reversible morphological changes. Since some effects of different factors can be similar, it is assumed that the ability to react to HHP comes from the cross-protection system HHP due to the fluidity of a complex response mechanism. Membrane fluidity among the critical factors is responsible for bacteria’s survival and growth under high-pressure conditions (Molina-Hoppner et al., 2003), where membrane lipids are stabilized by HHP and increase the melting points of lipids and transition the lipid bilayer to the gel state. Bacteria adapted to the HHP can adjust the phospholipid composition of the membrane by increasing the quantity of unsaturated fatty acids.

In bacteria, 90% of the genome represents genes; the rest contains small intergenic regions occupied by regulatory sites (Rocha, 2018). Genes are usually organized into operons and have only a few or no introns, and the insertions and deletions of genetic material tend to be determinants of gene expression (Rocha, 2018; Price et al., 2006). Many phenotypic variations among species are assigned to single nucleotide polymorphisms (SNPs) (Bailey et al., 2021; Hunt et al., 2009). Within the population, single base changes occur with a frequency greater than 1%. SNPs can be either synonymous (SS) when they do not cause changes in the amino acids (AA) or non-synonymous (NSS) when the AA structure is altered (Bailey et al., 2021; Lebeuf-Taylor et al., 2019; Hunt et al., 2009). It has been demonstrated that mechanisms altering the proteins’ structure, function, and expression level by affecting mRNA splicing, stability, structure, and protein folding are now better understood (López-González et al., 2018). NSS is more frequent and has a more substantial effect than SS mutations (Rocha, 2018). Since the effect is usually negative, those mutations are progressively removed from the population by the natural (purifying) selection, causing the low ratio of nonsynonymous (Ka) to synonymous (Ks) substitutions (ω - Ka/Ks). At the same time, synonymous substitutions can also be influenced by purifying selection, especially in fast-growing bacteria (Rocha, 2018). Many experimental studies prove that SS mutations can have positive solid fitness effects and drive adaptive evolution (Bailey et al., 2021; Liu et al., 2019).

It is crucial to provide insight into the genomic dynamics and polymorphism that characterize the physiological state of cells after exposure to stress, such as HHP, and understand the relationship between gene functions and phenotypic characteristics. Five LAB strains (two Loigolactobacillus backii, two Lactiplantibacillus plan-tarum, and one Lacticaseibacillus rhamnosus) were used. Strains were previously identified and analyzed (BuckaKolendo et al., 2020, 2021). This work described the proteomic and transcriptomic changes in selected Lactobacillus treated with HHP to determine the general adaptive response resulting from pressurization. As in previous studies, the HHP affects the proteome of treated LAB strains causing differences in the mass spectra profiles analyzed with MALDI-TOF MS (Bucka-Kolendo et al., 2020). Pressurization of the LAB strains can lead to changes in the expression patterns of stress-related genes (Bucka Kolendo et al., 2021). Among many functional genes associated with adaptation in LAB (Bucka-Kolendo et al., 2017), the three genes (dnaK, hrcA, and ctsR) previously described in the literature as stress-related (BuckaKolendo et al., 2017), were selected. The genes, like dnaK, hrcA, and ctsR, involved in the stress response (BuckaKolendo et al., 2021) can have different expressions under the stress factor, resulting in changes in the phenotype.

The aim was to relate the effect of the HHP on strains’ fitness through the changes in those gene expressions using SNP analysis of the partially sequenced genes. The hypothesis was that differences due to HHP might contribute to the bacteria in their gene mutations associated with the stress response and be elucidated with the genes’ phylogenetic clusters. The phenotypic and phylogenetic characterization of the LAB strains allows rising new insights into the adaptive abilities under the HHP. The fitness changes of the cellular response to stress factors can give an understanding of the individual strains’ responses to different factors (Douillard et al., 2016). As LAB are essential probiotics, starter, commensal, and pathogenic microorganisms, the in-depth research about the physiology of LAB stress has a significant meaning.

Materials and Methods

Schematic overview of the experiment

The schematic workflow of the study is presented in Figure 1. The diagram shows the overall process for the observation of the changes in the nucleotide sequences in the stress-related genes dnaK, ctsR, and hrcA under the HHP 300 MPa/5’. LAB strains were screened for possible adjustments to the changing environment. The first step was the isolation of strains from spoiled food products and the identification of the bacteria. The second step consists of applying HHP, molecular analysis with stress-related gene sequencing and expression, and growth analysis. The final step include the statistical and analytical analysis of SNP, phylogenesis, and PPI (protein-protein interaction).

Lactic acid bacteria

Five LAB strains (KKP 3565 Loigolactobacillus backii, KKP 3566 Loigolactobacillus backii, KKP 3568 Lactiplantibacillus plantarum, KKP 3569 Lactiplantibacillus plantarum, and KKP 3570 Lacticaseibacillus rhamnosus) were isolated from the food products, beer, tomato juice, and bread, and the process was performed according to the ISO 15214:2000 as described by Bucka-Kolendo et al. (2021). In addition, bacteria were identified using genetic and proteomic methods, as defined by Bucka-Kolendo et al. (2020, 2021). The isolated strains were deposited in the Culture Collection of Industrial Microorganisms-Microbiological Resource Center (IAFB, Warsaw, Poland), supported by the European Horizon 2020 research and innovation programme under grant agreement No 871129-IS_MIRRI21 Project. Strains were given the collection numbers, and the 16S rDNA sequences of each strain were deposited in the GenBank NCBI database (Table 1).

Bacteria growth conditions and count number analysis (PCA)

LAB were grown and incubated under anaerobic conditions at 30°C for 48–72 h on MRS agar (Lactobacillus Agar DeMan, Rogosa, and Sharpe, Merc KGaA, Darmstadt, Germany), as described by Bucka-Kolendo et al. (2021). For counting, plates containing less than 300 CFU/mL were selected according to ISO 4833-1:2013. The non-treated LAB were considered control strains.

Immediately after HHP processing, the viability of the LAB strains was evaluated by counting colony-forming units on MRS Agar. The difference between control and treated strains was counted based on the number of surviving bacteria. The analysis was performed in two replicates.

Application of HHP

Using U 4000/65 apparatus (Unipress, Warsaw, Poland), the stationary phase LAB were threatened with HHP, as mentioned previously by Bucka-Kolendo et al. (2020). Therefore, exposition to 300 MPa for 5 min was chosen based on the former analyses (data not shown). At the HHP higher than 300 MPa, a significant reduction of live cells was observed, and at lower HHP, there were no significant differences from the controls. Hence, the used parameters were selected to observe possible cell changes but did not cause the bacteria’s death.

The HHP chamber details were as follows: the 0.96 L working volume, 600 MPa of maximum working pressure, from −10°C to +80°C working temperature, and the pressure-transmitting fluid was (1:1, v/v) water-polypropylene glycol. The time needed to generate the 300 MPa pressure was 70–80 s, with a release time of 2–4 s. The pressurization times did not include the come-up and come-down times. The assays were performed under two independent processes, and unpressurized strains were used as a control.

DNA extraction

According to the instructions, DNA from the stationary phase of control and HHP-treated bacterial culture (1.5 mL) was extracted with the ExtractMe DNA Bacteria Kit (Blirt S.A.–DNA, Gdansk, Poland) as mentioned by Bucka-Kolendo et al. (2020). The quality and concentration of the obtained DNA were measured at the absorbance of 260 nm and the 260/280 nm ratio with a UV-Vis NanoDrop spectrophotometer (ThermoFisher Scientific, Waltham, MA, USA). The isolated DNA was stored at −20°C.

PCR amplification

PCR amplification of dnaK, ctsR, and hrcA sequences was performed using specific primers (Table 2). For dnaK, ctsR, and hrcA, a total PCR volume of 60 µL contained: 30 µL of Dream Taq PCR Master Mix (ThermoFisher Scientific, Waltham, MA, USA), 1 µl of each primer in the final concentration of 0.4 µM, and 10–20 ng of DNA. A peqSTAR 2X thermocycler (PeqLab, Germany) was used for the amplification run. PCR amplification was performed as described by Riccardi et al. (2012), with modifications. Reactions conditions were: initiating denaturation at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 60 s, annealing at 58°C for 60 s (for dnaK and hrcA), annealing at 57°C for 60 s (for ctsR), and elongation at 72°C for 90 s, with final elongation at 72°C for 2 min.

Electrophoresis in 1.5% (w/v) gel agarose was performed to analyze the PCR product size, and BioImaging Systems 06-2d.1-G: BOX (Syngen, UK) was used to visualize the product. Sequencing of the PCR products was performed in the 96-capillary 3730xl DNA Analyzer (Applied Biosystems-Life Technologies), and the obtained sequences were analyzed in NCBI BLAST databases. The obtained sequences were used to evaluate the DNA polymorphism, and the results are shown in Tables 3–5, respectively.

Bacterial mRNA extraction and RT-qPCR analysis

The total mRNA extraction and relative gene expression analysis were performed according to Bucka-Kolendo et al. (2021). The custom TaqMan gene expression assays (ThermoFisher Scientific, TFS) were used for dnaK, ctsR, and hrcA. The relative expression levels of analyzed genes were standardized to an endogenous control 16S rRNA gene. Endogenous control’s stability was evaluated for control and pressurized LAB using the ΔCT algorithm. For relative expression ratios in target genes, the 2^−ΔΔCT method was used. The analysis was performed as the mean of the two independent experiments.

Statistical data analysis

The sequences of each gene (dnaK, hrcA, ctsR) of each strain (control and pressurized) were trimmed, aligned, and analyzed. Multisequences’ alignment and phylogenetic analyses were performed using MEGA v. X (Kumar et al., 2016). Phylogenetical trees were created with the neighbor-joining evolution method based on sequences of 3 stress-related genes in control and pressurized strains. Evaluation of the number of polymorphic sites (S), nucleotide diversity (π) of the strains of variant obtained for control and HHP conditions, and the ratio of the nonsynonymous (Ks) to synonymous (Ka) mutations (ω) were calculated with the Dna SP. 5.1 (Rozas et al., 2017).

To evaluate known and predicted interactions between dnaK, ctsR, and hrcA proteins in Lactobacillus, the PPI network was created with the STRING database (string-db.org) (Szklarczyk et al., 2019). The PPI provided a new possibility for revealing molecular mechanisms.

Results

The DNA polymorphism analysis was performed to obtain preliminary insight into the effect of the HHP on the LAB and identify single mutations in stress-related dnaK, ctsR, and hrcA, (Tables 3–5, respectively). The gene sequences had a G+C% content of 46.5–47.8% (dnaK), 42.5–43.7% (ctsR), 45.6–49% (hrcA). The rate of Ka to the rate of Ks was calculated to determine the evolutionary pressure on protein-coding sequences. The Ka/Ks ratios (ω) for total control populations and total HHP treated vary between 0 and 1.12, with the number of mutations lower for pressurized strains (from 6 (hrcA) to 11 (dnaK)) than controlled (from 3 (ctsR) to 92 (hrcA)). ω ratios on the strains level were 0 to 2.678, where most were 0 or close to 1, suggesting little differences between synonymous and nonsynonymous mutations. DNA polymorphism analysis of bacteria under HHP suggests that strains are under neutral purifying selection pressure and acts on the genes in most cases. In all pressurized strains, the expression of genes was observed (Figures 2–4 respectively, for dnaK, ctsR, and hrcA), whereas, in control strains, the gene expression was detected in three of five strains. In most LAB, the SNPs do not significantly change gene expression or gene product function. Considering that all strains survived the pressurization (Bucka-Kolendo et al., 2021), it was assumed that although a low SNP level occurred in the genes, bacteria gained adaptive ability. However, defining the correlation between gene expression changes and the SS mutations’ fitness mutations was impossible. The expression levels for dnaK, ctsR, and hrcA under the 300 MPa/5’ did not correlate with the SNPs in those genes.

PCA analysis

Our experiment showed that the pressurization in 300 MPa’5 significantly affected the cell’s survival (P < 0.05). For strain KKP 3570 Lacticaseibacillus rhamnosus, the decrease level was 2 log (CFU/mL). For other strains, the reduction was around 1 log (CFU/mL) (Figure 5).

Phylogenetic analysis of dnaK, ctsR, and hrcA

Maximum likelihood clustering of aligned sequences of each gene was performed to assess the changes in the relationship of genes across LAB for controlled and HHP treated strains. Alignment was performed using CLUSTAL W, and results for unrooted neighbor- joining phylogenetic trees are shown in Figures 6–8 for dnaK, ctsR, and hrcA, respectively. In addition, control and pressured strains were compared for each gene to express whether mutations were reflected in the phylogeny.

All three genes were examined to display the discriminatory power and reproducibility in studied lactobacilli strains. For ctsR, the similarity was at 97% that mirrored in the trees’ topology of controlled and HHP-treated strains. For dnaK sequences obtained, tress had a high level of resemblance at 98% for controlled and 91% for pressurized strains. The strains grouped based on hrcA sequences showed the highest differentiation, wherein control strains L. backii KKP 3565 and KKP 3566 had 97% and were more distantly related (42%) to other studied strains. The hrcA sequences of strains after pressurization were highly related and had a 100% level of similarity. There was no correlation between the expression pattern and clustering of the strains based on the SNP sequences analysis, resulting from the small number of bacteria used or the complexity of the stress-response mechanisms in LAB. Phylogenetic trees of all three genes (dnaK, ctsR, hrcA) revealed a high degree of relatedness between LAB strains.

Protein-protein interactions and genes co-expression analysis

To forecast the network of PPI, the Search Tool for the Retrieval of Instances of Neighbouring Genes—STRING database was used (Szklarczyk et al., 2019). The PPI network is valuable for describing the molecular processes, and atypical PPI can be related to the many stress response changes. Therefore, genes with a high degree of connectedness were clustered into PPI based on the highest confidence of 0.9 and a maximum number of interactions ≤ 5. The study used a k-means clustering method with an average local grouping coefficient of 0.696. As a result, the network found three clusters containing eight nodes with an average node degree of 5.

Also, the co-expression analysis of the studied genes was performed with STRING (Figure 9). Results indicate the strength of data shown in experiments, describing the correlation of expression between two coding protein genes based on the patterns of RNA expression and protein coregulation. The PPI network analysis can give a network of molecular interactions formed between the protein products of the studied genes. As presented in Figure 9A, among studied genes, the highest degree of centrality had gene dnaK, indicating the importance of those proteins in the resulting PPI network. Figure 9B shows the functional relationship retrieved for genes belonging to the network, where the color indicates the similarity that correlates with the presence/absence of the interactions. The analysis suggests that the proteins can be biologically linked and interact more than expected (20 edges compared to 5 estimated).

Discussion

The present study investigated the polymorphism in three genes (dnaK, ctsR, hrcA) related to the HHP stress response in five LAB strains. The number of SNPs in our studies was low and did not display one specific outline of the tested conditions. Although SNPs and mutation rates of all three genes were comparable in all strains, it was observed that HHP trigged a strong response in LAB strains and induced stress-related genes expression (Figures 2–4), which can determine the HHP as a potent stress factor of those genes. Furthermore, our results showed the highest nucleotide diversity in the hrcA sequences (Table 5). However, we did not find significant diversity between control and pressurized LAB strains in dnaK, ctsR, and hrcA sequences (Tables 3–5, respectively).

Overall, the LAB strains adapt to the given pressure, as their decrease level was 2 log (CFU/mL) (Figure 5). Comparing the control with the pressurized strains, the impact of the HHP on the LAB strains based on the Ka/Ks ratio displays a low value. This suggests that NSS mutations have a small negative effect on fitness and can stay in the population long before being removed by natural selection. The diversity of the dnaK, ctsR, and hrcA sequences indicates small changes between studied LAB strains, which shows that strains were protected from functional mutations, or mutations that have occurred were almost neutral. These findings are in line with the reports of other authors (Bailey et al., 2014, 2021; Rocha, 2018). An observed high rate of nonsynonymous mutations suggests that adaptive selection favors different protein sequences depending on environmental changes. They can alert the function of the genes. Our studies indicate that the expression of the genes can be the defense mechanism of strains against HHP. However, among genes, dnaK, ctsR, and hrcA expression does not seem to play a significant role in stand bias and demonstrates genomic plasticity. As Bailey et al. (2014) noticed, the mutations indicate increased gene expression. However, in the less-performed codons, the molecular mechanisms are responsible for transcript regulation and are very important in LAB’s adaptation and evolutionary dynamics. The fitness effect of a SS mutation is not always due to the impact on the gene in which it occurs but rather through the changes in gene expression which has been studied in our previous work (Bucka-Kolendo et al., 2021).

Little is known regarding the fitness effect of these mutations and how they rely on stress-related factors like HHP. The changes in the environment, such as pressurization, may trigger the activation of adaptive mechanisms and thus maintain changes in the genome, proteome, transcriptome, or metabolome, leading to cell damage and death. Studies (Jeon et al., 2021) on SNPs related to cell wall synthesis in L. acidophilus under high temperatures suggested that changes in SNP can make the cell wall more rigid. Mutations that appear under natural selection increase the possibility of fixation and may often be adaptive. The cell morphology and division can be affected by HHP and impact the cytoskeleton and autolytic proteins (Molina-Hoppner et al., 2003). As Yang et al. (2021a) noticed, after the HHP treatment of 300–400 MPa, L. plantarum and L. curvatus acquire higher cell integrity, smoother cell surface, and uniformed protoplasm. These findings suggest that bacteria can increase their HHP resistance by modifying their structures.

The stress factor, like HHP, may impact the discriminatory power of the studied genes’ sequences. Thus, comparing the DNA polymorphism changes in stress-related genes with the phylogenetic trees assembled is required and can give a comprehensive understanding of the relationship between functional genes and phenotypic characteristics. The phylogenetic trees are highly valued tools that demonstrate the diversity of bacterial strains to develop more effective methods for their identification, prediction of gene functions, and underpin genetic research. The phylogenetic analysis of the unrooted trees demonstrates the relative relationship between strains and the impact of the HHP. Based on the maximum likelihood of aligned sequences (dnaK [Figure 6], ctsR [Figure 7], and hrcA [Figure 8]), the clustering did not reflect the phylogenetic relationships. It did not allow the separation into different clusters for studied LAB species. All strains had a high similarity. The PPI network (Figure 9) showed an insight into the molecular processes, and the functional relationship retrieved for the genes belonging to the network suggests that the proteins interactions can be biologically linked more than expected, and the whole process is more complicated. Although the studies showed significantly higher diversity of the hrcA sequences, it did not correspond with strain classification and their phylogenetical clustering. The highest diversity was obtained for the hrcA gene in control strains. Our studies confirmed the results of Guidone et al. (2015), on the alignment of the ctsR sequences, marking the gene’s taxonomic value for LAB classification. However, our results disagree with their results and that classification of the LAB species based on the hrcA gene sequence is a good taxonomic marker, as the hrcA sequences after pressurization had a 100% level of similarity. Contrary to the studies presented by Huang and Lee (2011), the dnaK gene is much more polymorphic than the housekeeping gene 16S rRNA in the LAB group and therefore has a discriminatory value in closely related species. Where Sharma et al. (2019), indicated in their research that, among other genes, dnaK showed no SNP while tracking in yogurt and probiotic powder, which confirmed the results of our work, where dnaK was characterized by low diversity.

Our current work, together with previous studies on proteome changes under HHP (Bucka-Kolendo et al., 2020), and changes in gene expression under HHP in LAB (Bucka-Kolendo et al., 2021), gives a more in-depth insight into the selected LAB strains stress response to HHP. Our findings are significant to the technological implications of LAB resistance to the food industry, host survival (in the case of probiotics), and bacteria stress behavior, where a complex regulatory network of genes regulates the bacteria’s response to HHP during food processing. Further studies on transcriptome and proteome are needed to confirm this hypothesis.

Additional studies on the larger group of Lactobacillus strains are needed to examine the impact of the HHP on the stress response mechanisms of those bacteria and how the HHP can contribute to the mutations in stress-related genes. As the mutations, particularly in the promoter region of the genes related to the stress, can have important implications for the fitness effect and the adaptive ability of the strains. The performance of the whole genome sequences (WGS) to determine the evolution of LAB genomes under HHP treatment would help to obtain knowledge of the bacteria with higher flexibility under stress conditions. With WGS methods, more insights can be gained, as it combines the determination of the strain’s similarities based on SNP and gene-by-gene approach. Thus, considering how the mutation impacts the whole strain, additional studies that affect other neighbor genes, especially those grouped in the operons and transcribed together (like dnaK and hrcA), can help understand SNP mechanisms adaptative evolution better.