Genetic Diversity of Leuconostoc mesenteroides Isolates from Traditional Montenegrin Brine Cheese

In many dairy products, Leuconostoc spp. is a natural part of non-starter lactic acid bacteria (NSLAB) accounting for flavor development. However, data on the genomic diversity of Leuconostoc spp. isolates obtained from cheese are still scarce. The focus of this study was the genomic characterization of Leuconostoc spp. obtained from different traditional Montenegrin brine cheeses with the aim to explore their diversity and provide genetic information as a basis for the selection of strains for future cheese production. In 2019, sixteen Leuconostoc spp. isolates were obtained from white brine cheeses from nine different producers located in three municipalities in the northern region of Montenegro. All isolates were identified as Ln. mesenteroides. Classical multilocus sequence tying (MLST) and core genome (cg) MLST revealed a high diversity of the Montenegrin Ln. mesenteroides cheese isolates. All isolates carried genes of the bacteriocin biosynthetic gene clusters, eight out of 16 strains carried the citCDEFG operon, 14 carried butA, and all 16 isolates carried alsS and ilv, genes involved in forming important aromas and flavor compounds. Safety evaluation indicated that isolates carried no pathogenic factors and no virulence factors. In conclusion, Ln. mesenteroides isolates from Montenegrin traditional cheeses displayed a high genetic diversity and were unrelated to strains deposited in GenBank.


Introduction
The genus Leuconostoc (Ln.) currently comprises 14 species and eight subspecies, all being Gram-positive, non-spore-forming, non-motile Lactic Acid Bacteria (LAB). Within the genus Leuconostoc, Ln. mesenteroides, Ln. pseudomesenteroides and Ln. lactis have their role in food fermentation and can be isolated from various food-related ecological niches, including beverages, meat and dairy products, and some plant materials, implying wide distribution and specialized adaptation to these diverse environments [1][2][3][4]. The first description of Ln. mesenteroides was by Van Tieghem in 1878 and was later proposed as the type strain [5]. In 1983, Ln. dextranicum and Ln. cremoris were reclassified as subspecies of Ln. mesenteroides due to the common properties they shared [6] and their high degree of relatedness shown by DNA-DNA hybridization [5].
In the food production, Leuconostoc spp. are usually applied as an adjunct culture in combination with the fast acid producing Lactococcus spp., as undefined mixed (DL)

Origin and Cultivation of Isolates
Bacteria were isolated from 13 different white brine ripened traditional cheeses collected during spring and summer 2019 from nine producers from the three municipalities Pljevlja, Šavnik, and Žabljak in the northern region of Montenegro (Supplementary Figure S1 Regions of Montengro), as described previously [30]. Briefly, for each sample, 20 g of cheese was added to 180 mL sterile 2% (w/v) sodiumcitrate solution and homogenized for 2 min in an Omni mixer (Omni International, Waterbury, CT, USA). For enumeration, serial dilutions in Ringer's solution were plated on de Man, Rogosa, Sharpe (MRS), and M17 agar plates (HiMedia, Mumbai, India) and incubated at 30 • C for 72 h. Isolates were biochemically characterized for their acidification and post-acidification ability, growth at different temperatures and NaCl concentrations, and their ability for lactose degradation. Isolates were stored in MRS and M17 medium supplemented with 15% glycerol (v/v) at −80 • C, and, for revitalization, subcultured three times in the respective broths (MRS and M17) at 30 • C overnight.
Up to ten single colonies were selected from MRS and/or M17 agar plates and subcultured for further processing. Species confirmation was carried out using matrix assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF-MS) [31] on a MALDI-TOF Microflex LT/SH with database MBT Compass IVD 4.2 according to the manufacturer's instructions (Bruker, Billerica, MA, USA).

Biochemical Testing 2.2.1. Growth at Different Temperatures
Overnight cultures grown at 30 • C were inoculated in MRS broth and incubated for 24 h at 4, 10, 15, 30 (optimal growth temperature, used as control), and 45 • C. The optical density was measured on a Jasco UV/VIS spectrophotometer V-730 (Jasco, Cremella, Italy) at 560 nm (OD 560 ).

Production of CO 2
The production of CO 2 from glucose was determined as described previously [32]. Briefly, bacterial cultures (50 µL) were added into a test tube containing 10 mL of suitable broth supplemented with 1% of glucose (Tokyo Chemical Industry, Tokyo, Japan) and a Durham's tube. After 24-48 h of incubation at 30 • C, gas production was assessed: if the gas accumulated in the Durham's tube to more than one third of its capacity, the result was considered positive.

