A Genome-Based Species Taxonomy of the Lactobacillus Genus Complex

We downloaded all 2,558 genomes that belong to the Lactobacillaceae and Leuconostocaceae families from GenBank (see Table S1 in the supplemental material) and identified single-copy core genes (SCGs) in this data set. Our strategy for SCG extraction consisted of four steps: (i) full gene family clustering on a random subset of 30 genomes, (ii) selection of candidate SCGs from these gene families using a very soft single-copy presence threshold, (iii) search for the candidate SCGs in all genomes, and (iv) filtering of the definitive SCGs from the candidates using a fixed single-copy presence threshold. After visual inspection of the single-copy presence of the candidate SCGs in the full genome data set (see Fig. S1 in the supplemental material), we retained candidate SCGs with >95% single-copy presence. This strategy resulted in 411 SCGs sensu lato (“sensu lato” because we did not enforce 100% single-copy presence). We performed quality control of all genomes based on this set of SCGs (Fig. 1). More specifically, we determined two quality measures for each genome: completeness, defined here as the percentage of sensu lato SCGs that are present in the genome, and redundancy, defined as the percentage of sensu lato SCGs showing two or more copies in the genome. The large majority of genomes had a completeness close to one and redundancy close to zero. Enforcing a minimum completeness of 90% and maximum redundancy of 10% resulted in 2,459 high-quality genomes, or 96.1% of all genomes.

Some of the genomes with low estimated completeness belonged to the genus Sharpea, which is part of the family Lactobacillaceae according to the NCBI taxonomy but seems to be phylogenetically quite distant to it (8). Other low-completeness genomes were flagged as “incomplete” in the NCBI database or were derived from metagenome data. Interestingly, some genomes belonging to Lactobacillus iners, the LGC species with the smallest genome, had very high completeness values. This suggests that genome reduction did not (significantly) influence our completeness estimates. As an extra check, we compared our completeness and redundancy estimates to the estimates available in the Genome Taxonomy Database (GTDB), where possible. There was a good correspondence between both estimates, with the exception of a small number of genomes (among which were the aforementioned Sharpea genomes) that had low completeness estimates by our pipeline, but high estimates in the GTDB (Table S1). This was not unexpected, since our pipeline uses LGC-specific core genes for quality estimation, while the GTDB uses a set of taxon-specific core genes.

We compared every genome to every other genome using the core nucleotide identity (CNI), which we calculated from the supermatrix of SCGs. While average nucleotide identity (ANI) is the most commonly used genome similarity measure, it suffers from conceptual and practical drawbacks. First, since it is based on all genomic regions that are orthologous between two genomes, these regions will be different for each pairwise comparison. Thus, the type and amount of information are different across genome comparisons. Second, ANI is slow to compute for large data sets because it scales quadratically with the number of input genomes. The recently introduced fastANI (9) is orders of magnitude faster, but it suffers from the same conceptual problem and remains an approximation. Thanks to our novel approach for the rapid extraction of SCGs, computing CNI similarities becomes feasible in almost linear time.

The density of these pairwise CNI values was strongly bimodal (Fig. 2A): there was one sharp peak with CNI values of >94%, and a flatter peak (including multiple subpeaks) with values of <90%. This indicated that two genomes were either very closely related (high CNI) or relatively distantly related (low CNI). Therefore, genomes were clustered into de novo species by grouping together genomes with similarities equal to or higher than 94%, using single-linkage clustering. This resulted in 239 genome clusters.

In single-linkage clustering, an object is added to a cluster when it exceeds a similarity threshold to at least one object within that cluster (as opposed to, for example, an average similarity threshold used in average linkage clustering). One problem that can arise with this strategy is that individual similarities within clusters can potentially be low, much lower than the clustering cutoff used. To assess whether this problem occurred for our clusters, for each cluster we calculated the minimum CNI between genomes within the cluster and the maximum CNI between genomes in the cluster and genomes outside of the cluster. Almost all within-cluster CNI values were higher than 94% (Fig. S2). There was only one exception: two genomes in cluster 231 had a CNI of 93.6%, but even this value was very close to the species boundary. Furthermore, all 239 genome clusters met the criterion of exclusivity: all their within-cluster CNI values were higher than their between-cluster CNI values. This included cluster 231 (later identified as Lactobacillus sakei): its highest between-cluster CNI value was 81.6% (to cluster 221, later identified as Lactobacillus curvatus).

