Antimicrobial Susceptibility of Lactobacillus delbrueckii subsp. lactis from Milk Products and Other Habitats

As components of many cheese starter cultures, strains of Lactobacillus delbrueckii subsp. lactis (LDL) must be tested for their antimicrobial susceptibility to avoid the potential horizontal transfer of antibiotic resistance (ABR) determinants in the human body or in the environment. To this end, a phenotypic test, as well as a screening for antibiotic resistance genes (ARGs) in genome sequences, is commonly performed. Historically, microbiological cutoffs (MCs), which are used to classify strains as either ‘sensitive’ or ‘resistant’ based on the minimal inhibitory concentrations (MICs) of a range of clinically-relevant antibiotics, have been defined for the whole group of the obligate homofermentative lactobacilli, which includes LDL among many other species. This often leads to inaccuracies in the appreciation of the ABR status of tested LDL strains and to false positive results. To define more accurate MCs for LDL, we analyzed the MIC profiles of strains originating from various habitats by using the broth microdilution method. These strains’ genomes were sequenced and used to complement our analysis involving a search for ARGs, as well as to assess the phylogenetic proximity between strains. Of LDL strains, 52.1% displayed MICs that were higher than the defined MCs for kanamycin, 9.9% for chloramphenicol, and 5.6% for tetracycline, but no ARG was conclusively detected. On the other hand, all strains displayed MICs below the defined MCs for ampicillin, gentamycin, erythromycin, and clindamycin. Considering our results, we propose the adaptation of the MCs for six of the tested clinically-relevant antibiotics to improve the accuracy of phenotypic antibiotic testing.


Introduction
The thermophilic lactic acid bacterium (LAB) Lactobacillus delbrueckii is important in many traditional fermented foods prepared worldwide (Table 1), and three subspecies (bulgaricus, delbrueckii, and lactis) have a long history of safe use [1]. To date, six subspecies have been described: L. delbrueckii subsp. bulgaricus (LDB), L. delbrueckii subsp. delbrueckii (LDD), L. delbrueckii subsp. indicus (LDI), L. delbrueckii subsp. jakobsenii (LDJ), L. delbrueckii subsp. lactis (LDL), and L. delbrueckii subsp. sunkii (LDS). In dairy fermentations, LDB is mainly used in yogurt making, whereas LDL is traditionally used in the production of cooked cheeses, owing its tolerance to the high temperatures during the early phases of cheese manufacturing. Both subspecies seem to have developed a similar adaptation towards optimized utilization of milk resources through reductive evolution and limited acquisition of particular functions. In its evolutionary path, LDB lost more functions than LDL and has thus further diverged from their common ancestor [2].

Bacterial Strains and Culture Conditions
A total of 101 strains of L. delbrueckii were selected for analysis ( Table 2). Identification using matrix-assisted laser desorption-ionization time of flight (MALDI-TOF) and genome comparison based on the calculated average nucleotide identity (ANI) values (see Identification using MALDI-TOF and Average Nucleotide Identities) revealed inaccuracies in the taxonomic assignments. One strain assigned to another species was excluded from this study. After being re-assigned at the subspecies level, the remaining 100 L. delbrueckii strains (80 LDL, 17 LDB, and three LDS) were further analyzed. Seventy-eight strains obtained from various dairy products or dairy starter cultures, two from distilleries, one from human urine, one from human saliva, one from dried calf stomachs, one from fermented vegetables (sunki), and 16 others were of unknown origin. Six strains were constituents of commercial starter cultures used in the cheese industry. The strains were stored at −80 • C in De Man, Rogosa, Sharpe (MRS) [22] broth with Tween 80 (Biolife Italiana Srl, Milan, Italy) containing sterile low-fat milk as the cryoprotectant. They were routinely cultured in MRS broth with Tween 80 at 37 • C under aerobic conditions. The absence of contaminants was confirmed by plating 10-fold serial dilutions in 10 mL 0.9% NaCl onto MRS agar plates and by observing the colony morphologies and the cells under a microscope.

