Using Core Genome Alignments To Assign Bacterial Species

Advantages of nucleotide alignments over protein alignments for bacterial species analyses.While core protein alignments are increasingly used for phylogenetics (39, 40), a core nucleotide alignment has more phylogenetically informative positions in the absence of substitution saturation (41), yielding a greater potential for phylogenetic signal. Nucleotide-based analyses outperform amino acid-based analyses at all time scales in terms of resolution, branch support, and congruence with independent evidence (42–44). Nucleotide and protein core genome alignments were constructed for 10 complete Wolbachia genomes (see Table S1 in the supplemental material), and the resulting phylogenies were compared.

There are a number of whole-genome aligners that can be used to generate nucleotide core genome alignments (for an overview, see reference 45). We generated core genome alignments for the 10 Wolbachia genomes using Mauve and Mugsy, currently two of the most commonly used whole-genome aligners (27, 28). Despite the different algorithms used by each of the two programs, the lengths of the generated core genome alignments are similar, with the Mauve and Mugsy core genome alignments being 682,949 and 579,495 bp in length, respectively. However, because Mugsy has been found to be better at handling larger genome data sets (28), Mugsy was used for all subsequent core genome alignment analyses discussed. Despite the presence of large-scale genome rearrangements present between the 10 Wolbachia genomes (Fig. S1), the resulting alignment is 45.9% of the average input genome length, indicating that collinear genomes are not necessarily required to construct core genome alignments. The lack of synteny results merely in more and smaller local colinear blocks.

FIG S1

Synteny between the 10 complete Wolbachia genomes. Syntenies between the 10 complete Wolbachia genomes were compared and visualized using the Artemis Comparative Tool. The red ribbons indicate conserved regions of >3 kbp between two genomes, while the blue ribbons indicate >3-kbp inverted conserved regions. Wolbachia supergroups A (●), B (▲), C (■), D (), E () and F () are represented in the genome subset. Shapes of the same color indicate that the multiple genomes are of the same species as determined using our determined CGASI cutoff of ≥96.8%. Download FIG S1, TIF file, 2.5 MB.

Copyright © 2018 Chung et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

Using the same 10 Wolbachia genomes and the available annotation at NCBI RefSeq, a core protein alignment was generated using FastOrtho, a reimplementation of the tool OrthoMCL (78). Despite the presence of different annotation methods used for gene calls, we were able to generate a protein alignment of 77,868 positions from the concatenated alignments of 152 shared genes between the 10 genomes using the available protein amino acid sequences (Table S1). In total, the core protein alignment contained 16,241 parsimony-informative positions compared to the 124,074 such positions observed with the core nucleotide alignment, indicating a 10-fold increase in potentially informative positions (Fig. 1A). When maximum-likelihood (ML) trees from the core protein alignment and the core nucleotide alignment are compared, the two trees are quite similar in topology and branch length. However, the nodes on the ML tree generated using the core nucleotide alignment consistently have higher bootstrap support values.

• Open in new tab
• Download powerpoint

FIG 1

Comparison of phylogenomic trees generated using the core nucleotide alignments versus protein-based alignments for 10 complete Wolbachia genomes. A maximum-likelihood phylogenomic tree generated from the core genome alignment (CGA) of 10 complete Wolbachia genomes was compared to a core protein alignment (CPA) containing 152 genes present in only one copy (A) and an alignment generated using PhyloPhlAn, containing 176 conserved proteins (B). Wolbachia supergroups A (●), B (▲), C (■), D (), E (), and F () are represented in the genome subset. Shapes of the same color indicate that the multiple genomes are of the same species as defined using our determined CGASI cutoff of ≥96.8%. When comparing the ML trees generated using the core nucleotide and core protein alignments, the trees are largely similar in both topology and branch length. In the comparison between the ML trees generated using the core nucleotide alignment and PhyloPhlAn, the trees are similar except for the relationship of wRi, which is sister to wHa in the core nucleotide alignment tree and sister to wAu + wMel in the PhyloPhlAn tree. Despite differences in clustering, the core nucleotide alignment ML tree consistently has higher bootstrap values than its core protein alignment or PhyloPhlAn counterpart.

