INTRODUCTION {#s1}
============

While acknowledging the disdain that some scientists have for taxonomy, Stephen Jay Gould frequently highlighted in his writings how classifications arising from a good taxonomy both reflect and direct our thinking, stating, "the way we order reflects the way we think. Historical changes in classification are the fossilized indicators of conceptual revolutions" ([@B1]). Historically, bacterial species delimitation relied on phenotypic, morphologic, and chemotaxonomic characterizations ([@B2][@B3][@B4]). The 1960s saw the introduction of molecular techniques in bacterial species delimitation through the use of GC content ([@B5]), DNA-DNA hybridization (DDH) ([@B6]), and 16S rRNA sequencing ([@B7], [@B8]). Currently, databases like SILVA ([@B9]) and Greengenes ([@B10]) use 16S rRNA sequencing to identify bacteria. However, 16S rRNA sequencing often fails to separate closely related taxa, and its utility for species-level identification is questionable ([@B10][@B11][@B12]). Multilocus sequence analysis (MLSA) has also been used to delineate species ([@B13]), as has the phylogenetic analysis of both rRNA and protein-coding genes ([@B3], [@B14], [@B15]). Nongenomic mass spectrometry-based approaches, in which expressed proteins and peptides are characterized, provide complementary data to phenotypic and genomic species delimitations ([@B16], [@B17]) and are used in clinical microbiology laboratories. However, DDH remains the "gold standard" of defining bacterial species ([@B18], [@B19]), despite the intensive labor involved and its inability to be applied to nonculturable organisms. A new genome-based bacterial species definition is attractive, given the increasing availability of bacterial genomes, rapid sequencing improvements with decreasing sequencing costs, and data standards and databases that enable data sharing.

Average nucleotide identity (ANI) and digital DDH (dDDH) were developed as genomic-era tools that allow for bacterial species classification with a high correlation to results obtained using wet lab DDH, while bypassing the associated difficulties ([@B20][@B21][@B22]). For ANI calculations, the genome of the query organism is split into 1-kbp fragments, which are then searched against the whole genome of a reference organism. The average sequence identity of all matches having \>60% overall sequence identity over \>70% of their length is defined as the ANI between the two organisms ([@B18]). dDDH uses the sequence similarity of conserved regions between two genomes of interest ([@B23]) to calculate genome-to-genome distances. These distances are converted to a dDDH value, which is intended to be analogous to DDH values obtained using traditional laboratory methods ([@B23]). There are three formulas for calculating dDDH values between two genomes using either (i) the length of all matches divided by the total genome length, (ii) the sum of all identities found in matches divided by the overall match length, and (iii) the sum of all identities found in matches divided by the total genome length, with the second formula being recommended for assigning species designations for draft genomes ([@B24], [@B25]). However, neither ANI nor dDDH reports the total length of fragments that match the reference genome, and problems arise when only a small number of fragments are unknowingly used. Additionally, neither method generates an alignment that can be used for complementary phylogenetic analyses, instead relying solely on strict cutoffs to define species, potentially leading to the formation of paraphyletic taxa.

Some recent phylogenomic analyses have shifted toward using the core proteome, a concatenated alignment constructed using the amino acid sequences of genes shared between the organisms of interest ([@B26]). However, differences in annotation that affect gene calls can add an unnecessary variable when deriving evolutionary relationships. Instead, we propose that such analyses should use nucleotide core genome alignments to infer phylogenetic relationships with a well-defined nucleotide identity threshold to identify where trees should be pruned to define species. Such an analysis would be enabled by genome aligners like Mauve ([@B27]) and Mugsy ([@B28]), which identify regions of nucleotide identity that are shared between genomes, which are referred to as locally collinear blocks (LCBs). For any given subset of genomes, a core genome alignment can be generated by concatenating these LCBs and retaining only positions present in all genomes. Using the length and sequence identity of the core genome alignment, paired with a phylogeny generated from this alignment, genus- and species-level taxonomic assignments can be developed. A core genome alignment-based method provides advantages to its protein-based counterpart in that it is of a higher resolution and independent of annotation, while also being transparent with respect to the data used in the calculations and very amenable to data sharing and deposition in data repositories.

The *Rickettsiales* are an order within the *Alphaproteobacteria* composed of obligate, intracellular bacteria where classic DNA-DNA hybridization is not possible and bacterial taxonomy is uneven, with each of the genera having its own criteria for assigning genus and species designations. Within the *Rickettsiales*, there are three major families, the *Anaplasmataceae*, *Midichloriaceae*, and *Rickettsiaceae*, with an abundance of genomic data being available for genera within the *Anaplasmataceae* and *Rickettsiaceae.* However, as noted above, species definitions in these families are inconsistent, and organisms in the *Wolbachia* genus lack community-supported species designations. While the *Wolbachia* community has vigorously discussed nomenclature at each of its biannual meetings over the past 20 years, thus far, the community has supported only supergroup designations. Proposed ANI- and dDDH-informed *Wolbachia* species definitions and nomenclature ([@B29]) lack support from the *Wolbachia* community for a number of reasons (see, e.g., reference [@B30]).

The last reorganization of the *Rickettsiaceae* taxonomy occurred in 2001 ([@B31]), a time when there were \<300 sequenced bacterial genomes ([@B32]), a limited number of *Rickettsia* genomes ([@B33], [@B34]), and no *Anaplasmataceae* genomes ([@B35][@B36][@B38]). As of 2014, there were \>14,000 bacterial genomes sequenced, and with this increase in available genomic information, more-informed decisions can be made with regard to taxonomic classification ([@B32]). By using core genome alignments, we can condense whole genomes into positions shared between the input genomes and use sequence identity to infer phylogenomic relationships. We can identify sequence identity thresholds for these alignments that are consistent with classically defined bacterial species and can be overlaid on phylogenetic trees generated from the same alignment. This combined approach, which we present in this study, enables a phylogenetically informed whole-genome approach to bacterial taxonomy. We apply this approach to organisms in the *Rickettsiales*, including *Rickettsia*, *Orientia*, *Ehrlichia*, *Neoehrlichia*, *Anaplasma*, *Neorickettsia*, and *Wolbachia*.

RESULTS {#s2}
=======

Advantages of nucleotide alignments over protein alignments for bacterial species analyses. {#s2.1}
-------------------------------------------------------------------------------------------

While core protein alignments are increasingly used for phylogenetics ([@B39], [@B40]), a core nucleotide alignment has more phylogenetically informative positions in the absence of substitution saturation ([@B41]), 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 ([@B42][@B43][@B44]). Nucleotide and protein core genome alignments were constructed for 10 complete *Wolbachia* genomes (see [Table S1](#tabS1){ref-type="supplementary-material"} in the supplemental material), and the resulting phylogenies were compared.

10.1128/mSystems.00236-18.7

Ten complete *Wolbachia* genomes used for PhyloPhlAn and nucleotide and protein core genome analyses. Download Table S1, XLSX file, 0.01 MB.

Copyright © 2018 Chung et al.

2018

Chung et al.

This content is distributed under the terms of the

Creative Commons Attribution 4.0 International license

.

There are a number of whole-genome aligners that can be used to generate nucleotide core genome alignments (for an overview, see reference [@B45]). We generated core genome alignments for the 10 *Wolbachia* genomes using Mauve and Mugsy, currently two of the most commonly used whole-genome aligners ([@B27], [@B28]). 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 ([@B28]), 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](#figS1){ref-type="supplementary-material"}), 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.

10.1128/mSystems.00236-18.1

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 (![](sys0061822900008.jpg)), E (![](sys0061822900009.jpg)) and F (![](sys0061822900010.jpg)) 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.

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 ([@B78]). 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](#tabS1){ref-type="supplementary-material"}). 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](#fig1){ref-type="fig"}). 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.

![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 (![](sys0061822900008.jpg)), E (![](sys0061822900009.jpg)), and F (![](sys0061822900010.jpg)) 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 *w*Ri, which is sister to *w*Ha in the core nucleotide alignment tree and sister to *w*Au + *w*Mel 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.](sys0061822900001){#fig1}

Similarly, we compared the *Wolbachia* core nucleotide alignment to a protein-based alignment generated using PhyloPhlAn ([@B39]). 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](#fig1){ref-type="fig"}). The one difference in clustering occurs within the supergroup A *Wolbachia*, in which *w*Ri clusters with *w*Ha in the core nucleotide alignment, while *w*Ri clusters with *w*Mel and *w*Au in the PhyloPhlAn ML tree. The *w*Ri and *w*Ha branch is supported by a bootstrap value of 100 in the core nucleotide alignment ML tree, while the *w*Ri, *w*Mel, and *w*Au 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? {#s2.2}
-----------------------------------------------------------------------

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](#tab1){ref-type="table"}). 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](#tab1){ref-type="table"}; [Table S2](#tabS2){ref-type="supplementary-material"}). 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](#tab1){ref-type="table"}). Each of these four core genome alignments are \>0.18 Mbp in length and contain \>10% of the average size of the input genomes.

###### 

Core genome alignment statistics

  --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Genus or genus and species                     Family                 Class                     No. of\    No. of\        Avg\         Core genome\   \% core\      No. of\      Minimum\
                                                                                                  genomes\   established\   genome\      alignment\     genome\       LCBs\        CGASI (%)
                                                                                                  analyzed   species\       size (Mbp)   size (Mbp)     composition   identified   
                                                                                                             represented                                                           
  ---------------------------------------------- ---------------------- ------------------------- ---------- -------------- ------------ -------------- ------------- ------------ -----------
  *Rickettsia*                                   *Rickettsiaceae*       *Alphaproteobacteria*     69         27             1.32         0.56           42.40         36,667       81.70

  *Orientia*                                     *Rickettsiaceae*       *Alphaproteobacteria*     3          1              2.04         0.97           47.50         12,101       96.30

  *Ehrlichia*                                    *Anaplasmataceae*      *Alphaproteobacteria*     16         4              1.23         0.49           39.80         1,423        81.60

  *Neoehrlichia*                                 *Anaplasmataceae*      *Alphaproteobacteria*     1          1              1.27                                                   

  *Neoehrlichia* with *Ehrlichia*                *Anaplasmataceae*      *Alphaproteobacteria*     17         5              1.24         0.11           8.90          1,661        80.10

  *Anaplasma*                                    *Anaplasmataceae*      *Alphaproteobacteria*     30         3              1.39         0.02           1.40          33,010       84.00

  A. phagocytophilum                             *Anaplasmataceae*      *Alphaproteobacteria*     20         1              1.5          1.25           83.30         31,062       98.90

  A. marginale or A. centrale                    *Anaplasmataceae*      *Alphaproteobacteria*     10         2              1.18         0.77           65.30         2,339        90.60

  *Neorickettsia*                                *Anaplasmataceae*      *Alphaproteobacteria*     4          3              0.87         0.02           2.30          55           85.70

  *Neorickettsia* excluding *N. helminthoeca*\   *Anaplasmataceae*      *Alphaproteobacteria*     3          2              0.87         0.76           87.40         31           85.70
  Oregon                                                                                                                                                                           

  *Wolbachia*                                    *Anaplasmataceae*      *Alphaproteobacteria*     23                        1.22         0.18           14.80         14,193       77.20

  *Wolbachia* excluding *w*Ppe                   *Anaplasmataceae*      *Alphaproteobacteria*     22                        1.23         0.54           43.90         14,146       80.10

  *Arcobacter*                                   *Campylobacteraceae*   *Epsilonproteobacteria*   44         12             2.27         0.43           18.90         34,089       78.90

  *Caulobacter*                                  *Caulobacteraceae*     *Alphaproteobacteria*     26         4              4.97         1.61           32.40         80,702       82.50

  *Erwinia*                                      *Enterobacteriaceae*   *Gammaproteobacteria*     22         10             4.42         0.6            13.60         17,290       78.80

  *Neisseria*                                    *Neisseriaceae*        *Betaproteobacteria*      66         10             2.24         0.33           14.70         64,633       80.50

  *Polaribacter*                                 *Flavobacteriaceae*    *Flavobacteriia*          24         11             3.55         0.8            22.50         20,905       77.90

  *Ralstonia*                                    *Burkholderiaceae*     *Betaproteobacteria*      21         3              5.44         1.24           22.80         27,112       82.20

  *Thermus*                                      *Thermaceae*           *Deinococci*              19         11             2.29         1.1            48.00         13,251       80.90
  --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

10.1128/mSystems.00236-18.8

Genomes used for ANI, dDDH, and CGASI analysis. Download Table S2, TXT file, 0.5 MB.

Copyright © 2018 Chung et al.

2018

Chung et al.

This content is distributed under the terms of the

Creative Commons Attribution 4.0 International license

.

Substitution saturation. {#s2.3}
------------------------

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 ([@B41]). 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^2^* \> 0.995), indicating that substitution saturation does not hinder the ability of the core genome alignments to represent evolutionary relationships ([@B46]) ([Fig. S2](#figS2){ref-type="supplementary-material"}). Given the fact that the analyzed genera span bacterial taxonomic classes, including *Alphaproteobacteria*, *Betaproteobacteria*, *Gammaproteobacteria*, *Epsilonproteobacteria*, *Flavobacteriia*, and *Deinococci* ([Table 1](#tab1){ref-type="table"}), we expect this result to be broadly applicable to bacterial core genome alignments.

10.1128/mSystems.00236-18.2

Assessing substitution saturation for core genome alignments. For each of the 14 core genome alignments that comprise ≥10% of the average input genome size, the uncorrected genetic distance between each of the members was plotted against the Tn*69*-model corrected genetic distance. The red line represents the best-fit line for each data set, while the black dotted line represents the identity line (*y* = *x*). In all cases, the relationship between the two distances are linear (*r*^2^ \> 0.995), indicating little substitution saturation in the core genome alignments. Download FIG S2, TIF file, 1.6 MB.

Copyright © 2018 Chung et al.

2018

Chung et al.

This content is distributed under the terms of the

Creative Commons Attribution 4.0 International license

.

Assessment of genus designations. {#s2.4}
---------------------------------

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. {#s2.5}
--------------------------------------------

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](#tabS3){ref-type="supplementary-material"}), of which 601 are between members of the same species. The ANI and CGASI follow a second-degree polynomial relationship (*r^2^* = 0.977) ([Fig. 2A](#fig2){ref-type="fig"}). 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^2^* = 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](#fig2){ref-type="fig"}).

![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^2^* = 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^2^* = 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.](sys0061822900002){#fig2}

10.1128/mSystems.00236-18.9

ANI, dDDH, and CGASI values for 7,264 interspecies comparisons of the genera *Rickettsia*, *Orientia*, *Ehrlichia*, *Anaplasma*, *Neorickettsia*, *Wolbachia*, *Arcobacter*, *Caulobacter*, *Erwinia*, *Neisseria*, *Polaribacter*, *Ralstonia*, and *Thermus*. Download Table S3, TXT file, 0.01 MB.

Copyright © 2018 Chung et al.

2018

Chung et al.

This content is distributed under the terms of the

Creative Commons Attribution 4.0 International license

.

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](#fig2){ref-type="fig"}). 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 ([@B47]), 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](#fig2){ref-type="fig"}). The dDDH-equivalent CGASI threshold of 97.6% failed to predict the classically defined taxa from 100 intraspecies comparisons ([Table S3](#tabS3){ref-type="supplementary-material"}) while the ANI-equivalent CGASI threshold of 96.8% failed to predict the classically defined taxa from 41 intraspecies comparisons ([Table S3](#tabS3){ref-type="supplementary-material"}).

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](#tabS3){ref-type="supplementary-material"}). 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](#tabS3){ref-type="supplementary-material"}).

*Rickettsiaceae* phylogenomic analyses. {#s2.6}
---------------------------------------

**(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 ([@B48]) compared to the 2.0- to 2.2-Mbp size of the *Orientia* genomes ([@B49]). 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) ([@B50]). 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 ([@B50]). However, the Fournier criteria have not yet been amended to classify the more recently established transitional *Rickettsia* group ([@B51]), 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](#fig3){ref-type="fig"}; [Table S2](#tabS2){ref-type="supplementary-material"}). 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](#tab1){ref-type="table"}). 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](#fig3){ref-type="fig"}), 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](#fig3){ref-type="fig"}), 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](#fig3){ref-type="fig"}). 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](#fig3){ref-type="fig"}).

![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 (![](sys0061822900008.jpg)), 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.](sys0061822900003){#fig3}

(ii) *Orientia*. {#s2.7}
----------------

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 ([@B52]).

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](#figS3){ref-type="supplementary-material"}). 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](#figS3){ref-type="supplementary-material"}). 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.

10.1128/mSystems.00236-18.3

Analysis of the ANI, dDDH, and CGASI values of three *Orientia* genomes. For 3 Orientia tsutsugamushi genomes, the ANI and dDDH values (A) and CGASI values (B) were calculated for each genome comparison and color-coded to illustrate the results with respect to cutoffs of an ANI of ≥95% and a dDDH value of ≥70%. Circles of the same colors next to the names of each genome indicate members of the same species as defined by a CGASI of ≥96.8%. Download FIG S3, TIF file, 2.0 MB.

Copyright © 2018 Chung et al.

2018

Chung et al.

This content is distributed under the terms of the

Creative Commons Attribution 4.0 International license

.

*Anaplasmataceae* phylogenomic analyses. {#s2.8}
----------------------------------------

**(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 *glt*A ([@B53]). 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* ([@B53][@B54][@B55]).

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](#tab1){ref-type="table"}). 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](#fig4){ref-type="fig"}). 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](#fig4){ref-type="fig"}). 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](#fig4){ref-type="fig"}).

![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.](sys0061822900004){#fig4}

(ii) *Anaplasma*. {#s2.9}
-----------------

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 ([@B56][@B57][@B59]). 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](#tab1){ref-type="table"}).

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](#fig5){ref-type="fig"}). This is supported by ANI, but dDDH again yields conflicting results ([Fig. 5](#fig5){ref-type="fig"}). 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 Anaplasmacentrale 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](#fig5){ref-type="fig"}).

![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%.](sys0061822900005){#fig5}

(iii) *Neorickettsia*. {#s2.10}
----------------------

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* ([@B60]). 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* ([@B31]). 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](#figS4){ref-type="supplementary-material"}). 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.

10.1128/mSystems.00236-18.4

Analysis of the ANI, dDDH, and CGASI values of 4 *Neorickettsia* genomes. (A) For the 4 *Neorickettsia* genomes, the ANI and dDDH values were calculated for each genome comparison and color-coded to illustrate the results with respect to cutoffs of an ANI of ≥95% and a dDDH value of ≥70%. (B) CGASI values were calculated and are illustrated using a core genome alignment that could only be constructed using 3 of the *Neorickettsia* genomes, excluding *N. helminthoeca* Oregon. Circles of the sample colors next to the names of each genome indicate members of the same species as defined by a CGASI of ≥96.8%. Download FIG S4, TIF file, 1.9 MB.

Copyright © 2018 Chung et al.

2018

Chung et al.

This content is distributed under the terms of the

Creative Commons Attribution 4.0 International license

.

(iv) *Wolbachia*. {#s2.11}
-----------------

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* ([@B61]). 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 *w*Oo and *w*Ov would be the same species, *w*Di Pavia should be considered a different species. Similarly, *w*Bm and *w*Wb would be considered the same species, while *w*Ls should be designated a separate species. The *w*Cle, *w*Fol, and *w*Ppe 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 *w*Inc SM, have CGASI values of ≥96.8 ([Fig. 6](#fig6){ref-type="fig"}). The genome of *w*Inc 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 *w*Inc SM would be a different species. Because \>10% of the positions in the *w*Inc 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](#fig6){ref-type="fig"}). However, a phylogenomic tree generated using the *Wolbachia* core genome alignment indicates that *w*Inc SM is nested within the other supergroup A *Wolbachia* taxa, a clade with 100% bootstrap support ([Fig. 6](#fig6){ref-type="fig"}). Similarly, the results between CGASI, ANI, and dDDH in supergroup B are discordant. Five of the six supergroup B *Wolbachia* genomes, all except for *w*Tpre, 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 *w*Tpre being considered a different species if the CGASI cutoff of 96.8% is used, a phylogenomic tree constructed using the core genome alignment shows *w*Tpre to be nested within the supergroup B *Wolbachia* organisms ([Fig. 6](#fig6){ref-type="fig"}).

![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 (![](sys0061822900008.jpg)), E (![](sys0061822900009.jpg)), F (![](sys0061822900010.jpg)), 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 *w*Pip_Pel, *w*Pip_JHB, *w*Aus, and *w*Stri to be polyphyletic.](sys0061822900006){#fig6}

In both cases, if *w*Inc SM or *w*Tpre 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 *w*Inc SM and *w*Tpre 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. {#s2.12}
-----------

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 Neisseriameningitidis, Neisseriagonorrhoeae, Neisseriaweaveri, and Neisseria lactamica, although one *N. lactamica* isolate, *N. lactamica* 338.rep1_NLAC, appears to be inaccurately assigned ([Fig. S5](#figS5){ref-type="supplementary-material"}). The CGASI values suggest that each Neisseriaelongata isolate is a distinct species, as are the Neisseriaflavescens 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 Neisseriamucosa 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](#figS5){ref-type="supplementary-material"}).

10.1128/mSystems.00236-18.5

Analysis of the ANI, dDDH, and CGASI values of 66 *Neisseria* genomes. For the 66 *Neisseria* genomes, the ANI and dDDH values (A) and CGASI values (B) were calculated for each genome comparison and color-coded to illustrate the results with respect to cutoffs of an ANI of ≥95% and a dDDH value of ≥70%. (C) An ML phylogenomic tree was generated using 1,000 bootstraps, with the 0.33-Mbp *Neisseria* core genome alignment. Circles of the sample colors next to the names of each genome indicate members of the same species as defined by a CGASI of ≥96.8. Download FIG S5, TIF file, 2.8 MB.

Copyright © 2018 Chung et al.

2018

Chung et al.

This content is distributed under the terms of the

Creative Commons Attribution 4.0 International license

.

DISCUSSION {#s3}
==========

Workflow for the taxonomic assignment of new genomes at the genus and species levels. {#s3.1}
-------------------------------------------------------------------------------------

These results highlight that while sequence identity cutoffs derived using nucleotide-based pairwise comparisons are important when delineating species, they should be coupled with an analysis of phylogenetic trees. Yet, phylogenetic trees alone are inadequate, as they cannot identify the level of branching where a species should be defined. However, the two results can be quite complementary when applied to core genome alignments and together make for a robust approach.

Given our results, we propose a workflow using three criteria to guide the taxonomic assignment of novel genomes at the genus and species levels ([Fig. 7](#fig7){ref-type="fig"}). Starting with a novel, sequenced query genome, the most closely related genus of the query genome should be identified using homology-based search algorithms (e.g., BLASTN searches of 16S rRNA). Once identified, a set of trusted genomes within the genus should be combined with this query genome to generate a core genome alignment using a tool such as Mauve ([@B27]) or Mugsy ([@B28]). These trusted genomes would be the type genomes, analogous to type strains. Unlike type strains, they are digitally recorded and thus highly unlikely to be lost. However, multiple organisms with different genome sequences could now serve in the capacity of type strains, eliminating the possibility that a single unusual type strain could unduly influence the taxonomy.

![Workflow for the taxonomic assignment of a novel genome at the genus and species levels. The proposed workflow for assigning genus- and species**-**level designations using core genome alignments is based on three criteria: (i) the length of the core genome alignment, (ii) the sequence identity of the core genome alignment, and (iii) phylogenomic analyses. Using a query genome and a set of trusted genomes from a single genus, a core genome alignment is generated. The length of the core genome alignment is the first criterion that is used as the genus-level cutoff, with a core genome alignment size of ≥10% of the average input genome size indicating that all genomes within the subset are of the same genus. Provided that all genomes in the subset are of the same genus, the second criterion is the sequence identity of the core genome alignment, with genomes sharing ≥96.8% similarity being designated the same species. The third and final criterion uses a phylogenomic tree generated using the core genome alignment to check for paraphyletic clades. If the query genome has \<96.8% core genome alignment sequence identity with any other genome and its designation as a new species does not form a paraphyletic clade, the genome should be considered a novel species.](sys0061822900007){#fig7}

From the core genome alignment, the genus assignment would be considered validated if the length of the core genome alignment is ≥10% of the average input genome size. Subsequently, a pairwise sequence identity matrix is calculated from the core genome alignment, identifying all genomes with ≥96.8% sequence identity. Finally, phylogenomic trees are generated with the core genome alignment and overlaid with the CGASI results. A query genome is deemed a novel species if it is \<96.8% identical to the genomes of the type strains and the clades remain monophyletic. Importantly, whole-genome sequences are not needed to assign an existing genus and species designation to a new isolate, merely to describe new genera or species. Community-derived species-specific MLST schemes may be well suited for assigning isolates to already-named and -described species.

The genomics and data era may necessitate increased flexibility or a modification of the rules in bacterial nomenclature. {#s3.2}
-------------------------------------------------------------------------------------------------------------------------

The International Code of Nomenclature of Prokaryotes ([@B62]) is a frequently updated governance that includes 65 rules that regulate bacterial nomenclature, with the laudable goal of maintaining stability in bacterial names. These rules are designed to assess the correctness of bacterial names, as well as procedures for creating new names, but are agnostic to the methods or criteria used to delineate genera or species ([@B62]).

However, modern methods (i.e., genome sequencing, bioinformatics, data sharing, and isolate repositories) have the potential to profoundly and negatively affect bacterial nomenclature in ways where the current code may need to adjust and adapt. Historically, isolates were frequently discovered, described, and named by the same individual or group of individuals. Conflicts arose over time (e.g., two named species turned out to be the same species), and many of the rules focus on resolving those inevitable conflicts. One principle in the code for resolving conflicts relies on using the first published name. However, with modern methods, it would be fairly simple to conduct a comprehensive bioinformatics analysis and be the first to publish names of organisms, without any input from the relevant research communities, including individuals who have spent their scientific careers working on unnamed organisms. In addition, with extensive genomic data, there have been increased efforts to harmonize taxonomy that will lead to large changes in bacterial nomenclature. While choosing the first published name may be effective at resolving small conflicts, it may not be ideal for larger taxonomic changes that may be necessary and warranted. In these cases, the code may need to prioritize large collaborative community-supported efforts that include all relevant stakeholders to reform the nomenclature in a meaningful and important way, particularly efforts that reflect the recent changes in the ways that we think about bacteria and are reflective of the modern conceptual revolutions to which Stephen Jay Gould referred.

Furthermore, a complete genome sequence fully describes an organism, including its phenotypic potential. As such, a genome-informed taxonomic scheme, particularly one that relies on groups of genome type strains, as described here, may render the "*Candidatus*" designation obsolete. "*Candidatus*" was originally proposed as a category or rank given to a potential new taxon that has only sequence data and a morphology determined through microscopy with a nucleic acid probe ([@B63]). At that time, this precluded the naming of a species with only a single gene sequence, like 16S rRNA, which could be a chimera. While "*Candidatus*" had an important role at that time and prior to the widespread use of whole-genome sequencing, it has led to extensive confusion among researchers, particularly those new to fields with organisms that use this naming scheme. It also places an undue burden on researchers that study obligate host-associated organisms that are often very highly characterized but unamenable to growth on laboratory media and long-term storage in culture collections. Additionally, among endosymbionts and other obligate intracellular organisms, it has been applied inconsistently, compounding this confusion. Therefore, it seems reasonable that its use should be limited to only organisms from which a complete genome sequence cannot be obtained consistently from the same host organism or host cell line.

Limitations. {#s3.3}
------------

A current issue of whole-genome aligners, such as Mauve or Mugsy, is the inability for the software to scale, being able to computationally handle only subsets of at most ∼80 genomes. While the aligner ParSNP in the Harvest suite can generate a core genome alignment from a larger subset of genomes, its use is limited to intraspecies comparisons ([@B45]), making it currently unamenable to this approach. However, for heavily sequenced genera, such as the *Rickettsia* genus, future species assignments for novel genomes do not require constructing a core genome alignment from every available genome. Instead, we recommend that for each genus, the relevant experts in the community establish and curate a set of trusted genomes with at least one representative of each named species that should be used for constructing core genome alignments. After an initial assessment, refinement could be made by using a core genome alignment with many more genomes from closely related species. As an example, in the case of *Rickettsia* species, the initial alignment would include representative genomes from each of the ancestral, typhus, spotted fever, and transitional groups. If the analysis indicated that the query was closest to typhus group genomes, a subsequent alignment would include a curated subset of genomes only from the typhus group. The first core genome alignment would serve to classify a genome into a subgroup, while the second would be used to assign a species designation.

Conclusions. {#s3.4}
------------

As noted with the development of MLST ([@B64]), sequence data have the advantage of being incredibly standardized and portable. While MLST methods allow for insight at the subspecies level of taxonomic classifications and are heavily relied upon during infectious disease outbreaks, they do not provide sufficient resolution to define species. Increased resolution is required for taxonomy and can be obtained through the construction of whole-genome alignments that maximize the number of evolutionarily informative positions.

By generating core genome alignments for different genera, we sought to identify and develop a universal, high-resolution method for species classification. Using the core genome alignment, a set of genomes can be assessed based on the length, sequence identity, and a phylogenomic tree of the same core genome alignment. The length of the core genome alignment, which reflects the ability of the genomes to be aligned, is informative in delineating genera. Meanwhile, the pairwise comparisons of core genome alignment sequence identity can be used to delineate species, which can in turn be further refined with a phylogenetic analysis to resolve paraphyletic clades in the taxonomy.

Through this work, we have identified criteria that reconstruct classical species definitions using a method that is transparent and portable. In the course of this work, we have identified modifications that need to be made to the species and genus designations of a number of organisms, particularly within the *Rickettsiales*. While we have identified these instances, we recommend that any changes in nomenclature be addressed by collaborative teams of experts in the respective communities.

MATERIALS AND METHODS {#s4}
=====================

Synteny plots. {#s4.1}
--------------

For the 10 complete *Wolbachia* genomes, NCBI BLAST v2.2.26 with the blastn algorithm was used to obtain the coordinates of the links for each of the pairwise genome comparisons. Synteny plots were generated and visualized using the Artemis Comparative Tool viewer ([@B65], [@B66]).

Core genome alignments. {#s4.2}
-----------------------

The exact commands for generating a core genome alignment from a directory containing all of the fasta files are provided in a bash script (see [Text S1](#textS1){ref-type="supplementary-material"} in the supplemental material); the local paths for the bin directors for Mugsy and mothur need to be specified in the script. More specifically, genomes used in the taxonomic analyses were downloaded from NCBI's GenBank ([@B67]). OrthoANIu v1.2 ([@B68]) and USEARCH v6.1.544 ([@B69]) were run with default settings and used for ANI calculations. GGDC v2.1 (ggdc.dsmz.de) ([@B25]) paired with the recommended BLAST+ alignment tool ([@B70]) was used to calculate dDDH values. For analyses in this paper, all dDDH calculations were performed using dDDH formula 2 due to the usage of draft genomes in taxonomic analyses ([@B24]). For each of the genome subsets, core genome alignments were generated using Mugsy v1.2 ([@B28]), with default settings, or Mauve v2.4 ([@B28]), with the progressive Mauve algorithm, retaining only LCBs of \>30 bp in length that are present in all genomes. mothur v1.22 was used to process each of the core genome alignments to remove positions not present in all organisms in the genome subset ([@B71]). Sequence identity matrices for the core genome alignments were created using BioEdit v7.2.5 (brownlab.mbio.ncsu.edu/JWB/papers/1999Hall1.pdf). ML phylogenomic trees with 1,000 bootstraps were calculated for each core genome alignment using IQ-TREE v1.6.2 ([@B72]) paired with ModelFinder ([@B73]) to select the best model of evolution and UFBoot2 ([@B74]) for fast bootstrap approximation. Trees were visualized and annotated using iTOL v4.1.1 ([@B75]). Construction of neighbor-network trees was done using the R packages ape ([@B81]) and phangorn ([@B77]).

10.1128/mSystems.00236-18.6

Bash script of commands used to generate a core genome alignment from a directory of whole-genome fasta files. The local paths for Mugsy and mothur must be provided in the bash script. Download Text S1, TXT file, 0.00 MB.

Copyright © 2018 Chung et al.

2018

Chung et al.

This content is distributed under the terms of the

Creative Commons Attribution 4.0 International license

.

Core protein alignments. {#s4.3}
------------------------

Orthologs between complete genomes of the same species were determined using FastOrtho, a reimplementation of OrthoMCL ([@B78]) that identifies orthologs using all-by-all BLAST searches, run with default settings. The amino acid sequences of proteins present in all organisms in only one copy were aligned with MAAFT v7.313 FFT-NS-2, using default settings ([@B79]). For every protein alignment, the best model of evolution was identified using ModelFinder ([@B73]), and phylogenomic trees were constructed using an edge-proportional partition model ([@B80]) with IQ-TREE v1.6.2 ([@B72]) and UFBoot2 ([@B74]) for fast bootstrap approximation. PhyloPhlAn v0.99 was run using default settings as an alternative method of assessing protein-based phylogenies ([@B39]). Comparative analyses of core protein and core nucleotide alignment trees were performed using the R packages ape ([@B76]) and phangorn ([@B77]).

We thank Joseph J. Gillespie for helpful discussions and encouragement.

This work was funded by the National Institute of Allergy and Infectious Diseases (grant U19AI110820).

[^1]: **Citation** Chung M, Munro JB, Tettelin H, Dunning Hotopp JC. 2018. Using core genome alignments to assign bacterial species. mSystems 3:e00236-18. <https://doi.org/10.1128/mSystems.00236-18>.