In addition to the CNI, we also calculated a second similarity measure between all our genomes: the fastANI similarity, which is a fast approximation of ANI. Just like the CNI distribution, the fastANI distribution was strongly bimodal (Fig. 2B). Its high-value peak covered the range 93% to 100%. A scatterplot of the CNI and fastANI similarities (Fig. 2C) showed that both measures were strongly correlated for closely related genomes (around or above the species cutoff) and that the fastANI values were generally lower than the CNIs in this range. For more distantly related genomes, with CNI of <80%, the values were clearly less correlated, and the fastANI values were (much) higher than their CNI counterparts. The CNI species delimitation cutoff of 94% roughly corresponded to a fastANI cutoff of 93%. At these cutoffs, both similarity measures were in strong agreement on whether two genomes belonged to the same species or not.

To assess whether the species level held properties absent from the other taxonomic ranks, we calculated the exclusivity of clusters at various CNI cutoff values. The results showed that for CNI clustering cutoffs between 88% and 96%, all clusters were exclusive (Fig. S3). This was not the case for most cutoffs outside of this range.

We then attempted to label all 239 genome clusters (de novo species) with the name of a published species. We did this in two stages: (i) identification of genomes belonging to type strains of published (sub)species (further called “type genomes”) and (ii) for the remaining clusters, extraction of 16S rRNA genes from the genomes and comparison to 16S rRNA genes from type strains of published (sub)species (further called “type 16S genes”).

In the first cluster-naming stage, we identified type genomes present in our data set (Table S2). First, we cross-referenced the strain names that were listed in the NCBI assembly database with type strain names of validly published (sub)species found on the List of Prokaryotic Names with Standing in Nomenclature (LPSN) (15, 16) or Prokaryotic Nomenclature Up-to-date (PNU) (17). This resulted in the identification of 367 type genomes within our data set, belonging to 226 unique names (species or subspecies) from 211 species. Several type genomes for the same name were often found because many type strains have been sequenced multiple times. Next, we supplemented this type genome list by manually scanning the literature for extra type genomes, including also species that were not published validly (i.e., not in an official taxonomic journal). This yielded type genomes for an additional 11 species. There were different reasons why these manually identified type genomes had not been detected by our automated approach. Some of them were not represented on LPSN and PNU because they were not validly published. For other type genomes, their species were validly published, but too recently for them to be already on LPSN or PNU. Then there were some special cases. Lactobacillus nuruki (18) was a validly published name but was not listed on LPSN or PNU, probably because its genus name was misspelled as “Lactobacilus” in the NCBI database and in the paper. Convivina intestini was missed by our automated approach because we did not include this genus in our search, while it does appear to be part of the LGC, more specifically of the family Leuconostocaceae (16). The Lactobacillus musae type genome was missed because its type strain name (strain 313) consists only of numbers and was consequently filtered out by our pipeline. Finally, Leuconostoc lactis had a type genome in the NCBI database that did not pass our quality control, but we could also identify a type genome of Leuconostoc argentinum, a later synonym of Leuconostoc lactis (Leuc. lactis). This type genome had not been identified by our automated pipeline because it filtered out rejected (sub)species names.

Of our set of 239 genome clusters, 210 contained one or more type genomes (Table S1). The large majority of these clusters (200 of them) contained type genomes from exactly one species and could thus be given the name of that species. Ten clusters contained type genomes from two or three different species, indicating that the type strains of these species were more closely related than the 94% CNI cutoff (Table 1). These clusters were cluster 135 with Lactobacillus amylotrophicus and Lactobacillus amylophilus, cluster 141 with Lactobacillus kimchii and Lactobacillus bobalius, cluster 173 with Lactobacillus micheneri, Lactobacillus timberlakei, and Lactobacillus kosoi, cluster 178 with Weissella jogaejeotgali and Weissella thailandensis, cluster 179 with Lactobacillus fructivorans and Lactobacillus homohiochii, cluster 218 with Pediococcus acidilactici and Pediococcus lolii, cluster 225 with Leuconostoc gelidum, Leuconostoc inhae, and Leuconostoc gasicomitatum, cluster 233 with Lactobacillus gasseri and Lactobacillus paragasseri, cluster 236 with Leuconostoc suionicum and Leuconostoc mesenteroides, and cluster 238 with Lactobacillus casei and Lactobacillus zeae. We assigned the name of the species that was described first to these clusters. In one situation, two subspecies of the same species had a type genome in different genome clusters: Lactobacillus aviarius subsp. aviarius and Lactobacillus aviarius subsp. araffinosus. We named the genome cluster with the latter type strain Lactobacillus araffinosus.

