We used MUMmer (Kurtz et al. 2004) to search for homologous genomic regions across five assembled, one partially assembled, and four short-read C. burnetii genomes (Table 1). Homologous regions were aligned and searched for SNPs using SolSNP (http://sourceforge.net/projects/solsnp/ last accessed June 19, 2013). We ensured site orthology by eliminating potential paralogs and regions that mapped to more than one genomic location on an assembled reference genome, as well as requiring alignment of other genomes to 100 bp flanking each side of the SNP. We minimized potential alignment errors by eliminating SNPs immediately adjacent to other polymorphisms. By requiring that all taxa share each site, we did not have to make any assumptions regarding likely character states before a site was lost. Such filtering left us with 11,386 SNP sites. SNPs were in pseudogenes (n = 1924), intergenic regions (n = 1288), ribosomal genes (n = 2), tRNA genes (n = 8) were non-synonymous (n = 5107), and synonymous (n = 3057). We used parsimony criteria and a heuristic search with default options using PAUP (Swofford 2003) to recover a single most parsimonious tree (Fig. 2). We performed 1000 bootstrap iterations to test the robustness of the tree, however as bootstrapping is a poor method for measuring accuracy for trees with little homoplasy (Felsenstein 1985), we also report homoplasies as a more appropriate and direct measure of accuracy (Archie 1996).

As trees are traditionally rooted using outgroup taxa, we initially attempted to root the C. burnetii tree with other species (Table 1). Using the SNP discovery pipeline described above, we included potential outgroup genomes. In order to increase the likelihood of identifying appropriate sites, we dropped the requirement that sites are shared among all C. burnetii genomes because some have deleted regions and may have regions that lacked sequencing coverage. Legionella pneumophila and R. grylli are reportedly the nearest neighbors to C. burnetii (Roux et al. 1997), however as comparisons to C. burnetii yielded few appropriate sites for rooting, we also included Pseudomonas syringae because of high levels of sequence similarity at the 16S locus identified when using C. burnetti as a query to search GenBank (see the Supplementary Material S2).

We found duplicated regions by searching for sequence similarities within the RSA493 genome (megablast of the RSA493 genome against itself). We discarded regions within 1 kb of insertion (IS) elements and transposases to minimize phylogenetic error as these regions are not likely to be evolutionarily stable. One of the resulting 12 duplicated regions longer than 900 bp was discarded as it was not a discrete duplication because it overlapped with itself and was part of another duplicated region. Sequences similar to the remaining 11 duplicated regions were identified in each of the five assembled and one partially assembled C. burnetii genomes using the BLASTN algorithm to determine their presence/absence and identify orthologs of each paralog. Four of the 11 duplicated regions were then eliminated from further analyses. The first eliminated region contained polymorphisms that were all parsimony uninformative and thus uninformative for rooting. In the second, multiple degenerative paralogs were identified, causing concerns of correctly identifying orthologs of each paralog in other genomes. The third region was part of another paralog, and a fourth region was discarded because it was not duplicated in other genomes and had a complex VNTR-like region, also causing concern for correctly identifying orthologs. Of the seven remaining duplicated regions, an additional paralog was identified for two of the regions, resulting in five duplicated and two triplicated regions (Fig. 3). Orthologs were identified by the best match and confirmed by comparing the flanking sequence around the region across the surveyed genomes. Alignments of orthologs and paralogs using ClustalW in BioEdit (http://www.mbio.ncsu.edu/bioedit/bioedit.html last accessed June 19, 2013) were used to find SNP polymorphisms among orthologs for each genome and identify the paralogous SNP site. Polymorphic paralogs for each SNP site were combined and the monomorphic paralogs assigned to a duplicated set of taxa that was used as an outgroup to root the resulting tree (Fig. 4). In our dataset, at least one of the polymorphic character states matched the state of a paralogous SNP. This was not a requirement, but avoids employing assumptions about the likelihood of different nucleotide changes. We used parsimony criteria as described above to recover a single most parsimonious tree. We also performed a Bayesian analysis using MrBayes (Ronquist and Huelsenbeck 2003); as only variable sites are included in the data matrix (Fig. 3d), we converted the nucleotide data into standard characters and employed the Mkv model using equal rates and state frequencies. The resulting tree showed an identical topology to the maximum parsimony tree with clade credibility values of 100% at all nodes.

Open reading frames (ORFs) were predicted in each genome using glimmer3 (Delcher et al. 2007). ORFs were de-replicated by clustering with USEARCH (Edgar 2010) at an ID of 0.9; ORFs were then translated with Transeq (Rice et al. 2000). Each peptide was aligned to each genome with TBLASTN (Altschul et al. 1997) and the bit score value was tabulated. The conservation of peptides in each genome was determined with a blast score ratio (BSR) analysis (Rasko et al. 2005; Sahl et al. 2013); the BSR value can range from 1.0 (exact peptide match) to 0.0 (no significant alignment). Lost genes were identified as those with a BSR value of ≤0.4, equivalent to ≤40% peptide identity over 100% of the peptide length. The functional annotation of selected peptides was performed by the Kegg automatic annotation server (http://www.genome.jp/tools/kaas/ last accessed June 19, 2013). For those peptides that could not be assigned to a Kegg orthology group, an alignment to the Genbank non-redundant (nr) database was performed; the annotation was then transferred from annotated genes to peptides identified in this study. Some genes (n = 32) were contaminants from host cells used for culturing (rabbit and chicken) or another bacterium known to be sequenced at the same time as M44 and thus excluded from further investigations. Genes deleted from the first two phylogenetic branches (Fig. 5a) were identified by mapping possible deletions as close to the root as possible.

Using 11,386 SNPs (Supplementary Material S1), of which 9096 were parsimony informative, we recovered a single maximum parsimony tree on which only 91 characters are inferred to be homoplastic. The consistency index (excluding parsimony uninformative characters) was therefore 0.9894 (Fig. 2). Relationships shown are congruent with other phylogenetic analyses of C. burnetii (Glazunova et al. 2005; Beare et al. 2006; Beare et al. 2009; Hornstra et al. 2011). The high consistency index suggests a completely clonal mode of genetic inheritance, with little, if any, evidence for lateral gene transfer among C. burnetii genomes. The high consistency index also suggests that the substitution rate is very low. Thus at this level of evolution, once a point mutation occurs, there is little likelihood that subsequent mutations will occur at that site and confound the phylogenetic signal (Keim et al. 2004; Pearson et al. 2004). Finally, the high consistency index (low homoplasy) among a large number of parsimony-informative characters distributed across the genome, as well as 100% bootstrap support for each node, suggests a highly accurate phylogenetic topology.

We identified three taxa that share portions of their genomes and could therefore potentially be used to root the C. burnetii phylogeny. For P. syringae, L. pneumophila, and R. grylli, comparisons to the 10 C. burnetii genomes yielded only a single parsimony-informative SNP out of the 9096 parsimony-informative SNPs found among the 10 C. burnetii genomes. Six other parsimony-informative SNPs were found when we dropped the requirement for sites to be present in all C. burnetii genomes. These seven parsimony-informative SNPs result in a tree with a consistency index of 1.0 and place the root with the RSA493 and AustQD genomes (Supplementary Material S2). The lack of homoplasy is suggestive of a well-supported tree, however the reliance on only seven SNPs provides little confidence in the resulting relationships (bootstrap support values range from 59 to 71). Indeed, given the relatively vast evolutionary distance to these outgroups, multiple substitutions per site are likely and, coupled with short and unequal branch lengths, would result in incorrect placement of the outgroup.

After filtering duplicated regions to exclude mobile elements, complex duplications, short duplications, and duplications without parsimony-informative SNPs (see “Finding and Assessing Duplicated Regions” in the Methods section), orthologs of the remaining five duplicated and two triplicated regions were identified (Fig. 3a–c). Orthologous regions were easily distinguished and not confused across paralogs, as orthologs are more similar to each other than they are to paralogs (Fig. 3b). Alignments of orthologs yielded 59 parsimony-informative SNP sites (Fig. 3d). For each SNP site, the polymorphic paralogs were merged and combined with a concatenation of the monomorphic paralogs assigned to a replicate set of taxa (Fig. 3d). In the resulting phylogenetic tree, the point where the two replicate groups intersect (Fig. 4a), indicates the root and either group can be used an outgroup to root the tree (Fig. 4b). The resulting maximum parsimony tree had a consistency index of 0.8551 and 100% bootstrap support for each node (Fig. 2b,c). Of the 49 non-homoplastic SNP sites, 33 placed the root on the branch leading to Q154 and Q321, with 16 on one side of the root and 14 on the other. Of the 10 homoplastic SNPs, 7 placed the root in 5 different phylogenetic positions. The best supported alternative root position was along the branch leading to Q212, albeit with only three SNPs. The topology and root position of trees inferred through Bayesian approaches were identical to the maximum parsimony tree. After determining the root position using paralogous regions, we were able to map this root position onto the more robust and better supported WGS SNP phylogeny (Fig. 4c). This tree and the 11,386 orthologous SNPs can be found in TreeBASE (http://purl.org/phylo/treebase/phylows/study/TB2:S14254).

Genome reduction.—With a well-supported root, the evolution of genetic characteristics in C. burnetii can be better determined. For example, genome reduction is often associated with the emergence of novel pathogens (Moran 2002; Losada et al. 2010) and may provide insights into genes that are superfluous or incompatible with virulence (Bliven and Maurelli 2012). The genome size of C. burnetii, coupled with the presence of pseudogenes, is suggestive of reductive evolution (Seshadri et al. 2003). The Dugway chromosome, however, is 163 kb larger than the smallest chromosome (RSA493), suggesting that genome reduction may not have occurred in all lineages or some genomes may have incorporated DNA from other species. To investigate this, we measured a pan genome size of 2148 genes and found 341 genes present in some genomes, but absent in others (Fig. 5).

Coxiella burnetii replicates in an exclusive intracellular niche (Howe et al. 2002) where there may be few or no opportunities to gain exogenous DNA, yet the differential presence of 116 genes among C. burnetii genomes is more parsimoniously ascribed to addition rather than deletion events. A majority of these genes (104/116) are found only in one genome (98 in Dugway and 6 in Q321). However, BLAST comparisons against the NCBI non-redundant database found only eight matches to other species with a fractional identity above 0.6 and none above 0.8. The lack of homologs in other species suggests that most, if not all of these genes were inherited from the last common C. burnetii ancestor and lost in multiple independent deletion (Fig. 5a) events rather than gained through lateral gene transfer with another species. It is important to note that while 341 genes are differentially present among C. burnetii genomes, many are adjacent to each other and would have been deleted from ancestral genomes in single events.

Functional characterization of the 341 C. burnetii accessory genes can provide insights into pathoadaptation (Bliven and Maurelli 2012), however extensive characterization is lacking with 56% being hypothetical and unknown (Fig. 5b). Many of the genes missing only from Q321 may be due to an assembly error as they are ribosomal genes thought to be essential for protein synthesis (Fig. 5a,b). Grouping the characterized genes into broad functional categories shows that many of these genes are associated with cellular processes and metabolism (Fig. 5b) and may have been required in a more diverse ancestral niche. Similarly, genes coding for enzymes and transcription/translation regulators may also be linked to metabolic functions and may be extraneous or deleterious once metabolic pathways are disrupted. However, some missing genes appear to be important in contemporary pathogen–host relationships, indicating functional redundancies. For example, ankyrin repeat domains may be involved with ensuring a stable replication niche during infections (Voth et al. 2009), LuxR—family regulators thought to control growth within a host (Minnick and Raghavan 2011) and OmpA-like transmembrane domain proteins known to be immunoreactive (Beare et al. 2008). Although the genes missing from different clades vary (Fig. 5a), there is a wide overlap in function (Fig. 5b), however known virulence and information processing genes were only missing from members of a single clade. Also, BSRs (Fig. 5a) and 3057 synonymous SNPs suggest varying degrees of polymorphisms and pseudogenization among genes found in all genomes. These patterns suggest that while broad evolutionary and adaptive patterns may be similar across clades and time, we can also expect lineage-specific fine scale adaptations.

Evolution of virulence.—Different virulence and pathogenicity attributes of C. burnetii may be an example of lineage-specific adaptivity and has been the subject of historical debate (Samuel et al. 1985; Moos and Hackstadt 1987; Stein and Raoult 1993; Thiele and Willems 1994). In humans, host genetic differences and underlying health conditions confound attempts to relate clinical observations to genetic groups. With the exception of the Dugway clade, strains belonging to the clades represented here have all been associated with disease in humans (Glazunova et al. 2005; Hornstra et al. 2011; Huijsmans et al. 2011), and different clinical manifestations have been reported among patients infected with similar strains. Also, due to the large number of potential hosts and lack of reported disease in most naturally infected animals, comparisons of virulence differences have been limited, making animal models particularly valuable. These models measure virulence of acute disease that is similar in humans, but do not model chronic disease. In these animal models, isolates from the clade consisting of Q321, Q154, and Q177, cause infections at low dose challenge in guinea pigs and mice, but show no significant pathology at low or high doses (Moos and Hackstadt 1987; Russell-Lodrigue et al. 2009). Similarly, isolates in the Dugway clade on the other side of the root, were originally reported as low or avirulent for acute disease in rodents and this was recently confirmed in mice and guinea pig models (Stoenner et al. 1959; Stoenner and Lackman 1960; Russell-Lodrigue et al. 2009). Several isolates in this group were obtained from chronic human cases or goat abortions, possibly representing a persistent (chronic) infection in goats as well as humans. Indeed, although variable, overall trends in human acute disease manifestations seem to largely follow predictions based on animal models. Since these two clades (the Q321/Q154/Q177 clade and the Dugway clade) straddle the root, we speculate that the ancestor of C. burnetii was more likely to cause a persistent or chronic, as opposed to acute, infection in humans, although no animal models are available to test virulence required for chronic disease. Isolates from the clade consisting of AustQD, RSA493, M44, and RSA331, however, are highly virulent and cause acute disease in mice and guinea pigs. Isolates from the clade between the Dugway and AustQD/RSA493/M44/RSA331 clades (consisting of Q212 and McMaster) caused an intermediate type of disease (Russell et al. 2009).Given the position of the root, the evolution of high acute virulence may have involved mutations before and after the bifurcation point leading to the Q212 and McMaster clade, leading to a gradual evolution of high acute virulence in the RSA493 clade and an alternate adaptive pathway involving intermediate acute virulence in the Q212 clade.

Geographic dispersal.—Despite knowing the root of the C. burnetii tree, we are still far from establishing the geographic origin of the species or even particular clades. Geographically diverse sample collections are small and rare and most genotyping methods are not amenable to comparisons across collections (but see Hornstra et al. 2011.Geographic assignment of isolates from existing collections suggest that many major clades contain isolates from different continents (Glazunova et al. 2005; Hornstra et al. 2011), hiding patterns of deep paraphyly of isolates from a particular region that could indicate geographic origins.