Comparative Genomic Analysis Determines the Functional Genes Related to Bile Salt Resistance in Lactobacillus salivarius

Lactobacillus salivarius has drawn attention because of its promising probiotic functions. Tolerance to the gastrointestinal tract condition is crucial for orally administrated probiotics to exert their functions. However, previous studies of L. salivarius have only focused on the bile salt resistance of particular strains, without uncovering the common molecular mechanisms of this species. Therefore, in this study, we expanded our research to 90 L. salivarius strains to explore their common functional genes for bile salt resistance. First, the survival rates of the 90 L. salivarius strains in 0.3% bile salt solutions were determined. Comparative genomics analysis was then performed to screen for the potential functional genes related to bile salt tolerance. Next, real-time polymerase chain reaction and gene knockout experiments were conducted to further verify the tolerance-related functional genes. The results indicated that the strain-dependent bile salt tolerance of L. salivarius was mainly associated with four peptidoglycan synthesis-related genes, seven phosphotransferase system-related genes, and one chaperone-encoding gene involved in the stress response. Among them, the GATase1-encoding gene showed the most significant association with bile salt tolerance. In addition, four genes related to DNA damage repair and substance transport were redundant in the strains with high bile salt tolerance. Besides, cluster analysis showed that bile salt hydrolases did not contribute to the bile salt tolerance of L. salivarius. In this study, we determined the global regulatory genes, including LSL_1568, LSL_1716 and LSL_1709, for bile salt tolerance in L. salivarius and provided a potential method for the rapid screening of bile salt-tolerant L. salivarius strains, based on PCR amplification of functional genes.


Introduction
Lactobacillus salivarius was first isolated from children's saliva in Maryland, USA [1]. Later studies have shown that it exists in various habitats, including the oral cavity, the gastrointestinal tract, and in the vagina of animals (including humans), as well as food materials [2][3][4].
In recent years, L. salivarius has drawn attention because of its promising probiotic functions, such as anti-aging, immunoregulation, and bacteriostasis [5]. A study indicated that L. salivarius FDB89, isolated from the feces of centenarians, extended the mean life span of Caenorhabditis elegans by up to 11.9% compared to that of the control group [6]. Furthermore, another L. salivarius strain, FXJCJ7-2, was shown to have significant antiinflammatory effects on lipopolysaccharide (LPS)-treated murine macrophages and mice