In the second cluster-naming stage, we manually named the 29 remaining genome clusters in which we had not found type genomes (Table 2). We could assign a species name to 13 genome clusters with relative certainty because two sources of information were in agreement: (i) the NCBI species labels attached to the genomes in the cluster and (ii) matches to our in-house-constructed database of 16S rRNA genes of validly published species not yet used to name a cluster via type genomes. For these genome clusters, we included a question mark in the species name to indicate that these names are not 100% certain but rather a “best guess.” Interestingly, we labeled eight clusters as “new species” because they yielded 16S rRNA gene sequences that showed no hits to our type 16S rRNA gene database, indicating that they were different from all known and validly published species. Each of these clusters consisted of exactly one genome. Finally, for the remaining nine clusters, the available information did not allow us to make a “best guess” nor to consider them novel species with certainty. We labeled those clusters “unidentified species.” An overview of all named genome clusters can be found in Table S3.

We inferred a maximum likelihood species tree of the LGC using one representative genome for each of the 239 genome clusters (Fig. 3). For representative genomes, we took the genomes with the largest number of SCGs, because the type strain genomes were not always of the best quality and because some genome clusters did not contain type strains. In this phylogeny, the eight new species and nine unidentified species were also included. Two of the unidentified species were situated in the Weissella genus; the other new and unidentified species could be found in the Lactobacillus genus. Some of the new/unidentified species had long “tip branches” in the tree, indicating that they were relatively far removed from the other species. We annotated the tree with the phylogroups defined by Zheng et al. (6) and predicted the lifestyle of the new species using the work of Duar et al. (2). Strikingly, the isolation source of each of the eight new species confirmed their predicted lifestyles, whether they were vertebrate-adapted, insect-adapted, or free-living bacterial species (Table 3). For nomadic phylogroups, this was trivially the case, since they can occur in any environment.

Of the 2,459 LGC genomes that passed quality control, 98 were unclassified at the species level in the NCBI database. As a direct result of our genome clustering and cluster naming pipelines, these genomes were automatically assigned to a species: either an existing one or a new or unidentified species (Fig. 4A). The most frequently identified species in this group was Lactobacillus rhamnosus, with 13 previously unclassified genomes assigned to it. Interestingly, several unclassified genomes belonged to clusters that we labeled as new or unidentified species.

Of the genomes that did have a species label in the NCBI database, 74 were reclassified to other species by our approach (Fig. 4B). Of these reclassifications, 21 were a direct consequence of the species mergers and splits we described. Next, 38 reclassifications were clear cases of genomes whose NCBI species label did not correspond to their most closely related type strain. We reclassified 20 genomes from Lactobacillus casei to Lactobacillus paracasei, 5 genomes from Lactobacillus paralimentarius to Lactobacillus kimchii, 3 genomes from Lactobacillus farciminis to Lactobacillus formosensis, 2 genomes from Lactobacillus gallinarum to Lactobacillus helveticus, 2 genomes from Lactobacillus plantarum to Lactobacillus pentosus, 1 genome from Lactobacillus brevis to Lactobacillus hilgardii, 1 genome from Lactobacillus buchneri to Lactobacillus hilgardii, 1 genome from Lactobacillus plantarum to Lactobacillus paraplantarum, and 1 genome from Leuconostoc citreum to Leuc. lactis. Finally, 15 genomes with a species label in the NCBI database ended up in genome clusters considered by us as new or unidentified species.

We constructed a genome-based species taxonomy of the Lactobacillus genus complex by performing single-linkage clustering on 2,459 high-quality genomes with a 94% core nucleotide identity cutoff. We found that the resulting species were discontinuous, fully exclusive, and almost fully transitive. On the basis of a comparison of the de novo species with published species, we proposed nine mergers of two or more species and one split of a species. Further, we discovered at least eight yet-to-be-named species that have not been published before, validly or otherwise. Finally, we have shown that our newly developed method to extract single-copy core genes allows for taxonomic studies on thousands of genomes on a desktop computer and is therefore applicable to other genera or even larger taxa. We believe that correcting the published species in the direction of our de novo taxonomy will lead to more meaningful species that show a consistent diversity.

All genome assemblies classified to the Lactobacillaceae family (NCBI taxonomy id 33958) or Leuconostocaceae family (NCBI taxonomy id 81850) were downloaded from GenBank using the script GenBank_get_genomes_by_taxon.py that is part of the tool pyani, version 0.2.7 (53). Gene prediction was performed on all of these assemblies using Prodigal version 2.6.3 (54).