Similarly, we compared the Wolbachia core nucleotide alignment to a protein-based alignment generated using PhyloPhlAn (39). PhyloPhlAn uses a set of 400 of the most conserved proteins across diverse bacterial taxa to infer phylogenetic relationships. Across the 10 complete Wolbachia genomes, 176 of these 400 genes were identified to be present in each of the genomes. A concatenated amino acid alignment of these 176 genes yields an alignment with 2,877 positions, of which 2,447 are parsimony informative, indicating a 50-fold decrease in the number of parsimony-informative positions from that of the core nucleotide alignment and an 8-fold decrease relative to that of the core protein alignment. As with the core protein alignment, the ML trees constructed using the Wolbachia core nucleotide alignment and PhyloPhlAn have similar topologies and relative branch lengths, with the core nucleotide alignment having higher support values on its nodes (Fig. 1B). The one difference in clustering occurs within the supergroup A Wolbachia, in which wRi clusters with wHa in the core nucleotide alignment, while wRi clusters with wMel and wAu in the PhyloPhlAn ML tree. The wRi and wHa branch is supported by a bootstrap value of 100 in the core nucleotide alignment ML tree, while the wRi, wMel, and wAu branch is supported by a bootstrap value of 98.9 in the PhyloPhlAn ML tree.

Collectively, these results demonstrate that while nucleotide and protein alignments give roughly the same result, there are key differences. Core nucleotide alignments have almost an order of magnitude more parsimony-informative positions, and this results in stronger branch support values in the corresponding phylogenies.

Can core genome alignments be constructed from bacteria within a genus?To compare differences between existing classically defined species designations by this method of using sequence identity thresholds with whole-genome phylogenies, core genome alignments were constructed for seven bacterial genera representative of six different bacterial taxonomic classes (Table 1). For these seven genera, which include Arcobacter, Caulobacter, Erwinia, Neisseria, Polaribacter, Ralstonia, and Thermus, each of the generated core genome alignments were found to be ≥0.33 Mbp in size, representative of ≥13.6% of the average input genome size (Table 1; Table S2). For this diverse group of genera, the core genome alignments comprise a large portion of their input genomes, indicating that this technique is applicable to a wide range of classically defined bacterial taxa. Additionally, core genome alignments were successfully constructed from the genomes of four of the six Rickettsiales genera, including 69 Rickettsia genomes, 3 Orientia genomes, 16 Ehrlichia genomes, and 23 Wolbachia genomes (Table 1). Each of these four core genome alignments are >0.18 Mbp in length and contain >10% of the average size of the input genomes.

Substitution saturation.A common caveat in nucleotide alignments is the presence of substitution saturation, an issue in which the phylogenetic distances reported by nucleotide alignments are underestimated due to an overabundance of reverse mutations within the organism subset (41). Compared to protein-based methods, substitution saturation more heavily impacts nucleotide-based phylogenetic distance measurements due to the higher rate of substitutions per position. However, for each core genome alignment, when the uncorrected pairwise genetic distances are plotted against the model-corrected distances, linear relationships are observed for all alignments (r² > 0.995), indicating that substitution saturation does not hinder the ability of the core genome alignments to represent evolutionary relationships (46) (Fig. S2). Given the fact that the analyzed genera span bacterial taxonomic classes, including Alphaproteobacteria, Betaproteobacteria, Gammaproteobacteria, Epsilonproteobacteria, Flavobacteriia, and Deinococci (Table 1), we expect this result to be broadly applicable to bacterial core genome alignments.