Genome Sequences of the L. delbrueckii Strains
Available assemblies were downloaded from the National Center for Biotechnology Information (NCBI). Strains with unavailable assemblies in open databases were sequenced de novo. Their DNA was extracted as described previously [24]. Quality control assessment of the extracted DNA, library generation, and sequencing runs were performed on the Next Generation Sequencing Platform, University of Bern, Switzerland. In brief, the libraries were prepared using a TruSeq DNA PCR-free  [25]. Quality control of raw data was performed with FastQC v.0.11.7 [26]. Adaptor removal and trimming of raw data were performed using fastp v.0.20.0 [27]. The trimmed reads were assembled with SPAdes v.3.14.0 [28] in the -isolate mode, and contigs shorter than 200 bp were removed from the final assembly. QUAST v.4.6 [29] and BUSCO v.4.0.6 [30] were used in the -auto-lineage-prok mode to assess the quality of the assemblies. Newly obtained assemblies were uploaded to NCBI. Taxonomic affiliation of the assemblies was performed at the species level using the GTDB-Tk v.0.3.2 'classify' workflow and reference data version r89 to confirm the previous identification. The Whole Genome Shotgun projects of the strains sequenced in this study have been deposited at GenBank under BioProject PRJNA777018 with GenBank WGS accessions JAJNSH000000000-JAJNVX000000000. The version described in this paper is version 01.

Strains' Phylogeny
Phylogenetic relationships of the isolates were assessed using PhaME v.1.0.2 [35] and IQ-TREE multicore v.2.0.3 [36], as described previously [19]. LDL and LDB were assessed separately. The resulting trees were visualized with the Interactive Tree Of Life (iTOL) v.5 [37]. Pairs of isolates with fewer than 50 single nucleotide polymorphisms (SNPs) in the SNP pairwise matrix from PhaME were considered to be possibly identical. In such cases, their origins and years of isolation were taken as an additional information to determine whether they should be considered different 'strains in the taxonomic sense' [38].

Antibiotic Susceptibility Testing and MC Determination
The existing MCs were evaluated based on the Standard Operating Procedure (SOP) 10.1 of the EUCAST [45]. The antibiotic susceptibility of the strains was tested at the Culture Collection of Switzerland (CCOS, Wädenswil, Switzerland) using a broth microdilution susceptibility method following the standard procedure of the International Organization for Standardization ISO 10932|IDF 223:2010 [46], as recommended by the FEEDAP [17].