To be able to rapidly extract single-copy core genes (SCGs) from our large genome data set, we devised a new approach based on seed genomes. Our strategy can be divided into four steps. In the first step, n seed genomes were selected completely at random, and their genes were clustered into gene families. The software program OrthoFinder (55) was used for gene family inference; it performs all-versus-all blastp (56) search, normalizes the scores for sequence length and evolutionary distances between the genomes, applies a score cutoff per gene based on its best bidirectional hits, and finally clusters the resulting weighted graph using MCL (57). This process resulted in gene families for the seed genomes only. In the second step, candidate SCGs were selected from all seed gene families by identifying gene families present in more than k seed genomes. The third step was the identification of these candidate SCGs in all genomes. For this purpose, the seed sequences for each candidate SCG were first aligned using MAFFT (58). Next, profile HMMs were constructed from those alignments using the hmmbuild tool of HMMER (hmmer.org). All genomes were then scanned for occurrences of the candidate SCGs with the hmmsearch tool of HMMER. Finally, a homology score cutoff was determined for each profile HMM (see the next paragraph for the details), and these cutoffs were applied to the HMMER hits to yield all “real” hits of the candidate SCGs in all genomes. In the fourth step, the final SCGs were determined by retaining only candidate SCGs with exactly one gene (copy) in at least p percentage of all genomes.

For the training of profile-specific hmmer score cutoffs, the following strategy was used. First, if a gene was a hit of multiple profiles (multiple candidate SCGs), only the highest scoring profile was kept. Next, profile-specific cutoffs were determined as follows. Imagine we are selecting a score cutoff for a given profile. For each genome that has at least one gene among the hits of the profile, the best-scoring gene is considered a “true” hit. All other genes (second, third, etc., best hits of genomes) are considered “false” hits. A score threshold is then determined that maximizes the F-measure of the hits; this is the harmonic mean of the precision and recall. In effect, this means that the threshold is set in such a way that as many as possible “best hits of genomes” are included in the results, while at the same time as many as possible results should be a “best hit of a genome.” The optimization was performed using the R package ROCR (59). This approach was applied separately to every profile HMM (each one representing a candidate SCG).

The SCG extraction strategy we described here was implemented in progenomics: a general, in-development toolkit for prokaryotic comparative genomics (github.com/SWittouck/progenomics). The parameter values used in this study were 30 for n (number of seed genomes), 25 for k (required number of seed genomes a candidate SCG should be present in), and 95 for p (required percentage of total genomes a final SCG should be present in). The software versions used were progenomics version 0.1.0, OrthoFinder version 2.1.2, blast version 2.6.0, MCL version 14-137, MAFFT version 7.407, HMMER version 3.1b2, and ROCR version 1.0.7.

For each genome, two quality measures were calculated based on the SCGs: the completeness (percentage of SCGs present in the genome) and redundancy (percentage of SCGs having more than one copy in the genome). Only genomes with >90% completeness and <10% redundancy were retained. To compare our quality estimates to those available on the Genome Taxonomy Database (GTDB), we downloaded the metadata of the bacterial genomes of the latest GTDB release on 16 July 2019.

For every unique pair of genomes, two similarity measures were computed: the core nucleotide identity (CNI) and average nucleotide identity (ANI). For the CNIs, multiple alignment was first performed for each SCG on the amino acid level using MAFFT version 7.407 (58) with the default parameters. Those amino acid alignments were then used as a reference to align the SCGs on the nucleotide level (reverse alignment). The aligned nucleotides were then concatenated into one large alignment (“supermatrix”). Next, pairwise nucleotide differences were calculated from this alignment using the distmat tool of EMBOSS version 6.6.0.0 (60), without correcting for multiple substitutions. Those differences were converted to similarities by subtracting them from 1. Pairwise ANIs were calculated using the tool fastANI (9), version 1.1, which quickly computes a relatively accurate approximation of ANI.

Genomes were clustered into de novo species by joining together all genome pairs with a CNI similarity higher than 0.94. One could think of this as “nonhierarchical single-linkage clustering.” A simple algorithm to perform this clustering was implemented in R version 3.5.1 (61).