Salt Tolerance
The evaluation of salt tolerance was performed as described previously [32,33]. Briefly, salt content of tubes containing 3 mL of MRS broth was adjusted to a final concentration of 2%, 3%, 4.5%, or 6.5% (w/v) NaCl. Tubes were inoculated with overnight cultures of the strains and incubated at 30 • C for 24 h. Bacteria grown in MRS broth without NaCl were used as controls. The optical density was measured on a Jasco UV/VIS spectrophotometer V-730 at 560 nm (OD 560 ).

Acidification and Post-Acidification Ability in Milk
Acidification and post acidification ability was tested as described previously by inoculating 50 mL ultra-high temperature (UHT) processed skimmed cow's milk (Imlek, Belgrade, Serbia) [34] with 50 µL of broth culture and incubated at 30 • C for 48 h. The resulting pH values, measured after 2, 4, 6, 8, 12, and 24 h post-acidification, were investigated by measuring the pH of the inoculated skimmed milk after 48 h of incubation [34]. pH of milk was measured with a WTW pH Meter inoLab ph 7110 (Xylem Analytics, Weilheim, Germany). Bacterial cells from a solid nutrient medium were transferred to empty Petri dishes. In addition, 1-2 drops of 3% hydrogen peroxide were dropped on the cells. Formation of air bubbles (production of oxygen due to the presence of catalase) was examined visually.
2.2.6. Ability for Formation of Exopolysaccharides (EPS) EPS formation ability was performed by identifying mucous-producing colonies grown on MRS agar plates. Additionally, mucoid colonies were stretched with an inoculation loop to detect the production of long filaments. If the mucus producing colonies or the formation of long filaments was missing, the strain was considered EPS negative.

Antimicrobial Susceptibility Testing
The Kirby-Bauer disk diffusion test was used to determine antimicrobial susceptibility of strain INF36 against erythromycin (15 µg, Bioanalyse, Ankara, Turkey) and of strain INF117 against ampicillin (10 µg) and penicillin (1U) (Bioanalyse, Ankara, Turkey). The disk diffusion patterns were evaluated according to the microbiological breakpoints for selected lactic acid bacteria as defined by EFSA [35,36].

DNA Extraction and Whole Genome Sequencing
For whole genome sequencing, high quality genomic DNA was isolated from overnight cultures grown in M17 or MRS using the MagAttract HMW DNA Kit (Qiagen, Hilden, Germany) according to the protocol for Gram-positive bacteria following the manufacturer's instructions (Qiagen). The amount of input DNA was quantified on a Lunatic instrument (Unchained Labs, Pleasanton, CA, USA). Ready to sequence libraries were prepared using the Nextera XT DNA library preparation kit (Illumina, San Diego, CA, USA) according to the instructions of the manufacturer. Paired-end sequencing with a read length of 2 × 300 bp using Reagent Kit v3 chemistry (Illumina, San Diego, CA, USA) was performed on a Miseq instrument (Illumina, San Diego, CA, USA) according to the instructions of the manufacturer (Illumina).
Subtyping of all 16 Ln. mesenteroides isolates was conducted in SeqSphere+ v7.2.3 using a recently published eight-gene multilocus sequence typing (MLST) scheme comprising the eight loci carB, groeL, murC, pheS, pyrG, recA, rpoB, and uvrC [25] and a newly defined ad hoc core genome multilocus sequence typing (cgMLST) scheme. For definition of the cgMLST scheme, a genome-wide gene-by-gene comparison was performed using the complete genome of strain Ln. mesenteroides subsp. mesenteroides ATCC 8293 as a reference genome and 18 complete genomes of Ln. mesenteroides available in GenBank as query genomes (Supplementary Table S1 cgMLST_scheme_data). The resulting ad hoc cgMLST scheme comprised 960 core genome target genes and 935 accessory genome target genes. Sixty genes of the ATCC 8293 genome were discarded since they did not fulfil  Table S1 cgMLST_scheme_data). Seven of the eight MLST targets [25] were also part of our cgMLST scheme while one target (uvrC) belonged to the accessory genome. For phylogenetic analysis, a minimum spanning tree (MST) was calculated based on the defined 960-target cgMLST scheme. Based on the 17 allelic differences observed between the query strains NBRC107766 and ATCC19254 (both are the same strain but originating from different repositories and sequenced at different laboratories), an arbitrary complex type threshold of 45 allelic differences was applied for detection of related isolates.