Taxonomic Assignment and ANIs
At the species level, the identity of most strains was confirmed using MALDI-TOF and GTDB-tk. Only strain CIRM BIA 879 was assigned to L. johnsonii and was excluded from the subsequent analyses.
Based on the ANI values, LDB formed a well-defined cluster, with pairs of strains displaying ANI similarities at 98.83-100% (median 99.19%), and it was evidently distinct from the other L. delbrueckii subspecies ( Figure S1). By contrast, the LDL strains formed a quite inhomogeneous cluster (ANI 97.59-100%, median 98.79%), although it was still evidently divergent from the other L. delbrueckii subspecies.
The taxonomic assignment at the subspecies level was re-evaluated for three strains (

SNP-Based Phylogenetic Relationships of the Strains
The SNP-based analysis revealed that the LDB strains were all phylogenetically distinct ( Figure 1). The closest pairs of strains displayed 34 (CIRM BIA 1379 and CIRM BIA 1581) and 94 (CIRM BIA 864 and CIRM BIA 1623) SNPs. As the origins of these isolates were different, they were considered closely related but different strains. dently divergent from the other L. delbrueckii subspecies.
The taxonomic assignment at the subspecies level was re-evaluated for three strains ( Table 2): CIP 101810 and CIRM BIA 266 (both originally assigned to LDL) displayed a high similarity (98.27% and 98.61%, respectively) with the LDS type strain DSM 24966 T , and CIRM BIA 1375 (originally LDB) displayed similarity with LDL DSM 20072 T (98.53%). NCIMB 702468 and NCIMB 702469, referenced in the NCIMB database as LDS and LDD, respectively, were assigned to LDL.

SNP-Based Phylogenetic Relationships of the Strains
The SNP-based analysis revealed that the LDB strains were all phylogenetically distinct ( Figure 1). The closest pairs of strains displayed 34 (CIRM BIA 1379 and CIRM BIA 1581) and 94 (CIRM BIA 864 and CIRM BIA 1623) SNPs. As the origins of these isolates were different, they were considered closely related but different strains. A more complex pattern was obtained for the LDL strains ( Figure 2). Twenty-eight strains formed a homogeneous phylogenetic group. All pairs of strains in this group, however, displayed more than 120 SNPs (median 1067 SNPs). All but one (CIRM BIA 1375) were isolated from Swiss dairy products, mainly from cheese or from undefined mixed starter cultures, mostly obtained from different cheese factories. Several other groups consisted of similar, though different, strains: FAM 1200, FAM 22091, FAM 22092, FAM 22093, FAM 22274, and FAM 22680 displayed more than 368 SNPs. These strains were isolated from Swiss whey (no information available for FAM 1200). NCIMB 700820 and NCIMB 700860 displayed more than 2000 SNPs. Both CIRM BIA 230 and CIRM BIA 265 were isolated from Finnish Emmental cheese in 1968, but they were considered different strains, as they displayed 138 SNPs. CIRM BIA 233 and CIRM BIA 234 were both isolated from a French yogurt factory in 1984, and they displayed 12 SNPs. Given the strong suspicion that they belong to the same strain, CIRM BIA 234 was not included in the antibiotic susceptibility profiles. CIRM BIA 269 displayed 12 SNPs when compared with CIRM BIA 233 and CIRM BIA 234, but it was still regarded as a different strain given its different origin. A more complex pattern was obtained for the LDL strains ( Figure 2). Twenty-eight strains formed a homogeneous phylogenetic group. All pairs of strains in this group, however, displayed more than 120 SNPs (median 1067 SNPs). All but one (CIRM BIA 1375) were isolated from Swiss dairy products, mainly from cheese or from undefined mixed starter cultures, mostly obtained from different cheese factories. Several other groups consisted of similar, though different, strains: FAM 1200, FAM 22091, FAM 22092, FAM 22093, FAM 22274, and FAM 22680 displayed more than 368 SNPs. These strains were isolated from Swiss whey (no information available for FAM 1200). NCIMB 700820 and NCIMB 700860 displayed more than 2000 SNPs. Both CIRM BIA 230 and CIRM BIA 265 were isolated from Finnish Emmental cheese in 1968, but they were considered different strains, as they displayed 138 SNPs. CIRM BIA 233 and CIRM BIA 234 were both isolated from a French yogurt factory in 1984, and they displayed 12 SNPs. Given the strong suspicion that they belong to the same strain, CIRM BIA 234 was not included in the antibiotic susceptibility profiles. CIRM BIA 269 displayed 12 SNPs when compared with CIRM BIA 233 and CIRM BIA 234, but it was still regarded as a different strain given its different origin. CIRM BIA 267 had the same origin as CIRM BIA 269, but these isolates displayed more than 150 SNPs. NCIMB 8183, NCIMB 8964, DSM 20355, DSM 20073, and DSM 20076 all displayed more than 100 SNPs. NCIMB 7278 and NCIMB 701437 displayed 36 SNPs. As information about their origin was unavailable, they were considered identical, and only NCIMB 701437 was included in the antibiotic susceptibility profiles. NCIMB 8011 and NCIMB 8170, although considerably closely related to NCIMB 701437, were still considered different strains, as they displayed more than 50 SNPs with respect to the latter strain and between each other. Finally, FAM 24849, FAM 24850, and FAM 24852 all displayed fewer than 12 SNPs. As these strains were isolated from different commercial starter cultures but from the same company, there is a high probability that the same strain was used in all three starter cultures. Therefore, only FAM 24850 was included in the susceptibility profiles.
NCIMB 8170, although considerably closely related to NCIMB 701437, were still considered different strains, as they displayed more than 50 SNPs with respect to the latter strain and between each other. Finally, FAM 24849, FAM 24850, and FAM 24852 all displayed fewer than 12 SNPs. As these strains were isolated from different commercial starter cultures but from the same company, there is a high probability that the same strain was used in all three starter cultures. Therefore, only FAM 24850 was included in the susceptibility profiles.

ARGs in the Genome Assemblies
The search for ARGs in the strain assemblies revealed a single match, namely strain FAM 22754. In the assembly of this strain, an ermB gene was detected in a small, low coverage, scaffold (763 bp). To confirm/infirm the presence of this gene in the strain FAM 22754, we performed a specific PCR targeting ermB, as described previously [50]. The presence of ermB could neither be confirmed using the DNA extracts used for sequencing, nor using new DNA extracts from the strain, and we concluded that its presence in the assembly was due to a contamination during DNA extraction, library preparation, or sequencing.  Figure 3). The MCs for all the investigated subspecies were higher than the FEEDAP MCs for KAN (256 mg/L instead of 16 mg/L), TET (16 mg/L instead of 4 mg/L), and CHL (8 mg/L instead of 4 mg/L).   As for LDL (Table 3, Figure 4), 37 strains out of 71 (52.1%) displayed MICs higher than the FEEDAP MC for KAN, seven (9.9%) for CHL, four (5.6%) for TET, and one (1.4%) each for VAN      No link could be drawn between the levels of resistance to antibiotics and the origins of the strains or their position in the SNP-based phylogenetic tree (Figure 2).