To identify genomes of type strains (“type genomes”), we first downloaded assembly metadata for all of our genomes from the NCBI database using the rentrez R package version 1.2.1 (62). One of those metadata fields was the strain name of the assembled genome. Next, we compiled a list of all type strain names of validly published species of the six genera that make up the LGC. Since type strains are deposited in at least two but often many more culture collections, all of them have at least two names. To collect those type strain names, we developed an R package called tidytypes (github.com/SWittouck/tidytypes, version 0.1.0), which searches the websites LPSN (15, 16), PNU (17), and StrainInfo (63). StrainInfo was used only to find extra synonyms of type strains of species and subspecies found on LPSN and PNU. By looking up those type strain names (including their many synonyms) in the assembly metadata from NCBI (which includes a strain name field), we were able to identify type genomes in our data set on a large scale and in an automated manner. We then supplemented this list by manually identifying extra type genomes in the following way. We gathered a list of species names that were present in the NCBI assembly metadata but for which we had not already found a type genome. For each of those species, we looked up the type strain names in the original paper that described the species. We then inspected the assembly metadata of the genomes classified to that species on NCBI in an attempt to find type genomes that were missed by the automated approach (for example because the species were not validly published). The result of the automated and manual approaches was a complete list of type genomes in our genome data set.

Two different strategies were employed to assign species names to the genome clusters: one based on the type genomes present in the clusters and one based on comparisons to 16S rRNA gene sequences of type strains. First, we used the list of type genomes to assign species names to genome clusters. If the cluster contained type genomes from exactly one species, we assigned this species name to the cluster. If it contained type genomes from more than one species, we concluded that those species were closely related enough to be considered heterotypic synonyms and assigned the oldest species name to the cluster. If two clusters contained type genomes from different subspecies of the same species, we proposed upgrading the cluster with the nontype subspecies (e.g., Lactobacillus casei subsp. tolerans) to species status (e.g., Lactobacillus tolerans).

For the genome clusters that did not contain one or more type genomes, we predicted 16S rRNA genes using barrnap version 0.9 (64). For all species and subspecies names on LPSN for which we did not find a type genome in our data set, we downloaded a 16S rRNA sequence from GenBank. We then scanned the 16S sequences extracted from our genomes against this reference 16S database using the blastn tool of blast version 2.6.0 (56). We then applied a percentage identity cutoff of 98% on those hits. When a genome cluster contained one or more genomes with a match to an official species through the 16S approach, and the same species was listed in the species labels of one or more genomes, we assigned this species name to the genome cluster as a “best guess.” When a genome cluster yielded 16S genes but none of them showed a match to the reference database, we considered it a new species. Finally, some genome clusters remained that could not be identified or be considered a new species; we labeled those “unidentified.”

For each genome cluster, the genome with the largest number of SCGs present was selected as a representative genome. We then selected the 100 SCGs that had the largest single-copy presence counts in those representative genomes. We then constructed a supermatrix with the selected genomes and SCGs, in the same way as described above in the section “Calculation of pairwise genome similarities” (including reverse alignment). Next, we trimmed columns that had gaps in more than 1% of the genomes using trimal version 1.4.rev15 (65). We inferred a maximum likelihood tree using RAxML version 8.2.11 (66) using the general time reversible model for nucleotide substitutions and the CAT method for modeling mutation rate heterogeneity across columns. We used the “-f a” option of RAxML, which combines rapid bootstrapping (67) with a slow search of the tree space starting from the bootstrap trees. The number of bootstrap trees was determined by the “N autoMRE” option, which will stop looking for bootstrap trees when they do not seem to add any extra information (68).

All processing and visualization of table-format data were done in R version 3.5.1 (61) using the tidyverse set of packages, version 1.2.1 (69). Phylogenetic tree visualization and annotation were done using ggtree version 1.12.7 (70).

NCBI assembly accession numbers for the genomes used in this study can be found in Table S1 in the supplemental material. The code used for this study was split over four repositories. The main pipeline (called “legen,” for Lactobacillus Evolutionary Genomics) is available at https://github.com/SWittouck/legen_pipeline (version 3.1). R Markdown scripts for the data analysis, including the creation of all figures and tables, are available at https://github.com/SWittouck/legen_data_analysis (version 3.0). This repository also contains the necessary data files and can thus be run on its own. Some of the steps of the pipeline, including the extraction of single-copy core genes, were implemented in the progenomics toolkit for prokaryotic comparative genomics for maximal reusability. This toolkit is available at https://github.com/SWittouck/progenomics (version 0.1.0). Finally, searching taxonomic websites to obtain type strain names of published prokaryotic species was implemented in an R package called tidytypes, which can be found at https://github.com/SWittouck/tidytypes (version 0.1.0).