In Silico Genome-Genome Hybridization Values Accurately and Precisely Predict Empirical DNA-DNA Hybridization Values for Classifying Prokaryotes

For nearly 50 years microbiologists have been determining prokaryotic genome relatedness by means of nucleic acid reassociation kinetics. These methods, however, are technically challenging, difficult to reproduce, and - given the time and resources it takes to generate a single data-point - not cost effective. In the post genomic era, with the cost of sequencing whole prokaryotic genomes no longer a limiting factor, we believed that computationally predicting the output value from a traditional DNA-DNA hybridization experiment using pair-wise comparisons of whole genome sequences to be of value. While other computational whole-genome classification methods exist, they predict values on widely different scales than DNA-DNA hybridization, introducing yet another metric into the polyphasic approach of defining microbial species. Our goal was to develop an in silico BLAST based pipeline that would predict with a high level of certainty the value of the wet lab-based DNA-DNA hybridization values. Here we report on one such method that produces estimates that are both accurate and precise with respect to the DNA-DNA hybridization values they are designed to emulate.

relating to microorganisms lead to the development of the polyphasic approach [4], whereby two microorganisms are combined into the same species if they share a certain number of phenotypic and genotypic traits. A principal criterion that the query and reference strains should meet to be considered the same species relates to the degree of their gross genome similarity. DNA-DNA hybridization (DDH) is a "wet-bench" molecular biological technique that allows for the direct experimental comparison of two genomes (for a modern high throughput version of this method see Loveland-Curtze et al. [5]). This method follows simple, and seemingly indisputable, logic: the closer two microorganisms are in gene content and sequence similarity, the higher the degree their measured DNA-DNA hybridization value must be. For example, if two microorganisms share half of their genes through homology, this sets the maximum level by which they can possibly be related at the genome level to 0.5. However, homologous sequences can, and usually do, vary from each other which acts to further reduce overall relatedness (i.e. sequence dissimilarity will decrease genome relatedness from 0.5). Thus, the DNA-DNA hybridization value is a direct measure of shared gene content weighted by shared sequence similarity.
The DDH method of classification was first used by McCarthy and Bolton in 1963 [6] to establish the genome level relatedness of bacteria and now serves as the "gold standard" for novel prokaryotic species classification [4,7]. The necessary empirical evidence to justify classifying one prokaryotic species as different from another is 70% genome similarity, or less, as measured by the DDH method [4,8]. Though this method is theoretically elegant, in practice it suffers from certain limitations. These methods are laborious and require a substantial expertise, which limits their application to a few specialized labs. In the post-genomics era, one naturally wonders if there is a straightforward easy-to-implement and inexpensive in silico alternative that could complement, if not replace, its wet-bench counterpart using whole genome sequences.
Indeed, a number of in silico genome similarity measures have been advanced in the last decade for the purpose of classifying microorganisms. An early attempt was the average nucleotide identity method (ANI) developed by the Tiedji group [9] later followed by the genome-genome distance (GGD) methods summarized by Auch et al. [10]. While these methods are capable of making reasonably accurate predictions of species assignments compared to experimental DDH values from the same pair of organisms, neither method produces a metric on the same scale as DDH and therefore requires comparison against some known set of DDH values which relies heavily on prior knowledge. In addition, both of these methods take a "distance-only" approach to determining genome similarity, likely missing out on the information provided by knowledge of shared gene content. Further, due to the great variety of GGD methods to choose from, it remains unclear which method to use under what circumstances, and why.
By reasoning, the information provided by the DDH method should be easily obtained computationally by creating a composite value that combines shared gene content as well as sequence similarity. Incidentally, this should advance a more polyphasic computational approach since shared gene content not only has been shown to recapitulate phylogeny [11], independent of sequence similarity, but it also serves as an indicator of shared physiology and ecology [12,13]. Here we propose a method of directly calculating DDH from whole genome sequences, which we term in silico genomegenome hybridization (isGGH). This method accurately and precisely predicts experimental DDH values and serves as a useful computational method that reproduces published experimental DDH values using modern genomics tools. Because it predicts the experimental results of DDH it is directly interpretable by microbiologists. isGGH is also capable of providing convincing classification evidence at different levels of classification hierarchy from incomplete genomes by allowing the analysis of the distribution of similarity scores produced, essentially reducing the method to distance-only approaches taken by others [9,10].
We validated this method by calculating pair-wise isGGH values for a set of finished bacterial genomes previously analyzed by Goris et al. [14] using ANI and again by Auch et al. [10] using GGD methods. We then compared isGGH values to the published ANI, GGD, and empirical DDH values. We also tested the robustness of the isGGH method to undersampling by sub-sampling portions of complete genomes and then comparing the results to finished reference genomes. All species delineations herein are based on genomic information alone and we make no attempt to include true phenotypic similarity and/or ecological uniqueness. Therefore species concepts relating strictly to the phylotype and ecotype are not considered.