Gene Presence and Absence
Detailed analysis of additional genomic information was performed for all 16 genomes and 12 reference strains: KMB608, ATCC8293, 406, 213M0, WC0331, DRC1506, ATCC19254, LN08, DSM20484, LN32, NBRC100495, and LbE15. The contigs of each assembly were filtered for a minimum length of 1000 nucleotides. Genes were predicted using prodigal v2.6.3 [49] with default parameters, and orthologous groups were calculated with OrthoFinder v1.4.2 [50]. Orthologous groups with differences in presence/absence were selected for visualization. Additionally, a species tree was created from the orthologous groups using default parameters for Orthofinder. Annotation of predicted genes was performed using NCBI BLAST+ v2.10.0 [51].

Biochemical Properties
Sixteen Ln. mesenteroides isolates were successfully isolated and identified by MALDI-TOF-MS from 13 Montenegrin brine cheeses (Supplementary Figure S1 Regions of Montenegro, Table 1). All isolates were Gram-positive, catalase negative, and had no ability for formation of exopolysaccharides (EPS). The mean number of viable bacterial cells was in the range between 41 to 92 cfu/g cheese. The acidification ability of the tested strains was in the range of pH 4.3-5.9 (24 h), while post acidification ability was between pH 4.2-5.7 (48 h). After 48 h of incubation at 30 • C, the pH in broths supplemented with 1% lactose was in the range of 4.4 to 5.5. Six tested isolates exhibited an ability to grow at 6.5% NaCl concentration, 11 strains at 4.5%, while three strains did not grow at 2% NaCl concentration. Both isolates obtained from producer A (INF2 and INF3b) showed growth ability at a salt concentration of 2%, but did not grow at higher salt concentrations. Isolate INF82b (producer D) was only growing at a salt concentration of 4.5%. Ten isolates showed an ability to grow at 45 • C, three isolates INF3b (producer A), INF36 (producer B), and INF46 (producer C) had the ability to grow at 4 • C (Table 1).

Detection of Bacteriocin and Secondary Metabolite Genes
All isolates but one carried genes of the bacteriocin biosynthetic gene clusters, ten out of 16 carried a betalactone biosynthesis gene, and all carried a type III polyketide synthase gene (T3PKS) ( Table 2). Eight out of our 16 Leuconostoc strains carried genes necessary for citrate metabolism (citCDEFG operon) (Figure 3), 14 isolates had the acetoin reductase gene butA and all isolates carried genes necessary for the biosynthesis of branched aminoacids acetolactate synthase alsS and isoleucine synthesis operon ilv.

Safety Evaluation
Safety evaluation indicated that none of the 16 isolates carried pathogenic and virulence factors. Five isolates carried plasmids. The CARD database identified a C656T transition in the 23S rRNA gene of isolate INF36, conferring resistance to erythromycin and clindamycin (100% identity, 10% coverage) using the CARD database. ResFinder did not detect this mutation in the 23S rRNA gene of isolate INF36. Disk diffusion testing revealed susceptibility of strain INF36 against erythromycin (Table 3). Isolate INF117 carried a partial blaTEM with 99.04% identity and 60.51% coverage to blaTEM-141 (and others) using ResFinder. CARD analysis (using the loose hits criteria) confirmed the presence of a partial blaTEM in isolate INF117 (97.5% identity and 55.94% coverage to several blaTEM) , and INF98 obtained from three cheeses of the same type from producer D differed by five to a maximum of 860 alleles in their core genome (Table 1, Figure 2). Isolates INF82a and INF82b obtained from the same cheese from producer D differed by 860 alleles (Figure 2). The cgMLST based comparison of the Montenegrin Ln. mesenteroides isolates with isolates available from GenBank revealed an allelic difference from 138 to a maximum of 646 alleles (Figure 2). The closest related strains available from GenBank were Ln. mesenteroides strain 406 and 213MO, two strains isolated from Mongolian mare milk differing by 138 alleles, respectively, and 165 alleles from strain INF103 from producer E (Figure 2). CgMLST revealed 295 to 865 allelic differences between the Montenegrin isolates and the Ln. mesenteroides subspecies reference strains (Supplementary Table S4 genome data). CgMLST results were concordant with TYGS and ANI analysis (Supplementary Table S2 Table S4 genome data, Figure 1).