Assessment of genus designations.We found our initial core genome alignments for Anaplasma and Neorickettsia to be considerably shorter than other core genome alignments, both being ∼20 kbp and accounting for ≤2.3% of the average input genome sizes. The small size of the core genome alignment indicates that the Anaplasma and Neorickettsia genome subsets each contain a diverse set of genomes that are likely not of the same genus. Using input genomes from different genera to construct a core genome alignment yields an alignment of an insufficient size to accurately represent the evolutionary distances between the input genomes. For example, when the genome of Neoehrlichia lotoris is supplemented with the genomes used to create the 0.49-Mbp Ehrlichia core genome alignment, the resultant core genome alignment is 0.11 Mbp, representing only 8.9% of the average input genome size, compared to the prior 39.8%. Therefore, we used subsets of species to test whether the Anaplasma and Neorickettsia genera are too broadly defined. A core genome alignment generated using only the 20 A. phagocytophilum genomes in the 30 Anaplasma genome set is 1.25 Mbp and represents 83.3% of the average input genome size, while a core genome alignment of the remaining 10 Anaplasma genomes is 0.77 Mbp, 65.3% of the average genome input size. This result suggests that the Anaplasma genus should be split into two separate genera. Similarly, when the genome of Neorickettsia helminthoeca Oregon is excluded from the Neorickettsia core genome alignment, a 0.77-Mbp Neorickettsia core genome alignment is generated, 87.4% of the average input genome size, suggesting that N. helminthoeca Oregon is not of the same genus as the other three Neorickettsia genomes. For the remainder of this paper, these genus reclassifications are used. Given these collective results, we recommend that the genus classification level can be defined as a group of genomes that together will yield a core genome alignment that represents ≥10% of the average input genome sizes.

Identifying a CGASI for species delineation.Within the Rickettsia, Orientia, Ehrlichia, Anaplasma, Neorickettsia, Wolbachia, Arcobacter, Caulobacter, Erwinia, Neisseria, Polaribacter, Ralstonia, and Thermus genome subsets, ANI, dDDH, and core genome alignment sequence identity (CGASI) values were calculated for 7,264 intragenus pairwise genome comparisons (Table S3), of which 601 are between members of the same species. The ANI and CGASI follow a second-degree polynomial relationship (r² = 0.977) (Fig. 2A). Using this model, the ANI species cutoff of ≥95% is analogous to a CGASI cutoff of ≥96.8%. dDDH and CGASI follow a third-degree polynomial model (r² = 0.978), with a dDDH of 70% being equivalent to a CGASI of 97.6%, indicating that the dDDH species cutoff is generally more stringent than the ANI species cutoff (Fig. 2B).

• Open in new tab
• Download powerpoint

FIG 2

ANI, dDDH, and CGASI correlation analysis. CGASI, ANI, and dDDH values were calculated for 7,264 intragenus pairwise comparisons of genomes for Rickettsia, Orientia, Ehrlichia, Anaplasma, Neorickettsia, Wolbachia, Caulobacter, Erwinia, Neisseria, Polaribacter, Ralstonia, and Thermus. (A) The CGASI and ANI values for the intragenus comparisons follow a second-degree polynomial model (r² = 0.977), with the ANI species cutoff of ≥95% being equivalent to a CGASI of 96.8%, indicated by the blue dashed box. (B) The CGASI and dDDH values for all pairwise comparisons follow a third-degree polynomial model (r² = 0.978), with the dDDH species cutoff of ≥70% being equivalent to a CGASI of 97.6%, indicated by the red dashed box. (C) To identify the optimal CGASI cutoff to use when classifying species, for each increment of the CGASI species cutoff plotted on the x axis, the percentage of intraspecies and interspecies comparisons correctly assigned was determined based on classically defined species designations. The ideal cutoff should maximize the prediction of classically defined species for both interspecies and intraspecies comparisons. The ANI-equivalent CGASI species cutoff is represented by the blue dashed line, while the dDDH-equivalent CGASI species cutoff is represented by the red dashed line.