Isolation and Genome Sequencing of L. salivarius
The 90 L. salivarius strains used in this study were isolated from 500 stool samples of Chinese populations in 17 provinces. The stool samples were incubated in the modified medium under an anaerobic condition [7]. Each single strain was then picked and conducted a PCR reaction with a 16S rRNA common primer pair (Table S1) [24]. Sanger sequencing of the PCR products was performed by the Genewiz Laboratory (Suzhou, Jiangsu, China). The sequencing data were identified via Blast online (https: //blast.ncbi.nlm.nih.gov/Blast.cgi (accessed on 10 November 2019)).
This study was approved by the Ethics Committee in Jiangnan University, China (SYXK 2012-0002). All donors or their legal guardians provided written informed consent before sample collection. The collection of fecal samples caused no harm or discomfort to the participants. Detailed information about the 90 L. salivarius strains used in this study is listed in Table S2. The isolated strains were preserved at −80 • C in glycerol.

Determination of the Tolerance of L. salivarius Strains in Bile Salt Solutions
Bile salt solutions were prepared according to a published method, modified by changing the final concentration of bile salts to 0.3% [9].
Strains were grown overnight in basic MRS broth at 37 • C for 12 h. After subculturing twice, a 2% (v/v) inoculum of each strain was inoculated into 10 mL of fresh MRS broth and was grown for 12 h. The OD 600 was measured, and the culture of each strain was diluted with MRS to the same absorbance. Cells in the diluted culture of each strain (1 mL) were collected via centrifugation (5000× g, 5 min) and washed twice with sterilized normal saline.
An aliquot (0.5 mL) of each washed cell suspension was used for the original viable count test. Cells in another 0.5 mL aliquot were collected via centrifugation (5000× g, 5 min) and mixed with a 0.3% bile salt solution (0.5 mL). The mixture was then incubated at 37 • C for 4 h. The total number of viable bacteria was determined after the tolerance assay. All experiments were conducted in triplicate. The survival rates of L. salivarius strains in bile salt solutions were calculated according to a previous method [28]. Briefly, the ratio of the total numbers of viable bacteria before and after incubation in the bile salt conditions was calculated for each strain.

Comparative Genomic Analysis
The top eight L. salivarius strains with high survival rates in the bile salt solutions (the tolerant group) and the bottom six strains with low survival rates (the non-tolerant group) were selected for comparative genomic analysis.
The core gene database of the tolerant group was built using OrthoMCL (https:// orthomcl.org/orthomcl/ (accessed on 10 February 2020)) [32] with the default parameters. Next, the whole genome of each non-tolerant strain was compared with the tolerant core database using PGAP (https://github.com/ncbi/pgap (accessed on 15 February 2020)) [33], picking genes deleted in at least two strains via the result of orthologs clusters. Then, copy numbers of the deleting genes in each strain were analyzed using BLAST with its default parameters. Similarly, the redundant genes were determined via comparing the non-tolerant gene database with the tolerance strains' genomes.

Quantitative RT-PCR
The strains with the highest and lowest survival rates were incubated in the bile salt solutions for 4 h after three subcultures, and a parallel culture grown in fresh MRS broth for the same duration was used as a control.
Cells in the culture of each strain (1 mL) were collected via centrifugation (8000× g, 10 min) in an enzyme-free tube and then mixed with 200 µL of lysozyme solution. The mixtures were incubated at 37 • C for 30 min to break the cell walls. Cells without cell walls were collected via centrifugation (12,000× g, 5 min, 4 • C). Total RNA from each strain was extracted using the TRIzol Plus RNA Purification Kit according to the manufacturer's protocol (Invitrogen, Carlsbad, CA, USA). The RNA concentrations were measured using an RNA assay with a NanoDrop spectrophotometer (Thermo, Shanghai, China). cDNA was synthesized in a 20 µL reaction using a HiFiScript gDNA Removal cDNA Synthesis Kit (CoWin Biosciences, Wuxi, Jiangsu, China). qRT-PCR was performed in a 96-well plate following the TransStart Tip Green qPCR SuperMix kit (Tiangen Biotech Co., Ltd., Beijing, China) instruction. The thermocycling program was set as previously described [34], only changing the annealing temperature of the primers. The specific primers were designed using Primer3Plus online (http://www.primer3plus.com/ (accessed on 7 May. 2020)) and were synthesized by the Sangon Biotechnology Laboratory (Shanghai, China) ( Table S1). The 16S rDNA was used as a reference gene [35]. All experiments were conducted in triplicate.

Gene Knockout of L. salivarius FWXBH36M1
A knockout experiment based on the Cre-lox system was performed according to a previously reported method, with minor modifications [36]. Briefly, 500-bp fragments of the upstream and downstream sequences of the target locus identified via RT-PCR were selected to construct gene-specific mutagenesis vectors ( Figure S1a). The designed mutagenesis vectors were synthesized by the Exsyn-Biotechnology Laboratory (Wuxi, Jiangsu, China) and were preserved in Escherichia coli Top10 cells.
E. coli Top10 cells harboring the constructed plasmid were grown aerobically in Luria-Bertani (LB) broth at 37 • C. The plasmids were extracted using the GeneJET Plasmid Miniprep Kit (Thermo, Shanghai, China), following the manufacturer's instructions. The plasmid concentrations were measured using a dsDNA assay with a NanoDrop spectrophotometer (Thermo, Shanghai, China).
Competent L. salivarius FWXBH36M1 cells were prepared according to a previously published method [37]. A 98-µL cell suspension and 2 µL of a plasmid DNA solution (containing approximately 2 µg plasmid) were electroporated in cuvettes with a 1-mm electroporation gap at a voltage of 2.2 kV, using an Electroporator 2510 system (Eppendorf, Shanghai, China). After transformation, 1 mL recovery buffer (MRS containing 1 mM MgCl 2 , 0.1 mM CaCl 2 , and 75.2 mM sorbitol) was added, and the cells were incubated for 4 h at 37 • C. Bacteria were grown on MRS plates containing 20 µg/mL chloramphenicol as a selective marker. The transformants were confirmed via PCR using vector-specific primers (Table S1) and agarose gel electrophoresis. Sanger sequencing of the amplified fragments was performed by the Genewiz Laboratory (Suzhou, Jiangsu, China).
To excise the chloramphenicol gene in the double-crossover mutants, the Cre expression vector with erythromycin-resistant gene was transformed into the mutants. The in-frame deletion mutants were selected as previously described with a minor change [36]. Briefly, after the transformation and recovery, bacteria were grown on MRS plates containing 15 µg/mL erythromycin as a selective marker. Then, the Em r colonies were confirmed via PCR using erythromycin-specific primers (Table S1). The Cre expression vector was cured of Em r colonies by growth without erythromycin selection pressure for 6 generations. Then, the mutants were confirmed via PCR using strain-specific primers (Table S1 and Figure S1b) and agarose gel electrophoresis. Sanger sequencing of the amplified fragments was performed by the Genewiz Laboratory (Suzhou, Jiangsu, China).

Statistical Analysis
Relative gene expression levels were calculated using the 2 −∆∆Ct method [38]. Differences between groups were analyzed via one-way ANOVA and the Kruskal-Wallis test using SPSS Statistics 25 (IBM, Armonk, NY, USA) and were visualized using GraphPad 8 software (GraphPad Software, San Diego, CA, USA). The heatmap package in R-4.0.2 was used to perform a hierarchical cluster analysis and visualize the different functional genes between the two groups.

Different L. salivarius Strains Showed Diverse Survival Rates in Bile Salt Solutions
To compare the tolerance of different L. salivarius strains in bile salt solutions, their survival rates under the condition were analyzed. The results showed differences in intrastrain survival rates (Table S4). Among the 90 L. salivarius strains tested, the survival rates of eight strains exceeded 9.5%, where L. salivarius FWXBH36M1 showed the maximum survival rate (17.76%). The survival rates of six strains were ≤0.7%, and L. salivarius FYNDL5M1 exhibited the lowest survival rate (0.02%).
Notably, this phenotypic variation was significant when the analysis was performed according to the isolated region of each strain ( Figure 1, p < 0.10). Strains from the Guangdong, Gansu, Anhui, and Jiangsu provinces showed higher survival rates. In contrast, strains from the Fujian, Chongqing, Yunnan, and Zhejiang provinces showed lower survival rates. This tendency suggests that the tolerance of strains to bile salt solutions varied with their isolated regions. However, the area-dependent influence disappeared when the strains were classified according to the isolated parts of China ( Figure 1, p > 0.10). trast, strains from the Fujian, Chongqing, Yunnan, and Zhejiang provinces s lower survival rates. This tendency suggests that the tolerance of strains to b tions varied with their isolated regions. However, the area-dependent influe peared when the strains were classified according to the isolated parts of Ch 1, p > 0.10).

Strain Tolerance in Bile Salt Solutions Was Not Associated with BSH Variatio
Previous studies have indicated that BSHs play a critical role in bact tolerance. To determine the contribution of different BSHs to the survival r varius strains in bile salt solutions, the different BSH types per strain were and phylogenetic trees of the three BSH subtypes were constructed.
The results showed that 15 L. salivarius strains had all three BSH subt maining 71 strains had only two types of BSHs (69 with BSHA and BSHB, on and BSHC, and one with BSHB and BSHC). Among them, FWXBH9M2 and had double BSHB-encoding genes. The other four strains only had one typ coding gene (three with BSHA only and one with BSHB only). However, the types of BSHs in each strain did not show region-dependent diversity (Figu Hierarchical clustering indicated that changes in the number of BSHs ge a similar trend as the survival rates of the strains in bile salt solutions (Figu each BSH subtype was further grouped according to its evolutionary dis BSHB, and BSHC had four, seven, and three groups, respectively. Howev analysis combining the BSH groups based on their evolutionary distances an rates of the strains in bile salt solutions did not reveal any similar patterns and Table S4). Hence, the intra-strain diversity of bile salt tolerance in L. saliv related to the variations in BSHs.

Strain Tolerance in Bile Salt Solutions Was Not Associated with BSH Variations
Previous studies have indicated that BSHs play a critical role in bacterial bile salt tolerance. To determine the contribution of different BSHs to the survival rates of L. salivarius strains in bile salt solutions, the different BSH types per strain were determined, and phylogenetic trees of the three BSH subtypes were constructed.
The results showed that 15 L. salivarius strains had all three BSH subtypes. The remaining 71 strains had only two types of BSHs (69 with BSHA and BSHB, one with BSHA and BSHC, and one with BSHB and BSHC). Among them, FWXBH9M2 and FSDHZD3L5 had double BSHB-encoding genes. The other four strains only had one type of BSHencoding gene (three with BSHA only and one with BSHB only). However, the numbers and types of BSHs in each strain did not show region-dependent diversity (Figure 2a).
Hierarchical clustering indicated that changes in the number of BSHs genes followed a similar trend as the survival rates of the strains in bile salt solutions ( Figure 2a). Then, each BSH subtype was further grouped according to its evolutionary distance. BSHA, BSHB, and BSHC had four, seven, and three groups, respectively. However, a further analysis combining the BSH groups based on their evolutionary distances and the survival rates of the strains in bile salt solutions did not reveal any similar patterns (Figure 2b-d and  Table S4). Hence, the intra-strain diversity of bile salt tolerance in L. salivarius was not related to the variations in BSHs. The phylogenetic trees of the strains' BSHA, BSHB, and BSHC, respectively, using L. johnsonii PF01 as an outgroup. Strain names in red represent strains with survival rates higher than 9.5%; strain names in green represent strains with survival rates lower than 0.7%.

Comparative Genomic Analysis Determined the potential L. salivarius Genes Responsible for Bile Salt Tolerance
To determine the potential mechanisms of how functional genes were involved in the bile salt tolerance of L. salivarius, a comparative genomic analysis of eight tolerant (survival rate ≥ 9.5%) and six non-tolerant (survival rate ≤ 0.7%) strains was performed.
The results showed that there were 15 tolerant core genes deleting in at least two nontolerant strains (Figure 3a). Among them, four were involved in carbohydrate transport, three encoded for different glycoside hydrolases, and six (two for each following function) were responsible for glycosyltransferase, transcription regulation, and hypothetical protein-encoding, respectively. The remaining two genes contributed to amidotransferase and phosphorylase encoding, respectively ( Table 2). Further analysis via CD-BLAST revealed that the gene encoding amidotransferase (LSL_1709) was highly homologous with the Ydr533c protein-encoding genes in Saccharomyces cerevisiae, which are upregulated in response to various stress conditions [39].
Notably, further analysis of the gene functions showed that two functional clusters existed in the different genes, including four genes involved in peptidoglycan synthesis and seven genes in the phosphotransferase system (PTS) (Figure 3b). In the former cluster, three glycosyltransferase-related genes were involved in the acylation of N-acetylmuramic acid and N-acetylglucosamine, and one gene encoded β-N-acetylglucosaminidase, participating in the hydrolysis of chitin (Table 2). KEGG analysis identified that the latter cluster was enriched in mannose-and N-acetylglucosamine-specific PTS pathways. In particular, two genes (LSL_1716 and LSL_1715) encoding the PTS system transporter subunits IIB and IIC, respectively, were involved in mannose transmembrane transport. Two genes (LSL_1714 and LSL_1713), which encode for the PTS system transporter subunits IIABC and IIBC, respectively, are involved in N-acetylglucosamine transmembrane transportation (Table 2, Figure S3). The phylogenetic trees of the strains' BSHA, BSHB, and BSHC, respectively, using L. johnsonii PF01 as an outgroup. Strain names in red represent strains with survival rates higher than 9.5%; strain names in green represent strains with survival rates lower than 0.7%.

Comparative Genomic Analysis Determined the Potential L. salivarius Genes Responsible for Bile Salt Tolerance
To determine the potential mechanisms of how functional genes were involved in the bile salt tolerance of L. salivarius, a comparative genomic analysis of eight tolerant (survival rate ≥ 9.5%) and six non-tolerant (survival rate ≤ 0.7%) strains was performed.
The results showed that there were 15 tolerant core genes deleting in at least two nontolerant strains (Figure 3a). Among them, four were involved in carbohydrate transport, three encoded for different glycoside hydrolases, and six (two for each following function) were responsible for glycosyltransferase, transcription regulation, and hypothetical proteinencoding, respectively. The remaining two genes contributed to amidotransferase and phosphorylase encoding, respectively ( Table 2). Further analysis via CD-BLAST revealed that the gene encoding amidotransferase (LSL_1709) was highly homologous with the Ydr533c protein-encoding genes in Saccharomyces cerevisiae, which are upregulated in response to various stress conditions [39].
Notably, further analysis of the gene functions showed that two functional clusters existed in the different genes, including four genes involved in peptidoglycan synthesis and seven genes in the phosphotransferase system (PTS) (Figure 3b). In the former cluster, three glycosyltransferase-related genes were involved in the acylation of N-acetylmuramic acid and N-acetylglucosamine, and one gene encoded β-N-acetylglucosaminidase, participating in the hydrolysis of chitin (Table 2). KEGG analysis identified that the latter cluster was enriched in mannose-and N-acetylglucosamine-specific PTS pathways. In particular, two genes (LSL_1716 and LSL_1715) encoding the PTS system transporter subunits IIB and IIC, respectively, were involved in mannose transmembrane transport. Two genes (LSL_1714 and LSL_1713), which encode for the PTS system transporter subunits IIABC and IIBC, respectively, are involved in N-acetylglucosamine transmembrane transportation (Table 2, Figure S3). Microorganisms 2021, 9, x FOR PEER REVIEW 9 of 16

Possible Redundant L. salivarius Genes for Bile Salt Tolerance Were Identified via Comparative Genomic Analysis
Next, the whole genomes of the tolerant strains were compared with the core genes of the non-tolerant group using PGAP to determine the redundant genes under bile salt solutions. The results showed seven different genes between the two groups, including one TetR/AcrR family transcriptional regulator-encoding gene, one site-specific integraseencoding gene, one aspartate carbamoyltransferase-encoding gene, and four genes that  Next, the whole genomes of the tolerant strains were compared with the core genes of the non-tolerant group using PGAP to determine the redundant genes under bile salt solutions. The results showed seven different genes between the two groups, including one TetR/AcrR family transcriptional regulator-encoding gene, one site-specific integraseencoding gene, one aspartate carbamoyltransferase-encoding gene, and four genes that encode hypothetical proteins ( Table 2). The copy numbers of all seven genes were reduced in the bile salt-tolerant strains (Figure 3a).
After combining the annotations of the upstream and downstream genes of the unknown genes, it could be inferred that LSL_0261, LSL_0166, and LSL_0167 were related to substance metabolism and transportation ( Figure S4). In addition, CD-BLAST showed that the LSL_0252-encoding protein belongs to the HTH_17 superfamily, involved in DNA recombination.

qRT-PCR Showed the Transcriptional Changes of the Possible Functional and Redundant Genes Related to Tolerance in Bile Salt Solutions
qRT-PCR was conducted on bile salt-tolerant (FWXBH36M1 and FJLHD9M1) and nontolerant strains (FYNDL5M1 and FCQHC3ML6) to determine the transcriptional changes of the functional and redundant genes identified above upon culture in bile salt solutions. According to their genome coordinates and functions, the tolerance-related genes could be classified into five groups (Table S3), and one target gene from each group was selected for qRT-PCR. qRT-PCR for the redundant genes was performed after excluding genes with high copy number changes in both groups (pyrB and xerC) (Figure 3a).
For the tolerance-related genes, the results showed that the relative expression level of LSL_1568, LSL_1716, LSL_1709, and LSL_0951 were increased in both strains when strains were grown in bile salt solutions (Figure 4a,b). However, the relative expression level of rfaG showed a slight decrease in FWXBH36M1 (0.18 on average); in FJLHD9M1, it showed a maximum increase (41.12 on average).

Knockout of the Functional Genes in L. salivarius FWXBH36M1 Verified Their Contributions to the Bile Salt Tolerance
To prove that the functional genes obtained above influenced the tolerance of L. salivarius to bile salt solutions, in-frame deletions of the top three most highly relatively expressed genes in FWXBH36M1 were performed. Confirmation of positive transformants via DNA electrophoresis is shown in Figure S5.
Compared with that of the wild-type strain, the growth profile of the three mutant strains in MRS did not show any significant differences (Figure 5a, p > 0.05), indicating that gene knockout did not affect the growth of the strains under normal conditions.
However, the survival rates of the three mutant strains were significantly lower than that of the wild-type strain under bile salt conditions (Figure 5b, p < 0.05). Additionally, the survival rates of the three mutant strains showed a consistent trend with the relative expression of the corresponding genes. Specifically, the knockout of the most highly expressed gene in FWXBH36M1 also led to the most significant drop-in survival rate.
Therefore, the functional genes obtained through the comparative genomic analysis were related to the survival rates of L. salivarius strains in bile salt solutions. The four genes with increased relative expression levels were changed to different extents in the two tested strains. In FWXBH36M1, the relative expression of LSL_1568, LSL_1716, and LSL_1709 showed remarkable changes compared with that of LSL_0951. There was an opposite trend in the relative expression level of these functional genes in FJLHD9M1, where LSL_0951 demonstrated a higher expression level than others. The differences in the fold-changes of the functional genes upon growth in a bile salt system may explain the intra-strain diversity in the bile salt tolerance of L. salivarius.
In addition, an inconsistency in the relative expression of the redundant genes was found in strains FYNDL5M1 and FCQHC3ML6, despite the noticeable growth of all four genes (Figure 4c,d). This inconsistency may have resulted in the fluctuations in the survival rates in the non-tolerant strains.

Knockout of the Functional Genes in L. salivarius FWXBH36M1 Verified Their Contributions to the Bile Salt Tolerance
To prove that the functional genes obtained above influenced the tolerance of L. salivarius to bile salt solutions, in-frame deletions of the top three most highly relatively expressed genes in FWXBH36M1 were performed. Confirmation of positive transformants via DNA electrophoresis is shown in Figure S5.
Compared with that of the wild-type strain, the growth profile of the three mutant strains in MRS did not show any significant differences (Figure 5a, p > 0.05), indicating that gene knockout did not affect the growth of the strains under normal conditions.

Discussion
In this study, we identified the potential functional genes that confer bile salt tolerance in L. salivarius via comparative genomic analysis and verified the results via qRT-PCR and gene knockout experiments.
Ninety L. salivarius strains showed intra-strain differences in survival rates in bile salt solutions, indicating that intrinsic bile salt tolerance is strain-specific. Similar results have been reported in previous studies on Bifidobacterium [40] and L. acidophilus strains [41]. Considering the open pan-genome of L. salivarius [42], phenotypic differences might have resulted from gene variations. The region-level analysis demonstrated that the phenotypic differences were area-dependent, which might be attributed to the niche-specific adaptive evolution of Lactobacillus [43]. Notably, as the grouping area expanded further, this influence disappeared. Similar results have been reported for L. plantarum [44], thus suggesting that the grouping level of strain habitats greatly influences the results of genetic evolution analysis.
Further analysis of BSH variations did not reveal any associations between the BSHs and the bile salt resistance of L. salivarius strains, thus indicating that intra-strain phenotypic differences in bile salt tolerance were induced by mechanisms other than the presence of the BSHs. Similar results have been reported in previous studies on L. johnsonii However, the survival rates of the three mutant strains were significantly lower than that of the wild-type strain under bile salt conditions (Figure 5b, p < 0.05). Additionally, the survival rates of the three mutant strains showed a consistent trend with the relative expression of the corresponding genes. Specifically, the knockout of the most highly expressed gene in FWXBH36M1 also led to the most significant drop-in survival rate.
Therefore, the functional genes obtained through the comparative genomic analysis were related to the survival rates of L. salivarius strains in bile salt solutions.

Discussion
In this study, we identified the potential functional genes that confer bile salt tolerance in L. salivarius via comparative genomic analysis and verified the results via qRT-PCR and gene knockout experiments.
Ninety L. salivarius strains showed intra-strain differences in survival rates in bile salt solutions, indicating that intrinsic bile salt tolerance is strain-specific. Similar results have been reported in previous studies on Bifidobacterium [40] and L. acidophilus strains [41]. Considering the open pan-genome of L. salivarius [42], phenotypic differences might have resulted from gene variations. The region-level analysis demonstrated that the phenotypic differences were area-dependent, which might be attributed to the niche-specific adaptive evolution of Lactobacillus [43]. Notably, as the grouping area expanded further, this influence disappeared. Similar results have been reported for L. plantarum [44], thus suggesting that the grouping level of strain habitats greatly influences the results of genetic evolution analysis.
Further analysis of BSH variations did not reveal any associations between the BSHs and the bile salt resistance of L. salivarius strains, thus indicating that intra-strain phenotypic differences in bile salt tolerance were induced by mechanisms other than the presence of the BSHs. Similar results have been reported in previous studies on L. johnsonii [18] and L. salivarius [6], indicating that bacterial bile salt resistance is a multifactorial phenomenon [10].
Comparative genomic analysis based on the whole bacterial genome showed that the bile salt resistance of L. salivarius strains was associated with the peptidoglycan synthesis process. Peptidoglycan is the main component of the cell wall [45]. The cell wall is the first line of defense for microorganisms against the external environment [10]. Therefore, strains with a higher abundance of genes related to cell wall synthesis potentially have stronger survival abilities in bile salt solutions. Several studies on L. fermentum [19], L. johnsonii [18], and L. plantarum [11] found similar results regarding the transcription level of these genes under bile salt conditions.
Comparative genomic analysis and KEGG enrichment analysis indicated that bacterial PTS also played an important role in the bile salt resistance of L. salivarius strains. PTS is a complex kinase system that catalyzes sugar transport and phosphorylation [46]. It is composed of the enzyme I, enzyme II, and HPr (heat-stable). Enzyme II is further composed of the A, B, and C subunits. Our results showed that the bile salt-tolerant strains had a mannose-specific translocator-encoding gene. Notably, mannose is inefficient in terms of energy supply, so most bacteria do not use this sugar for energy supplementation [47]. In a stressful environment, however, bacteria tend to extensively use various carbon sources in the environment to meet their survival needs [20]. Therefore, strains enriched with energy intake genes have a potential survival advantage.
Additionally, the relative expression of LSL_1709 in the strain with maximum stress tolerance increased dramatically in the bile salt condition, and the survival rate of LSL_1709mutant strain significantly lowered in the bile salt solutions. All of these indicated this gene has a strong association with the strains' bile salt-tolerant ability. Considering the chaperone activity of LSL_1709, its encoding product may involve the strains' bile salt tolerance by inhibiting the misfolding and denaturation of protein caused by the stress condition [48].
Analysis of the strains' redundant genes related to bile salt tolerance showed that four potential redundant genes were related to the repair of DNA damage, substance transportation and transcription regulation. All of these can involve repairing bacterial damage caused by bile salts [49,50]. This may explain the up-regulated expression of the four genes in the bile salt condition. However, for the strains with high bile salt tolerance, the four genes exhibited significant variations, which implies that these genes are involved in the downstream pathways of bile resistance regulation. Also, considering the redundant metabolism-related genes were involved in the amino acid usage, which is a priority energy utilization strategy. Therefore, all these genes were not necessary for tolerant strains to keep these genes.
Interestingly, the relative expressions of the possible functional and redundant genes varied from strain to strain, indicating that even the same gene has different responses to bile salt conditions in different strain's genomic backgrounds. This phenomenon can also be observed in S. cerevisiae [51] and Bacillus mycoides [52], which may be the reason for the huge differences in the phenotypes of strains. Therefore, in future functional genome studies, we should fully consider the contribution weight of different functional genes to the specific phenotype in different bacterial genome backgrounds and develop an appropriate phenotypic prediction scheme based on this, which will effectively improve the accuracy of the results.
However, it should be noticed that this research was based on in vitro experiments only, lacking the corresponding in vivo confirmation. Previous studies have shown that the maximum concentrations of bile salts tolerated by probiotics differ in vivo and in vitro. Besides, the types of bile salts used and the bile conditions in the in vitro experiments were different from those in the human intestine [53]. Therefore, it is necessary to assess the strain's tolerant ability in the host condition where they are used to precisely uncover the host-specific tolerant mechanisms of probiotic strains.
Overall, in this study, we applied comparative genomic analyses to determine the mechanism of how the gene variations involved in the bile salt tolerance of L. salivarius. The results indicated that the bile salt resistance of L. salivarius strains was mainly related to the chaperones, PTS and peptidoglycan synthesis. Our work also provided a potential method for rapidly screening the bile-salt-tolerant strains of L. salivarius based on PCR of the functional genes.