Detection of Bacteriocin and Secondary Metabolite Genes
All isolates but one carried genes of the bacteriocin biosynthetic gene clusters, ten out of 16 carried a betalactone biosynthesis gene, and all carried a type III polyketide synthase gene (T3PKS) ( Table 2). Eight out of our 16 Leuconostoc strains carried genes necessary for citrate metabolism (citCDEFG operon) (Figure 3), 14 isolates had the acetoin reductase gene butA and all isolates carried genes necessary for the biosynthesis of branched aminoacids acetolactate synthase alsS and isoleucine synthesis operon ilv.

Orthofinder Analysis
Orthofinder analysis assigned 57,889,803 genes (97.9% of total) to 3080 orthologous groups. The number of orthologous groups identified for our isolates ranged from 1873 to 2111 (mean: 1979). There were 1280 orthogroups present in all Leuconostoc subspecies and 1091 of these orthologous groups consisted entirely of single-copy genes ( Figure 3).

Discussion
Traditional food products are an acknowledged part of the identity, culture, and heritage of a country and, in addition, they may serve as a valuable resource for future food production [16,17,52]. To explore the diversity of Leuconostoc spp. Isolates, we characterized in this study 16 Leuconostoc mesenteroides (Ln. mesenteroides) isolates obtained from 13 traditional Montenegrin brine cheeses from nine different producers by whole genome sequencing (WGS). To the best of our knowledge, our study comprises the largest set of Ln. mesenteroides strains from traditional cheese characterized biochemically and by WGS.
The biochemical abilities of strains were in the range of expected technological significance. Growth at low/high salt and or low/high temperatures of single isolates could not be linked to genotypes.
Nevertheless, genome sequencing is considered as a suitable tool for the selection of particularly suitable strains used as adjunct cultures for the production of traditional cheese [15,21,22,26]. Leuconostoc species were "Generally Regarded As Safe" (GRAS) organisms by the US Food and Drug Administration (FDA) and the European Food Safety Authority (EFSA) for food production [53,54]. Reports about clinical infections caused by Ln. mesenteroides [55] in combination with the role of LAB as a reservoir for antimicrobial resistances have changed the uncritical classification of LAB as GRAS. Antimicrobial resistance and multi-resistance have been described for several Ln. mesenteroides isolates [36,56]. Safety evaluation of the Montenegrin brine cheese isolates revealed that they carried no pathogenic factors and no virulence genes. However, two isolates carried antimicrobial resistance determinants. One isolate carried a fragment of the blaTEM, which also exists in the genomes of several other Ln. mesenteroides strains including reference strains ATCC 8293. For the second strain, we found for the first time in a Leuconostoc spp. isolate a C656T point mutation in the 23S rRNA gene. This mutation was detected for the first time in Clostridioides difficile (Clostridium difficile) conferring high resistance to erythromycin and low resistance to clindamycin [57]. However, this point mutation was not found via ResFinder analysis. Phenotypic testing revealed that both strains were susceptible against the respective antibiotics, which can be explained by the fact that the detected antimicrobial resistance gene fragments are non-functional.
Analysis of genes essential for the biosynthesis of secondary metabolites revealed that all but one strain carried genes of the bacteriocin biosynthetic gene clusters, which is an important natural feature in cheese production through inhibiting the growth of several important pathogenic bacteria. All strains carried a type III polyketide synthase gene essential for the biosynthesis of important secondary metabolites [43] and all isolates carried genes essential for the biosynthesis of branched aminoacids that contribute to flavor. Interestingly, only 50% of our Montenegrin cheese isolates carried the cítrate operon, which is essential for citrate metabolism and a characteristic of dairy Ln. mesenteroides strains [20].
To investigate strain diversity, a previously published MLST scheme [25] and our novel ad hoc core cgMLST scheme were used. Our ad hoc cgMLST scheme comprised 960 core targets in comparison to a recently published cgMLST scheme, which comprised 999 core genes [19]. We used a so-called "hard core genome" definition and 19 instead of 17 complete genomes for core genome generation, which resulted in the observed lower number of 960 core genome targets. All tested isolates except four had more than 97% good core genome targets. Type Strain Server Analysis (TYGS) [38] identified most of our isolated strains as subspecies Ln. mesenteroides mesenteroides, confirming recent findings that subspecies Ln. mesenteroides mesenteroides is better adapted to cheese production than other subspecies and also than subspecies Ln. mesenteroides cremoris, which is mainly used in commercial starter cultures [3]. In addition to subspecies Ln. mesenteroides mesenteroides, TYGS analysis assigned the remaining six Montenegrin brine cheese isolates to subspecies Ln. mesenteroides cremoris, Ln. mesenteroides jonggajibkimchii, and Ln. mesenteroides dextranicum. Isolates identified as Ln. mesenteroides mesenteroides and Ln. mesenteroides dextranicum by TYGS clustered appropriately with the respective subspecies available from GenBank when using dDDH, ANI, or our cgMLST approach. Four isolates were not accurately assignable by TYGS, ANI, or cgMLST to a Leuconostoc mesenteroides subspecies. All these methods placed them between different Leuconostoc subspecies. Congruent with TYGS and ANI, cgMLST analysis placed two of the Montenegrin isolates between strains previously identified as Ln. mesenteroides mesenteroides [20] and Ln. mesenteroides cremoris [26]. Two other Montenegrin isolates were placed between Ln. mesenteroides dextranicum strains isolated from Italian cheese [26] and Ln. mesenteroides mesenteroides strains isolated from sheep milk in Slovakia [unpublished, https://www.ncbi.nlm.nih.gov/biosample/SAMN09398940] (accessed on 25 February 2020). From this Italian study [26], diverse Ln. mesenteroides subspecies were recently reassigned to Ln. mesenteroides lineage M4 and Ln. mesenteroides lineage M1 [20], which confirmed our core genome results, i.e., clustering of these strains on the same branches. In contrast, all other Ln. mesenteroides cremoris strains available from GenBank clustered together and were unrelated to the Italian and our isolates. This awkward subspecies assignment indicates that Ln. mesenteroides lineage M1 and M4 strains have some specific features from diverse subspecies [20]. The four Montenegrin isolates that were not assignable to a subspecies (and related to lineage M1 and M4 isolates according to Frantzen et al., 2017) [20] had, in comparison to the other Montenegrin cheese isolates, a lower percentage of good core genome targets (≤97%). This lower number of core genes is also characteristic for Ln. mesenteroides cremoris ATCC19254 and Ln. mesenteroides jonggajibkimchii DRC1506 due to the loss of genes of these two Leuconostoc subspecies as a consequence of adaptation to their ecological niches [19,20]. Our results support recent findings that a taxonomic revision of the genus Leuconostoc should be considered and that some of the isolates available in GenBank have been assigned to the wrong subspecies due to the use of phenotypic properties [20,58].
MLST and cgMLST analysis revealed that the 16 Montenegrin Ln. mesenteroides cheese isolates showed a high diversity, which contrasts with isolates obtained from commercial starter cultures, which show a low genetic diversity [20,25]. This high diversity underlines the importance of traditional food products as a valuable source for strains with unique and interesting features to be used in the dairy industry as novel starter cultures or for production of functional dairy foods [15,18,52]. The observed high diversity of our Montenegrin Ln. mesenteroides isolates and, in addition, their non-relatedness to genomes of strains available from GenBank reflects that no commercial starter cultures were used to produce the traditional Montenegrin brine cheeses included in this study. Although we investigated only a small number of isolates, the observed diversity corresponds with results from previous studies showing a high strain diversity of Ln. mesenteroides obtained from a variety of food sources [23][24][25]. Isolation of different strains from the same cheese product of the same producer and the diversity of strains obtained from the same cheeses also indicate a variable composition of the Leuconostoc population in traditional cheeses, which may also affect the sensory characteristics and quality of the final products. However, more extensive studies are necessary to prove these findings. The differences between isolates from the same producer as well as from the neighboring regions might be explained by the specific natural plant biodiversity, which is strongly supported by the fact that, in Montenegro, the feeding of animals is still to the highest extent based on natural grazing, with a relatively low percentage of concentrate feeding [52]. Thus, genome-based characterization of Leuconostoc spp. may present an innovative tool to prove the origin of strains and to ensure consistent manufacture of high quality and safe traditional food products [19,20]. In the future, with the ongoing progress in multi-omics technologies, the prediction of gene functions, metabolic pathways, and inter-microbe interactionsparticularly within starter cultures-may allow a selection of strains based on specific marker genes more efficiently [28]. Currently, the accurate prediction of metabolic patterns and the comparison of metabolic capacities of organisms using genome-scale metabolic models is still a major challenge or even impossible due to the lack of complete genomes [8].
In conclusion, Ln. mesenteroides isolates from the Montenegrin traditional brine cheeses show a high genetic diversity, which can be explained by the high level of plant biodiversity in Montenegro and the fact that the feeding of animals is mainly based on natural grazing. WGS based strain-typing proves the Montenegrin origin and endogenous status of strains. With the ongoing progress in predictive microbiology, indigenous cheese isolates with specific properties will become selectable and become available for standardization of the production, to preserve designation of origin as well as to create added-value sensory characteristics.