The ideal CGASI threshold for species delineation would maximize the prediction of classically defined species, neither creating nor destroying the majority of the classically defined species. To examine this, all 7,264 intragenus pairwise comparisons were classified as either intraspecies or interspecies. Intraspecies comparisons are comparisons between genomes with the same classically defined species designation, while interspecies comparisons are between genomes within the same genus but with different classically defined species designations. Every possible CGASI threshold value was then tested for the ability to recapitulate these classically defined taxonomic classifications (Fig. 2C). In all cases, an abnormally high number of interspecies Rickettsia comparisons were found above both the established ANI and the dDDH species thresholds, consistent with previous observations that guidelines for establishing novel Rickettsia species are too lax (47), and as such they were excluded from this specific analysis. Below a CGASI of 97%, classically defined species begin to be separated, while organisms classically defined as different species begin to be collapsed. This coincides with the above-calculated ANI-equivalent threshold but differs from the above-calculated dDDH equivalent (Fig. 2C). The dDDH-equivalent CGASI threshold of 97.6% failed to predict the classically defined taxa from 100 intraspecies comparisons (Table S3) while the ANI-equivalent CGASI threshold of 96.8% failed to predict the classically defined taxa from 41 intraspecies comparisons (Table S3).

Given these results, we selected the ANI-equivalent CGASI value of ≥96.8% to further analyze these taxa. Overall, there were 41 pairwise comparisons of organisms across diverse taxa, including Arcobacter, Caulobacter, Neisseria, Orientia, Ralstonia, and Thermus, that are classically defined as the same species but would be classified as distinct species when assessed using our suggested CGASI threshold (Table S3). Additionally, there were 782 pairwise comparisons of organisms classically defined as different species that this threshold suggests should be the same species, of which only 10 were not in the genus Rickettsia, instead belonging to the Arcobacter and Caulobacter genera (Table S3).

Rickettsiaceae phylogenomic analyses.(i) Rickettsia. The Rickettsiaceae family includes two genera, the Rickettsia and the Orientia, and while both genera are obligate intracellular bacteria, Rickettsia genomes have undergone more reductive evolution, having a genome size ranging from 1.1 to 1.5 Mbp (48) compared to the 2.0- to 2.2-Mbp size of the Orientia genomes (49). Of the Rickettsiales, the Rickettsia genus contains the greatest number of sequenced genomes and named species, containing 69 genome assemblies in ≤100 contigs representing 27 unique species. Rickettsia genomes are currently classified based on the Fournier criteria, an MLST approach established in 2003 based on the sequence similarity of five conserved genes: the 16S rRNA gene, citrate synthase gene (gltA), and three surface-exposed protein antigen genes (ompA, ompB, and gene D) (50). To be considered a Rickettsia species, an isolate must have a sequence identity of ≥98.1% to the 16S rRNA and ≥86.5% to gltA in at least one preexisting Rickettsia species. Within the Rickettsia, using ompA, ompB, and gene D sequence identities, the Fournier criteria also support the further classification of Rickettsia species into three groups: the typhus group, the spotted fever group (SFG), and the ancestral group (50). However, the Fournier criteria have not yet been amended to classify the more recently established transitional Rickettsia group (51), indicating a need to update the Rickettsia taxonomic scheme.

A total of 69 Rickettsia genomes representative of 27 different established species were used for ANI, dDDH, and core genome alignment analyses. Regardless of the method used, a major reclassification is justified (Fig. 3; Table S2). A core genome alignment constructed using the 69 Rickettsia species genomes yielded a core genome alignment size of ∼0.56 Mbp, 42.4% of the average lengths of the input Rickettsia genomes (Table 1). For 42 of the 44 SFG Rickettsia genomes, the CGASI between any two genomes is ≥98.2%, well within the proposed CGASI species cutoff of ≥96.8% (Fig. 3), while the CGASI is ≤97.2% in the ancestral and typhus groups. If a CGASI cutoff of 96.8% is used to reclassify the Rickettsia species, all but two of the SFG Rickettsia genomes would be classified as the same species (Fig. 3), with the two remaining SFG Rickettsia genomes, Rickettsia monacensis IrR Munich and Rickettsia sp. strain Humboldt, being designated as the same species. This is consistent with ANI results as well, while dDDH yields conflicting results (Fig. 3). For the transitional group Rickettsia, Rickettsia akari and Rickettsia australis would be collapsed into a single species due to having CGASI values of 97.2% in a comparison with one another. Similarly, Rickettsia asembonensis, Rickettsia felis, and Rickettsia hoogstraalii, all classified as transitional group Rickettsia, would be collapsed into another species, all having CGASI values of 97.2% in comparisons with one another. This is consistent with a phylogenomic tree generated from the Rickettsia core genome alignments, where the SFG Rickettsia genomes have far less sequence divergence than the rest of the Rickettsia genomes (Fig. 3).