Materials and methods:
ANI values with corresponding DDH values were reported in Goris et al. [14], and were subsequently reused by Auch et al. [10] to test the GGD methods. To facilitate a three-way comparison of these methods, the same organismal genome comparisons and DDH values were used again in this study of isGGH. A total of 75 finished bacterial genomes were downloaded from the NCBI ftp website, and one from The Sanger Institute ftp website (see Supporting Table). An all-BLAST-all search was conducted on coding sequences (CDS) using the BLASTN algorithm [15], v. 2.2.25+, with default settings. Initial filtering of the BLASTN results removed all self-hits from the data set as well as all comparisons where the length of the pairwise alignment was less than 70% of the average length of the query and subject sequences. A final filter was implemented which removed all instances of multiple hits of a query sequence to a single subject genome, and vice versa. This final filter was necessary to reduce the great number of hits resulting from transposases, which have been shown to be the most ubiquitous of all genes, particularly in prokaryotes [16,17], greatly inflating the number of BLAST hits found for any pair of genomes in our data set.
The relatedness of an arbitrary pair of genomes is first determined by finding the number of common genes relative to the the number of CDSs of each the query genome and the reference genome, respectively: This provides two values in the same way DDH provides two values, depending on which is used as the test (query) and which the reference (subject). Common CDSs between any pair of species were determined from the filtered BLAST dataset. Similarity for each pair of CDSs from the BLAST dataset was determined using the H-value of Fukiya [18] (see also Ou [19]; here denoted as "h-value" so as not to confuse with Shannon's entropy H): This h-value is the percent identity (p id ) of a pair of sequences along the entire length of the alignment (l aln ) normalized by the length of the query sequence (l query ). The percent identity, alignment length, and query length were obtained from the BLAST dataset. Since the underlying distribution of similarity scores between species is unknown and apparently asymmetrical we used the median h-value to weight the fraction of common genes shared by the two organisms to produce the isGGH values.

Results and Discussion:
Here, we attempted to show that the isGGH method is, at the least, comparable to other computational methods at classifying microorganisms as either the same or different species. In order to compare ANI and GGD on the DDH/isGGH scale each value was transformed linearly based on a regression equation obtained from comparison with published DDH values [10,14]. Table1 shows that the regression equations for ANI and GGD are variable from genera to genera. Because of this clade specific variability, the regression equation used was generated from the total of all points, as recommended by Auch et al [20]. In Figure 1 we show the relationship between the computed isGGH, ANI, and GGD values, after necessary linear transformations using the regression equations, and the corresponding empirically derived DDH value (ANI and DDH values taken from Goris et al. [14], GGD from Auch et al. [10]); all three computational methods explain a large portion of the variance in DDH values (R 2 as a measure of precision), while only isGGH has a slope near unity (a measure of accuracy) -note that without the linear transformation the slopes of ANI and GGD are 3.12 and -126.99, respectively. isGGH maintains similar statistics when the number of comparisons is increased to 71 across a more diverse set of organisms ( Figure 2) [21]. When predicting species assignments corresponding to a 70% DDH cutoff, ANI and isGGH incorrectly classified six of 51 pairs of organisms and GGD misclassified eight of 51 pairs. Of the six errors made by isGGH, four were attributed to E. We are pleased to note that the distribution of h-values (distance only measures) produced has diagnostic capabilities which are very robust to incomplete/undersampled genomes. In Figure 3 we show the shapes of the resulting distributions arising from haphazardly chosen query genomes, from various levels of classification, against a single reference genome (Escherichia coli K-12 MG1655).
We tested how robust histogram shape is to undersampling/missing information by sub-sampling various portions of the query genome without replacement, and found that random removal of up to 95% of CDSs from the query genome had an insignificant effect on the shape of the histogram (p = 0.09) when compared to a complete reference genome (Figure 4). Pearson's chi-square was used to determine the consistency of shape between all histograms, following the normalization steps outlined in Porter [23].

Conclusions:
In precisely estimate empirical DDH values shows that it is a viable alternative to ANI and GGD, which do not mirror DDH. As a result, isGGH has the potential to be integrated seamlessly in to the field of prokaryotic classification and wherever DDH values are used. The isGGH method (blue-squares) is compared against ANI (red-diamonds) and GGD (yellow-triangles) for a set of 51 pairs of genomes obtained from previous studies [10,14]. To facilitate comparison on the DDH scale, which isGGH computes naturally, ANI and GGD datapoints (from the NCBI BLASTN trNF 6 F3 method) were transformed using the slope and intercept from Table 1   Robustness of the isGGH value to genome completeness was tested for 7 genome pairs.
Pairs of genomes were chosen so that their values were spread from low-relatedness to highrelatedness. For each pair, one genome was chosen as the query and the other as the subject. Portions of the query genome were randomly subsampled and compared to the complete subject genome and the isGGH value calculated.