The present project aims to investigate the evolutionary relationships within the genus Leuconostoc through a phylogenomic approach and to verify the consistency of phylogenetic relationships with the current taxonomy. This genus holds a group of Lactic Acid Bacteria (LAB) belonging, together with the genera Convivina, Fructobacillus, Oenococcus, and Weissella, to the family of Leuconostocaceae (www.bacterio.net; Nieminen et al., 2014; Bello et al., 2022), in its turn included in the order Lactobacillales, based on the average amino acid identity of core proteins (cAAI) (Zheng et al., 2020). The genus Leuconostoc encompasses 17 species according to the List of Prokaryotic names with Standing in Nomenclature (https://lpsn.dsmz.de/genus/leuconostoc) (Parte et al., 2020). They are saccharolytic bacteria that catabolize carbohydrates to lactic acid through heterolactic fermentation and inhabit a variety of niches where carbohydrate-based substrates are available, such as plants, plant-derived matrices, silage, fermented foods (e.g., dairy products, fermented dough, milk, vegetables, and meats), spoiled foods, and sewage (Dellaglio et al., 1995). In some cases, species are associated with a specific habitat, such as vegetables (Yu et al., 2020), meats (Candeliere et al., 2021), or other foods (Vedamuthu, 1994). The genus is considered safe products and has been accorded the status of “generally recognized as safe” (GRAS) (Ogier et al., 2008). Fermentation of suitable substrates has attracted attention in recent years to produce antimicrobial extracts (Venegas-Ortega et al., 2019; Ahmadi-Ashtiani et al., 2020; Costa et al., 2021). Strains of Leuconostoc spp. have found promising industrial applications in the preparation of biopreservative systems by fermentation of different substrates. Leuconostoc are used as starters in food and beverage fermentation to improve the nutritional and sensorial properties and to extend the shelf life (Shin and Han, 2015). Moreover, bioactive extracts obtained from vegetables fermented with Leuconostoc spp. found application in the formulation of innovative cosmetics (INCI Name: Leuconostoc/Radish Root Ferment Filtrate) (Ahmadi-Ashtiani et al., 2020).

The early taxonomy of LAB was based on phenotypic and morphological features, then chemotaxonomic criteria such as DNA–DNA hybridization and G+C content became the reference for species assignment (Vandamme et al., 1996). With the advent of phylogenetic taxonomy (Woese and Fox, 1977), the sequence of the gene encoding 16S rRNA turned into the gold standard for taxonomic and phylogenetic analysis (Lane et al., 1985; Stackebrandt and Goebel, 1994; Stackebrandt and Ebers, 2006), and allowed the systematic study of evolutionary relationships among prokaryotes. A phylogenetic tree of a conserved gene that is assumed to be vertically inherited is expected to produce the statistical trend or the real genealogy of evolving entities. However, due to the limited size of individual genes, multiple substitutions, parallel, convergent, or reversal events, and horizontal transfer of DNA, the strength of the phylogenetic signal in single molecules is often too low to infer the evolutionary relationships underlying species differentiation. Therefore, multiple phylogenetic markers have been applied to obtain a well-resolved and informative tree, either by the concatenation of genes, aiming to average their phylogenetic information, or by the corroboration of individual phylogenetic signals (Whelan and Morrison, 2017).

With the advent of high-throughput sequencing, an increasing amount of whole bacterial genomes had been accumulated, thus plenty of information became available to improve the resolution of bacterial diversity and the accuracy of phylogenetic reconstruction. In this context, core genome phylogenesis is widely used to infer phyletic lines in the evolutionary history of prokaryotic species (Stott and Bobay, 2020). The average nucleotide identity (ANI) of the genes shared between two genomes was introduced as the gold standard for the delineation of bacterial species (Richter and Rosselló-Móra, 2009; Chun and Rainey, 2014) and has been recently proposed as a tool to determine statistically supported phylogenies (Gosselin et al., 2022).

In the present study, the taxonomy of 221 Leuconostoc genome sequences belonging to 17 species has been preliminarily investigated based on ANI to measure genome similarity. ANI, which provides an average measure of similarity across homologous regions shared by a pair of genomes, is a major metric used for this purpose (Palmer et al., 2020), considering the ANI threshold of 95% for species delineation (Richter and Rosselló-Móra, 2009). The subspecies delimitations were also investigated, considering a threshold of 98% (Minias et al., 2020; Pearce et al., 2021), although it is not recognized as a standard in the taxonomy of prokaryotes. For the strains presenting low ANI values, the Average Aminoacidic Identity (AAI) was utilized to confirm the inclusion in the genus Leuconostoc and the different putative species, considering the threshold of 55–60% and 85–90%, respectively (Rodriguez-R and Konstantinidis, 2014). The evolutive genealogy of the genomes was reconstructed utilizing a minimum evolution tree of ANI, computed according to Gosselin et al. (2022), and the core genome alignment.

Phylogenetic relationships were also inferred using the housekeeping genes encoding phenylalanyl-tRNA synthase alpha-subunit (pheS) and RNA polymerase alpha-subunit (rpoA) (Das et al., 2014), in addition to the 16S rRNA gene. rpoA was previously exploited for the identification of L. falkenbergense strains at the species level (Wu and Gu, 2021a), while pheS is one of the targets in the Multilocus Sequence Analysis (MLSA), successfully applied for differentiation of species of the genus Leuconostoc (Rahkila et al., 2014). As a whole, multiple phylogenomic approaches have been applied to investigate evolutionary relationships within the genus Leuconostoc, following a preliminary reclassification of the strains in putative species according to ANI comparison. Split decomposition analysis has been carried out to check the ambiguities in the inferred phylogenesis, evaluate consistency among the output of the diverse analysis, and determine the approach that provides the strongest phylogenetic signal (Whelan and Morrison, 2017).

In total, 221 Leuconostoc genome sequences available on 1 October 2021 were retrieved from the NCBI database. The accession numbers are reported in Supplementary Table S1. Only genomes obtained from pure culture sequencing were used, while metagenome-assembled genomes were discarded. The genomes were inspected for completeness and contamination with CheckM (Parks et al., 2015). Genomes were also checked for the presence of the molecular signatures characterizing the genus Leuconostoc (Bello et al., 2022). The genomes were annotated with Prokka (Seemann, 2014) to calculate the pangenome using Roary (Page et al., 2015), with the minimum percentage identity parameter set at 80%.

The genes encoding rpoA, pheS, and 16S rRNA were extracted from each genome and aligned with Clustal Omega on EMBL-EBI website (Sievers et al., 2011; Madeira et al., 2019). The Total Average Nucleotide Identity (ANI) was calculated following the method described by Gosselin et al. (2022). This method allowed us to calculate the bootstrap values from 100 ANI distance matrices replicates, thanks to a process of genome segmentation and random segment selection utilized to calculate the desired number of replicates. The AAI between all the genomes was calculated with AAI Calculator (http://enve-omics.ce.gatech.edu/aai/) by all-against-all approach (Rodriguez-R and Konstantinidis, 2016).

The alignments of 16S rRNA, rpoA, and pheS genes and the 247 core genes identified by Roary out of a total of 25,377 genes (Supplementary Spreadsheet 1) were utilized to generate phylogenetic trees constructed with the maximum likelihood method using the RAxML tool with 100 bootstrap replicates (Stamatakis, 2014). The ANI distance matrix was used to build a tree using the script provided by Gosselin et al. (2022). In this script, the balanced minimum evolution algorithm implemented in the FastME function of the R package APE (Paradis et al., 2004) was applied to generate phylogenies for each distance matrix (Desper and Gascuel, 2002), whereas the function plotBS of the R package Phangorn (Schliep, 2011) was exploited to map support values onto the tree. The AAI distance matrix was used to compute an Unweighted Pair Group Method with the Arithmetic mean (UPGMA) unrooted phylogenetic tree with the DendroUPGMA tool (Garcia-Vallvé et al., 1999). The phylogeny was further inferred with SplitsTree v. 4.18.2 (Huson and Bryant, 2006) with a neighbor net drawing and Jukes–Cantor correction for alignment-derived trees (Bandelt and Dress, 1992; Huson and Bryant, 2006).

The tree of ANI was compared with those constructed from the alignment of 16S rRNA, pheS, rpoA genes, or of the core genome, utilizing the following “generalized” Robinson–Foulds metrics: the Jaccard–Robinson–Foulds (JRF), computed with k = 1 (Nye et al., 2006), and the information-based measures of Mutual Clustering Information (MCI) and Shared Phylogenetic Information (SPI) (Smith, 2020). JRF, MCI, and SPI were computed with the R package “TreeDist,” archived at https://dx.doi.org/10.5281/zenodo.3528123. Comparison relied on matching each informative split within a tree (i.e., an internal branch of the tree, with at least two leaves at each extremity) with an informative split within the other tree. For each metric, the comparison yielded a normalized similarity score in the range 0–1 between two trees and a score for each match of paired tree splits.

Bacterial species can be defined within populations of similar strains on the basis of maximum likelihood methods that determine the point of transition of the evolutionary processes. As a result of common evolution, the strains belonging to the same species share genetic features and thus biological properties. Species delimitations constitute an important challenge in biodiversity studies, mostly for genera such as Leuconostoc that have an important role in several industrial and food fermentation processes (Buckenhüskes, 1993; Caplice and Fitzgerald, 1999; Steinkraus, 2002), to produce fermented sausages, fermented vegetables and cereals, and dairy products (e.g., butter, cream, fresh and raw milk, cheeses), and to obtain novel ingredients for cosmeceutical formulations (Ahmadi-Ashtiani et al., 2020).

Taxonomic classifications of the 216 Leuconostoc strains based on the core genome sequences indicated that many strains attributed to the same Leuconostoc species laying onto paraphyletic branches necessitate reclassification and that the taxonomy of the genus is not entirely resolved. During the revision of the manuscript, a study suggesting the reorganization of Leuconostoc taxonomy was published (Kumar et al., 2022). In that manuscript, presenting some overlapping with our survey also in terms of procedures, the proposal of a taxonomic update of the genus Leuconostoc was based on ANI and core genome analyses and is consistent with the outcome of our investigation. The phylogenomic reconstruction of the genus herein presented, integrating the information of ANI, core genome, and housekeeping genes complemented the taxonomic delineation with solid information on the phylogenetic lineages evolved within the genus Leuconostoc.

In this study, the phylogenomic trees built aligning housekeeping single genes corroborated ANI results. The topology of the phylogenetic trees consistently indicates the remote diversification of most of the species, with a few exceptions. Split decomposition analysis and comparison of the ANI tree with the other ones revealed robust consistency for most of the branches identified, with strong signals for the majority of the nodes identified by ANI, core genome, AAI, rpoA, and pheS genes, that generally presented very high bootstrap values. 16S rRNA gene yielded weak signals and thus was not suitable to define taxonomy and phylogeny within this genus. According to the consistent results of ANI and phylogenomic trees, the strains currently attributed to L. mesenteroides lay onto paraphyletic branches and need to be split into two separate species, corresponding to G16 and G18. Interestingly, these groups are already distinct in GTDB where they are referred to as Leuconostoc mesenteroides_B and Leuconostoc mesenteroides, respectively, and have been delineated also by Kumar et al. (2022). G18 encompasses most of the strains currently assigned to L. mesenteroides, including the type strain, and all the strains belonging to the subspecies L. mesenteroides subsp. cremoris, L. mesenteroides subsp. dextranicum, L. mesenteroides subsp. jonggajibkimchii, and L. mesenteroides subsp. mesenteroides. Nonetheless, the strains belonging to different L. mesenteroides subspecies present a genome similarity always >99.1%, higher than the threshold of 98% utilized to delineate subspecies in other bacterial taxa (Minias et al., 2020; Pearce et al., 2021). Accordingly, both in this study and in previous ones the conserved genes failed to differentiate between the subspecies of L. mesenteroides (Ricciardi et al., 2020).

G18 is closely related to G17, which includes strains of the species L. suionicum, previously recognized as subspecies L. mesenteroides subsp. suionicum (Jeon et al., 2017). The two strains belonging to G16, currently assigned to L. mesenteroides, should form a new species that would require characterization and formal definition (Das et al., 2014; Ramasamy et al., 2014). The strains currently assigned to L. pseudomesenteroides should be split into two species, G13 and G14. The former harbors the type strain L. pseudomesenteroides NCDO 768^T, and the latter encompasses also the strains of L. falkenbergense, including the type strain L. falkenbergense LMG 10779^T. L. falkenbergense has been already recognized as a different species phylogenetically related to L. pseudomesenteroides (Wu and Gu, 2021a) and recently added in GTDB.

The strains currently attributed to L. lactis should be split into two separate species, corresponding to the groups G9 and G10. Consistently, in GTDB the groups are already distinct and referred to as Leuconostoc lactis_A and Leuconostoc lactis. G10 encompasses the type strain L. lactis JCM 6123^T, the strain KFRI01 named L. garlicum, a species that is not formally recognized, and L. citreum 1300_LCIT. G9 harbors L. lactis KCTC 3773, previously identified as the type strain of L. argentinum, species that is not anymore accepted (Dicks et al., 1993; Vancanneyt et al., 2006). Within both G10 and G9, subgroups of strains putatively ascribable to different subspecies can be identified, that present ANI scores between each other lower than 98%. Therefore, further phenotypic and genotypic studies could investigate whether subspecies have to be created within G10 and G9 groups. The strains of L. gelidum and L. gasicomitatum fall into G6 and G7, respectively, which include the type strains L. gelidum subsp. gelidum KCTC 3527^T and L. gasicomitatum LMG_18811^T. This classification confirms the observations by Wu and Gu (2021b) that rejected the proposal to recognize L. gasicomitatum as a subspecies of L. gelidum (L. gelidum subsp. gasicomitatum) (Rahkila et al., 2014). Moreover, G7 includes also the type strain L. inhae KCTC 3774^T, suggesting that L. inhae and L. gasicomitatum could belong to the same species. No clear correlation could be established between the phylogenetic relationships and the source of isolation of the strains, except for G2 (L. carnosum) and G14 (L. falkenbergense) that encompassed mostly strains isolated from meat and dairy products, respectively (Supplementary Table S1) (Raimondi et al., 2018; Wu and Gu, 2021a). The evolutive force that shaped the speciation of Leuconostoc deserves deeper comparative genomics and functional analysis.

The results of this study confirmed and deepened the evidence on the evolutionary relationship among Leuconostoc species and may provide a basis for a possible future reorganization of the genus, as summarized in Supplementary Table S1. However, any update in the organization of Leuconostoc genus, such as the creation of new species and/or subspecies, would require, all the biochemical and physiological data based on strains survey that support the creation of new species and the reclassification of others.

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author/s.

Conceptualization: SR, AA, and MR. Investigation: SR, FC, and GS. Data curation: FC, AA, SV, and SC. Formal Analysis: AA. Visualization: SR, FC, and AA. Writing—original draft: All authors. Funding acquisition: MR and SV. Supervision: MR. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.