• Open in new tab
• Download powerpoint

FIG 3

Analysis of the ANI, dDDH, and CGASI values of 69 Rickettsia genomes. Across 69 Rickettsia genomes, the ANI and dDDH values (A) and CGASI values (B) were calculated for each pairwise genome comparison. The shape next to each Rickettsia genome represents whether the genome originates from an ancestral (⚫), transitional (), typhus group (▲), or spotted fever group (■) Rickettsia species, while the colors of the shapes on the axes represent species designations as determined by a CGASI cutoff of ≥96.8%. (C) An ML phylogenomic tree with 1,000 bootstraps was generated using the core genome alignment. (D) The relationships in the green box in panel C cannot be adequately visualized at the necessary scale, so they are illustrated separately with a different scale. For both trees, red branches represent branches with <100 bootstrap support.

(ii) Orientia.The organisms within Orientia have no standardized criteria to define novel species. There are far fewer high-quality Orientia genomes than Rickettsia genomes, which is due partly to the large number of repeat elements found in Orientia genomes, with the genome of Orientia tsutsugamushi being the most highly repetitive sequenced bacterial genome to date, with ∼42% of its genome being comprised of short repetitive sequences and transposable elements (52).

The Orientia core genome alignment was constructed using three O. tsutsugamushi genomes and is 0.97 Mbp in size, ∼47.6% of the average input genome size, with CGASI values ranging from 96.3 to 97.2% (Fig. S3). Reclassifying the Orientia genomes using a CGASI species cutoff of 96.8% would result in O. tsutsugamushi Gillliam and O. tsutsugamushi Ikeda being classified as species separate from O. tsutsugamushi Boryong (Fig. S3). In this case, this reclassification would not be consistent with recommendations from using ANI or dDDH. We suspect that ANI is strongly influenced by the large number of repeats in the genome due to ANI calculations being based off the sequence identity of 1-kbp query genome fragments. In comparison, we do not anticipate that whole-genome alignments would be confounded by the repeats. While the LCBs may be fragmented by the repeats, creating smaller syntenic blocks, the nonphylogenetically informative repeats are eliminated from an LCB-based analysis.

Anaplasmataceae phylogenomic analyses.(i) Ehrlichia. Within the Anaplasmataceae, species designations are frequently assigned based on sequence similarity and clustering patterns from phylogenetic analyses generated using the sequences of genes such as 16S rRNA, groEL, and gltA (53). For example, the species designations for Ehrlichia khabarensis and Ehrlichia ornithorhynchi are justified based on having a lower sequence similarity for 16S rRNA, groEL, and gltA (53–55).

A core genome alignment constructed using 16 Ehrlichia genomes, representative of four defined species, yields a 0.49-Mbp alignment, which equates to 39.8% of the average Ehrlichia genome size (1.25 Mbp) (Table 1). Using a CGASI species cutoff of 96.8%, the Ehrlichia chaffeensis and Ehrlichia ruminantium genomes were recovered as monophyletic and well-supported species, which is consistent with ANI and dDDH results (Fig. 4). The two Ehrlichia muris and Ehrlichia sp. strain Wisconsin_h genomes have CGASI values of >97.8%, indicating that the three genomes represent one species, which is consistent with ANI but not dDDH results (Fig. 4). The genomes of Ehrlichia sp. strain HF and E. canis Jake do not have CGASI values of ≥96.8% with any other species, confirming their status as individual species, identical to the results found with ANI and dDDH (Fig. 4).

• Open in new tab
• Download powerpoint

FIG 4

Analysis of the ANI, dDDH, and CGASI values of 16 Ehrlichia genomes. For 16 Ehrlichia genomes, the ANI and dDDH values (A) and CGASI values (B) were calculated for each pairwise genome comparison. The colors of the shapes next to each Ehrlichia genome represent species designations as determined by a CGASI cutoff of ≥96.8%. (C) An ML phylogenomic tree with 1,000 bootstraps was generated using the core genome, with red branches representing branches with <100 bootstrap support.