Discussion
Used alongside S. salivarius subsp. thermophilus, LDL is the major L. delbrueckii subspecies added as a thermophilic starter in the production of cooked Swiss cheese varieties. The development of starter or adjunct cultures for the cheese market involves a thorough screening of strain characteristics. MIC values determined through broth microdilution testing are commonly used as an easy criterion for the early acceptance or rejection of candidate strains prior to the in-depth characterization of their properties to determine their suitability for cheese production. In the case of LDL, the ABR status of new candidates is often subject to debate, as values regularly exceed the defined FEEDAP MCs, particularly that for KAN (Agroscope, data not shown). The FEEDAP MCs have been defined for the entire group of obligate homofermentative lactobacilli, consisting of at least 31 species [18]. Here, we questioned the relevance of the current FEEDAP MCs for L. delbrueckii, more specifically for LDL.
An unambiguous identification of a strain is a prerequisite for any further safety assessment [17]. Controlling the taxonomic affiliation of our strain assortment led to the re-assignment of one and five strains at the species and subspecies levels, respectively. In addition, during the preliminary selection of strains suitable for this study, six strains in the Agroscope Culture Collection were taxonomically re-assigned at the subspecies level (data not shown). Here, we assessed the taxonomic affiliation at the subspecies level based on the ANI values. According to previous research, strains whose genomes display an ANI of ≥96.5% and an alignment fraction of ≥0.6 can be grouped into a single species [51]. At the subspecies level, however, no similar cutoffs have been defined. We therefore assessed visually the positions of the strains and their similarities in the ANI-based tree relative to different type strains. The subspecies of strains in some culture collections have often been defined based solely on physiological properties, and modern molecular techniques have not been systematically applied. Today, as NGS has become an affordable approach, taxonomic affiliations can be reassessed in the light of whole genome comparisons. Considering the taxonomic assessment performed here and the subsequent re-classifications, we strongly recommend that the taxonomic affiliation of the examined strains, even those deposited in culture collections, be systematically questioned. Comparisons at the genomic level are a good approach; however, more rapid methods that require less computational analysis (e.g., MALDI-TOF) may also provide a good evaluation of certain taxa. In our study, MALDI-TOF allowed for the discrimination of LDL and LDB, but not of LDS, as this subspecies is not yet referenced in the BDAL v11.0 library used here.
Most strains selected for analysis in our study were dairy isolates (obtained from milk, dairy products, natural whey, and other starter cultures) from Switzerland, as most strains were provided by the Agroscope Culture Collection, which includes strains mainly obtained from Swiss dairy products. To determine the MICs of the non-dairy isolates, we broadened the range of origins by including strains obtained from other habitats, such as human saliva, calf stomachs, or distilleries. Moreover, care was taken in expanding the strain selection to include those obtained from heterogeneous dairy products and those with different geographical origins. Furthermore, the selected strains were isolated within a period covering at least 69 years (1950-2019), within which the golden age of antibiotics and the subsequent dissemination of ABR determinants is included [52]. On the one hand, our strain selection allowed to avoid potential habitat-related biases. Although strains in some environments may be more exposed to selective pressure exerted by antibiotics or other antimicrobials, the spread of ABR is less likely to occur between habitats than within a single habitat [53,54]. Therefore, the overall elevated resistance towards KAN observed here is most likely an intrinsic characteristic of L. delbrueckii rather than a resistance that has been transmitted across habitats and time. On the other hand, selecting strains of different origins increases the chance of getting as many different strains as possible.
MCs were calculated using MIC distributions based on the susceptibility profiles of distinct strains. Here, potential duplicate strains were discarded after empirically examining the strains' relationships in a SNP-based phylogenetic tree and after determining the number of SNPs between pairs of strains. This step is crucial to avoid duplicated MIC results that would lead to biased MIC distributions. The number of SNPs below which two strains were considered identical was defined subjectively based on previous SNPs analyses of different strains (data not shown) and may be adapted. Although the use of SNPs alone to discriminate between different strains is not standard, the threshold defined in this study is relatively consistent with results of other studies investigating clinical isolates [55][56][57]. Furthermore, the defined threshold of 50 SNPs appears reasonable considering the origins of the analyzed strains.
No association between the origins of the strains and their antimicrobial susceptibility patterns could be evidenced, suggesting that no particular antibiotic susceptibilities were selected in various ecological niches. Similarly, no link could be drawn between the year of isolation and the antimicrobial susceptibility patterns, suggesting that the selection pressure exerted by antibiotics used in veterinary medicine and farming towards L. delbrueckii has remained relatively moderate.
In this study, we postulated that the observed phenotypic resistances in LDL are actually false positives due to the inaccurate MCs. Our results tend to confirm our hypothesis. Indeed, although more than half of the tested LDL strains were phenotypically defined as 'resistant' to KAN according to the FEEDAP criteria, none of them displayed any known ARG. Similar observations were previously reported, with, e.g., 25% of tested LDB having MICs of KAN above the EFSA MC but no detected ARG [58]. The proportion of L. delbrueckii strains with MICs above the MCs for chloramphenicol (8.9%) in our study is higher than previously reported [58,59]. Conversely, the proportion of TET-resistant strains was lower in our study (4.4%) than in a previous study carried out on 11 L. delbrueckii strains [59]. The range of MICs of TET seems to be relatively broad by L. delbrueckii, as shown here, but also considering results from other studies. Indeed, the proportion of TET-resistant strains found here was lower (4.4%) than in a previous study carried out on 11 L. delbrueckii strains [59] but higher than in another one based on 4 LDL strains from natural whey starters, where no resistance to GEN, ERY, TET, and VAN was detected [60]. In addition to the fact that we detected no ARGs, all MIC distributions of LDL were unimodal, while a bimodal distribution would be expected in case of an acquired antibiotic resistance. On the other hand, the high susceptibilities to AMP, ERY, and CLI observed in most of our investigated strains represent a common feature of Lactobacillus delbrueckii that has already been reported by others [58,59,61]. For these three antibiotics, it rather seems that the MCs have been defined too high, which may open the door to false negative results and consequently the potential introduction of antibiotic resistant strains into the food chain.
Based on our observations, we therefore conclude that the MCs defined for the entire group of homofermentative lactobacilli should be redefined at the species or even at the subspecies level. Our strain selection may be too narrow in order that new MCs could be proposed for the entire species, as we have not tested some of the subspecies for their antimicrobial susceptibilities (LDD, LDI, and LDJ). However, our calculations based on strains of the subspecies bulgaricus, lactis, and sunkii suggest that the MCs for KAN, TET and CHL for L. delbrueckii should be increased, whereas the MCs for AMP, ERY and CLI should be decreased. For LDL, we suggest that the MCs for KAN, TET and CHL be increased to at least 128, 16, and 16 mg/L, respectively. Conversely, the MCs for AMP, ERY and CLI could be reduced to 0.25, 0.125, and 0.063 mg/L, respectively.
Recently, Stefanska et al. detected 32% of KAN-resistant strains in a wider range of LAB species using broth microdilution [62]. Of these strains, only three (4.6%) carried aph(3 )-IIIa, a gene conferring resistance to this antibiotic. The issue raised in our study and that has already been reported for other species may thus be widespread in other LAB species, too [19,23].
The facts that a non-negligible number of strains were unable to grow in the test medium and that some variations occur in the phenotypic measures support the dual approach proposed by the FEEDAP, that is, phenotypic testing complemented by a search of the WGS for known ARGs [17]. The phenotypic approach may allow for the detection of yet unknown ABR mechanisms and is relatively advantageous in terms of cost and time. It can thus be used as a primary screening approach, along with the use of accurate MCs to avoid false positives that may drastically reduce the number of candidates [19]. Meanwhile, the WGS approach provides an appropriate method to evaluate strains that do not grow on test media and to confirm the resistances observed using the phenotypic test.