(ii) Anaplasma.A novel Anaplasma species is currently defined based on phylogenetic analyses involving 16S rRNA, gltA, and groEL, with a new species having a lower sequence identity and a divergent phylogenetic position relative to established Anaplasma species (56–59). A core genome alignment constructed using 30 Anaplasma genomes yielded a 20-kbp alignment, 1.4% of the average input genome size. As noted before, such low values are indicative of more than one genus being represented in the taxa that were included in the analysis. Thus, CGASI analyses for Anaplasma species were done on two Anaplasma genome subsets, one containing the 20 Anaplasma phagocytophilum genomes and the other containing 9 Anaplasma marginale genomes and 1 Anaplasma centrale genome (Table 1).

A 1.25-Mbp core genome alignment, consisting of 39.8% of the average input genome size, was constructed using the 20 A. phagocytophilum genomes. All 20 genomes have CGASI values of ≥96.8%, supporting their designation as members of a single species (Fig. 5). This is supported by ANI, but dDDH again yields conflicting results (Fig. 5). The core genome alignment generated using the remaining 10 Anaplasma genomes yields a 0.77-Mbp core genome alignment, 65.3% of the average input genome size. While the A. marginale genomes all have CGASI values of ≥96.8% in comparisons with one another, the Anaplasma centrale genome has CGASI values ranging from 90.7 to 91.0% compared to the nine A. marginale genomes, supporting A. centrale as a separate species from A. marginale, consistent with existing taxonomy and with the ANI and dDDH species cutoffs (Fig. 5).

• Open in new tab
• Download powerpoint

FIG 5

Analysis of the ANI, dDDH, and CGASI values of 30 Anaplasma genomes. (A) For 30 Anaplasma genomes, the ANI and dDDH values were calculated for each genome comparison and color-coded to illustrate the results with respect to ANI cutoffs of ≥95% and dDDH cutoffs of ≥70%. (B) When we attempted to construct a core genome alignment using all 30 Anaplasma genomes, only a 20-kbp alignment was generated, accounting for <1% of the average Anaplasma genome size. Therefore, CGASI values were calculated after the Anaplasma genomes were split into two subsets containing 20 A. phagocytophilum genomes and the remaining 10 Anaplasma genomes. In all panels, the shape next to each genome denotes genus designations as determined by the size of their core genome alignments, while the color of the shape denotes the species as defined by a CGASI cutoff of ≥96.8%.

(iii) Neorickettsia.The Neorickettsia genus contains four genome assemblies: Neorickettsia helminthoeca Oregon, Neorickettsia risticii Illinois, Neorickettsia sennetsu Miyayama, and Neorickettsia sp. strain 179522. The genus was first established in 1954 with the discovery of N. helminthoeca (60). In 2001, N. risticii and N. sennetsu, both initially classified as Ehrlichia strains, were added to the Neorickettsia based on phylogenetic analyses of 16S rRNA and groESL (31). A core genome alignment constructed using all four Neorickettsia genomes yields a 20-kbp alignment, 2.3% of the average input genome size. When excluding the genome of N. helminthoeca Oregon, the three remaining Neorickettsia genomes form a core genome alignment of 0.76 Mbp in size, 87.4% of the average input genome size. The sequence identity of the core genome alignment indicates that N. risticii Illinois, N. sennetsu Miyayama, and Neorickettsia sp. 179522 are three distinct species within the same genus, while the length of the core genome alignment indicates that N. helminthoeca Oregon is of a separate genus (Fig. S4). When assessing the Neorickettsia species designations using ANI and dDDH cutoffs, the four Neorickettsia genomes can only be determined to be different species, as the two techniques are unable to delineate phylogenomic relationships at the genus level.

(iv) Wolbachia.The current Wolbachia classification system lacks traditional species designations and instead groups organisms by supergroup designations using an MLST system consisting of 450- to 500-bp internal fragments of five genes: gatB, coxA, hcpA, ftsZ, and fbpA (61). A core genome alignment generated using 23 Wolbachia genomes yields a 0.18-Mbp alignment, amounting to 14.8% of the average input genome size.

Among filarial Wolbachia supergroups C and D, the CGASI cutoff of 96.8% would split each of the traditionally recognized supergroups into two groups each, which is also supported by ANI and dDDH. While wOo and wOv would be the same species, wDi Pavia should be considered a different species. Similarly, wBm and wWb would be considered the same species, while wLs should be designated a separate species. The wCle, wFol, and wPpe endosymbionts from supergroups E, F, and L, respectively, would all be considered distinct species using CGASI, ANI, or dDDH.

The CGASI values between the supergroup A Wolbachia organisms, apart from wInc SM, have CGASI values of ≥96.8 (Fig. 6). The genome of wInc SM has CGASI values ranging from 94.6% to 95.9% compared to other supergroup A Wolbachia genomes, indicating that if species assignments were dependent solely on a CGASI cutoff, then wInc SM would be a different species. Because >10% of the positions in the wInc SM genome are ambiguous nucleotide positions, this leads to a large penalty in sequence identity scores, as seen with CGASI, ANI, and DDH (Fig. 6). However, a phylogenomic tree generated using the Wolbachia core genome alignment indicates that wInc SM is nested within the other supergroup A Wolbachia taxa, a clade with 100% bootstrap support (Fig. 6). Similarly, the results between CGASI, ANI, and dDDH in supergroup B are discordant. Five of the six supergroup B Wolbachia genomes, all except for wTpre, have CGASI values of ≥96.8% compared to one of the other five, indicating the five genomes are likely of the same species. However, despite wTpre being considered a different species if the CGASI cutoff of 96.8% is used, a phylogenomic tree constructed using the core genome alignment shows wTpre to be nested within the supergroup B Wolbachia organisms (Fig. 6).

• Open in new tab
• Download powerpoint

FIG 6

Analysis of the ANI, dDDH, and CGASI values of 23 Wolbachia genomes. For 23 Wolbachia genomes, the ANI and dDDH values (A) and CGASI values (B) were calculated for each pairwise genome comparison. The shape next to each Wolbachia genome represents supergroup A (⚫), B (▲), C (■), D (), E (), F (), and L (◆) designations, while the color of the shape indicates species as defined by a CGASI cutoff of ≥96.8%. (C) An ML phylogenomic tree was generated using the Wolbachia core genome alignment constructed using 23 Wolbachia genomes, with red branches representing branches with <100 bootstrap support. The species designations in supergroup B show the CGASI-designated species clusters of wPip_Pel, wPip_JHB, wAus, and wStri to be polyphyletic.

In both cases, if wInc SM or wTpre were designated distinct species from the other supergroup A and B Wolbachia, respectively, paraphyletic clades would be created. This indicates the importance of using a phylogenomic analysis in addition to a threshold. In this case, the core genome alignment can easily be used with robust phylogenetic algorithms to generate both network trees and ML trees. For both cases, despite wInc SM and wTpre having <96.8% CGASI, the phylogenies reveal that these organisms should be considered the same species as the other supergroup A and B Wolbachia species, respectively.

Other taxa.The power of an approach that combines phylogenetics with a nucleotide identity threshold is similarly observed with Neisseria. A 0.33-Mbp core genome alignment was generated from 66 Neisseria genomes, accounting for 14.7% of the average input genome size. The core genome alignment largely supports the classically defined species, including Neisseria meningitidis, Neisseria gonorrhoeae, Neisseria weaveri, and Neisseria lactamica, although one N. lactamica isolate, N. lactamica 338.rep1_NLAC, appears to be inaccurately assigned (Fig. S5). The CGASI values suggest that each Neisseria elongata isolate is a distinct species, as are the Neisseria flavescens isolates. Like what was observed in Wolbachia species, a paraphyletic clade is observed when using a CGASI cutoff of 96.8%, with the genome of Neisseria sp. strain HMSC061B04 being nested within a clade of two Neisseria mucosa genomes, N. lactamica 338.rep1_NLAC, and several unnamed Neisseria taxa while not having a CGASI of ≥96.8% compared to any other Neisseria genome (Fig. S5).