Phylogenetic inference of Coxiella burnetii by 16S rRNA gene sequencing

Coxiella burnetii is a human pathogen that causes the serious zoonotic disease Q fever. It is ubiquitous in the environment and due to its wide host range, long-range dispersal potential and classification as a bioterrorism agent, this microorganism is considered an HHS Select Agent. In the event of an outbreak or intentional release, laboratory strain typing methods can contribute to epidemiological investigations, law enforcement investigation and the public health response by providing critical information about the relatedness between C. burnetii isolates collected from different sources. Laboratory cultivation of C. burnetii is both time-consuming and challenging. Availability of strain collections is often limited and while several strain typing methods have been described over the years, a true gold-standard method is still elusive. Building upon epidemiological knowledge from limited, historical strain collections and typing data is essential to more accurately infer C. burnetii phylogeny. Harmonization of auspicious high-resolution laboratory typing techniques is critical to support epidemiological and law enforcement investigation. The single nucleotide polymorphism (SNP) -based genotyping approach offers simplicity, rapidity and robustness. Herein, we demonstrate SNPs identified within 16S rRNA gene sequences can differentiate C. burnetii strains. Using this method, 55 isolates were assigned to six groups based on six polymorphisms. These 16S rRNA SNP-based genotyping results were largely congruent with those obtained by analyzing restriction-endonuclease (RE)-digested DNA separated by SDS-PAGE and by the high-resolution approach based on SNPs within multispacer sequence typing (MST) loci. The SNPs identified within the 16S rRNA gene can be used as targets for the development of additional SNP-based genotyping assays for C. burnetii.


Introduction
Coxiella burnetii is the etiological agent of Q fever, a zoonotic disease which can result in large outbreaks because of its low infectious dose. The main route of transmission to humans is inhalation of contaminated aerosols which are capable of traveling on wind currents for miles [1]. Although less common, outbreaks have been reported following ingestion of unpasteurized milk and dairy products as well as through contact with contaminated clothing [2,3]. Clinical manifestations of the acute form can range from flu-like symptoms to pneumonia and in some cases result in death caused by respiratory distress [4]. As the acute form presents with non-specific symptoms, healthcare providers are often challenged with accurate diagnosis for Q fever. Centers for Disease Control and Prevention (CDC) recommends a combination of serology testing complemented with polymerase chain reaction (PCR) of clinical specimens for definitive laboratory diagnosis of early stage acute Q fever [5]. Persistent focalized C. burnetii infections, the chronic form, can present as endocarditis [4,[6][7][8], vascular or osteoarticular infections [9,10] and are primarily diagnosed by serological testing. The low infectious dose (less than 10 bacteria), high resistance to environmental stress in its nonreplicating small-cell variant (SCV) form, and potential for deliberate dissemination are the basis for classifying C. burnetii as a select agent regulated by the Federal Select Agent Program (https://www. selectagents.gov/SelectAgentsandToxinsList.html) [11]. This highly virulent pathogen has been weaponized by both the U.S. and former U.S.S.R. under various biological warfare programs [12] and is considered a bioterrorism concern for both civilians and military personnel.
The epidemiology of Q fever is complex due to its world-wide geographic distribution and diverse reservoirs that include wild and domestic animals, livestock, arthropods, and humans [4]. Epidemics of this zoonosis have been reported in residential areas as a result of dissemination from nearby farms [13], among troops during wartime [14][15][16], and can be hyperendemic in certain regions with sporadic cases spanning several decades [17]. The largest reported global outbreak occurred in the Netherlands between 2007 and 2010 [18]. Despite the high prevalence of C. burnetii in the environment, epidemiological data are limited because isolate collections are usually small and are governed by select agent regulations [19]. The advent of host cell-free anexic media and the contribution of genomics were major advances in the study of Q fever. Whole genome sequencing (WGS) of C. burnetii has revealed genetic diversity between strains and has enabled the development of new molecular typing methods that will help bridge gaps in epidemiological knowledge and define phylogenic relatedness. It has been demonstrated that the epidemiological features and clinical presentation of acute Q fever can be strain dependent [20]. Simple and rapid laboratory genotyping methods would aid in traceback investigations during natural or intentional outbreaks or other public health events.
Numerous molecular typing methods for C. burnetii have been described in the literature and applied to epidemiological investigations involving naturally occurring outbreaks and a suspected intentional release [21]. Initial typing methods used before 2005 were based on plasmid types [22], pulsed-field gel electrophoresis (PFGE) [23], restriction fragment length polymorphisms (RFLP) analyzed with SDS-PAGE [24], and sequence analysis of genes including com1, mucZ, icd and 16S/23S [25][26][27][28]. These methods were shown to have varying degrees of discriminatory power and their routine use could not be established as a result of limited reproducibility within and between laboratories [21]. Recently, more high resolution sequence-based genotyping methods such as multispacer sequence typing (MST) [29] and multiple locus variable number of tandem repeats (VNTR) analysis (MLVA) [30] have become the most commonly used approaches. Although the new genotyping methods have increased discriminatory power between C. burnetii strains, converting to these new methods frequently reduces the amount of information originally obtained by more traditional methods, as data may be discarded or cannot be directly compared [19]. According to Hornstra et al. [19], new typing methods should be unambiguous, simply transferrable and allow for C. burnetii typing results to be compared to existing collections. In 2011, Hornstra and colleagues used the publicly available MST scheme from the largest and most diverse C. burnetii collection in the world [29] to design SNP assays that could be used to type isolates based on previously described phylogenetic groups.
Phylogenetic analysis, mainly based on 16S ribosomal RNA (rRNA) gene sequences, has placed Coxiella in the Class Gammaproteobacteria along with Legionella, Francisella, and Rickettsiella [27,31]. It has also assigned all Coxiella isolates to the single C. burnetii species [27,32]. While 16S rRNA gene sequencing has confirmed the phylogenic homogeneity of C. burnetii strains, reports suggest this method does not have sufficient discriminatory power to be used for genotyping [27]. It has also been demonstrated that other classical genotyping methods based on 16S-23S internal transcribed spacer sequencing and RNA polymerase's β-subunit (rpoB) sequencing are not useful for epidemiological studies due to high levels of sequence homology [28,33]. Here, we demonstrate 16S rRNA sequencing does indeed differentiate C. burnetii strains and can be used as a simple, unambiguous, transferrable genotyping approach. To maintain synchronization between typing methods, strains used in this study that overlapped with strains genotyped by other methods were directly compared.

Materials and methods
Bacterial strains and DNA extraction C. burnetii strains included in this study and respective epidemiological information are listed in Table 1. To synchronize results between typing methods we compared 16 strains that overlapped with the Hendrix et al. [24] study and 39 strains that overlapped with the Hornstra et al. [19] study. Genomic DNA was isolated using the QIAamp DNA Mini Kit (Qiagen, Valencia, CA) following the manufacturer's "DNA Purification from Tissues" protocol with overnight proteinase K lysis at 56˚C.
16S rRNA gene amplification and Sanger sequencing 16S rRNA genes from 55 C. burnetii strains were amplified by PCR using the universal, bacteria-specific primers E8F and E1541R to obtain an amplicon of approximately 1500 bp [34]. PCR amplification reactions were performed as outlined in Hakovirta et al. [35] and products were confirmed by gel electrophoresis using a 0.8% E-Gel1 electrophoresis system (Life Technologies, Eugene, OR). ExoSAP-IT (ThermoFisher Scientific, Pittsburgh, PA) was used to hydrolyze any remaining PCR primers in the amplification reactions and DNA cycle sequencing reactions were performed with the Big Dye1 Terminator v 3.1 cycle sequencing kit (Ther-moFisher Scientific, Pittsburgh, PA). 16S rRNA gene sequences were generated using five previously described, bi-directional oligonucleotide primers [34,35] and two additional primers, E341R (TGCIGCCICCCGTAGG) and E1115F (CAACGAGCGCAACCCT), designed to increase sequencing coverage at the 3' and 5' ends. Reaction products generated from cycle sequencing were purified using the BigDye1 XTerminator Purification kit (Life Technologies, Eugene, OR) and DNA was sequenced using an Applied Biosystems 3500xL Genetic Analyzer and KB Basecaller Software v1.4.1 (Life Technologies, Eugene, OR).

DNA sequence analysis and identification of 16S rRNA SNPs
16S rRNA contig assemblies, alignments, and identification of SNPs were performed with Sequencher™ v5.0 software (Gene Codes Corporation, Ann Arbor, MI). Prior to assembling a   16S rRNA consensus sequence, the electropherogram for each DNA sequence was visually examined and trimmed for quality. Fifty-five 16S rRNA consensus sequences generated in this study were compared with other 16S rRNA gene sequences in Ribosomal Database Project (RDP, http://rdp.cme.msu.edu) and in GenBank of NCBI (https://blast.ncbi.nlm.nih.gov) using the Basic local alignment search tool (BLAST) [36] to confirm the C. burnetii species. Gene sequences for all strains generated herein, along with two additional sequences obtained from NCBI for C. burnetii Z3055 (Accession number NZ_LK937696) and 3262 (Accession number CP013667), were then aligned to the 16S rRNA reference sequence of C. burnetii Nine Mile. Position numbers for the five SNPs identified in this study were assigned based on alignment to this reference sequence. A discriminatory power calculator (http://insilico.ehu.es/ mini_tools/discriminatory_power/index.php) was used to calculate discriminatory power (D). D is a single numerical index of discrimination expressed by the formula of Simpson's index of diversity and is based on the probability that two unrelated strains sampled from the test population will be placed in different typing groups [37]. These values were calculated based on 36 strains and 6 types for the Hendrix et al. [24] study, 50 strains and 10 types for the Hornstra et al. [19] study and 50 strains and 6 types for this study. As epidemiological data is limited for many of the historical C. burnetii isolates, relatedness could not be established between all strains. In our study, the six known epidemiologically related strains, GP-MT2, ES-WA1, GP-WA1, ES-MT1, ES-WA2 and GS-MT1, were considered as one for the purpose of calculating D.

SNP-based genotyping assay
Hornstra et al. [19] identified polymorphisms within MST loci and designed a genotyping assay based on 14 SNPs. SNP data for 39 of the 55 C. burnetii strains examined in this work have previously been published [19]. SNP assay results for the remaining 16 strains were obtained as part of this work. Briefly, 12 SNPs were used to develop Melt-Mismatch Amplification Mutation Assays (MAMA), and melt curves were analyzed according to Vogler et al. [38]. Two SNPs were used to develop TaqMan minor groove binding dual-probe assays, and results were analyzed as described by Easterday et al. [39]. For all 14 assays, 1 μL of DNA was used in a total reaction volume of 10 μL. 1 x SYBR Green PCR Master Mix and 1 x TaqMan Genotyping Master Mix (both by Life Technologies, Foster City, CA) were used for the Melt-MAMA assays and TaqMan minor groove-binding dual-probe assays, respectively. Thermal cycling conditions, allele-specific primers, instruments and software details are all outlined in Hornstra et al. [19].
Phylogenetic tree assembly 16S rRNA gene sequences were aligned using the multiple sequence alignment method MUS-CLE and overhanging ends were trimmed in the software suite MEGA v6.06. Pairwise distances between all C. burnetii strains were computed for both the 16S rRNA and genomic SNP sites, and neighbor-joining was used to infer phylogeny from the observed pairwise distances. Stability of the 16S rRNA tree was tested by bootstrapping 500 times.

Results
Identification of 16S rRNA gene single nucleotide polymorphisms C. burnetii strains included in this study and the SNPs within their respective 16S rRNA genes are listed in Table 1. A diverse set of 55 C. burnetii strains isolated from various geographical locations and host reservoirs were chosen for analysis in this study. 16S rRNA gene sequences from these strains were determined for amplified product with lengths greater than or equal to 1500 nucleotides. When sequences were aligned and analyzed, six polymorphisms were identified among the 55 strains; five SNPs and one nucleotide insertion. Nucleotide position numbers for these six polymorphisms were assigned based on sequence alignment to the C. burnetii Nine Mile reference strain and their locations are illustrated in green in Fig 1. These six positions are all located within or directly flanking one of nine determined variable regions based on a variability map produced for the Escherichia coli 16S rRNA gene [34,40]. The nine regions within the 16S rRNA gene of the C. burnetii Nine Mile reference strain corresponding to the variable regions reported in E.coli are highlighted in yellow (Fig 1). Three of the five SNPs identified at nucleotide positions 164, 622 and 1022 are pyrimidine transitions with either a thymine or cytosine present in all strains. The SNP located at position 64 represents a purine transition where a guanine was confirmed in 39 of the 55 strains and an adenine in the remaining 16. A transversion mutation involving an adenine/cytosine exchange represents the SNP at position 820. C. burnetii Mauriet was the only strain in our study in which an insertion was identified at nucleotide position 213 (Fig 2). This indel was detected on DNA sequences obtained using both the forward and reverse sequence that cover this region. Visual inspection of the two electropherograms generated from these DNA sequences also confirmed this insertion (Fig 2A). The 16S rRNA sequence of the Mauriet strain was then compared to those found in the Gen-Bank database of NCBI using BLAST. Two C. burnetii strains, Z3055 (accession number NZ_LK937696) and 3262 (accession number CP013667), were shown to possess an A inserted at position 213 which is absent in the Nine Mile reference strain sequence (Fig 2B). Although not detected using the BLAST search, this indel was also identified in the 16S rRNA gene sequence of C. burnetii Cb109 (data not shown).

Genotyping by 16S rRNA SNP analysis
The neighbor-joining clustering method was used to create a phylogenetic tree from pairwise distances observed from SNP sites between all C. burnetii strains (Fig 3). The five SNPs and one indel identified among the 55 16S rRNA gene sequences gave rise to six genotypes. The largest number of isolates, including the Nine Mile reference strain, was assigned to Group 1 and possesses the SNP signature GCCAC at positions 64, 164, 622, 830, and 1022, respectively. Strains clustered in this group have worldwide distribution originating in North America, Africa, Central America, Europe, and Australia ( Table 1). The 16S rRNA gene of C. burnetii Mauriet contains the same five SNPs as Group 1, and the indel at position 213 situates this strain separately in Group 6. The five SNPs in the 16S rRNA sequences of strains Z3055 and 3262 obtained from GenBank were determined to be identical to those confirmed in the Mauriet strain and would also place them in Group 6 (data not shown).
One SNP separates both Groups 2 and 5 from Group 1 where GCTAC is observed for Group 2 and GCCAT for Group 5. The geographic distributions of strains assigned to these two groups are more concentrated compared to Group 1. Group 2 includes 7 strains isolated in the mid-20 th century from Europe and Russia whereas strains in Group 5 originate from Nova Scotia and the northern United States, all of which were isolated in the 1980s ( Table 1). Strains assigned to Group 3 form a separate cluster at the initial bifurcation point of the phylogenetic tree and contain the SNP signature ACCCC (Fig 3). Members of this Group are associated with Australia [41] and the HHV-WA1 strain was isolated from a person with travel history to that country and who was likely exposed there. Thirteen C. burnetii strains occupy Group 4 on the next closest branch and differ by only a T transition at position 164. The majority of strains associated with Group 4 were isolated from goats or environments from the north western part of the United States. The three human heart valve isolates in this group, K, P, and HHV-WA2, were also found in the same region of the United States located in Oregon, California, and Washington states, respectively.

Typing method comparisons
To evaluate the discriminatory power of the 16S rRNA SNP-based genotyping method, we compared our inference of C. burnetii phylogeny to inferences obtained using other molecular typing methods. Sixteen isolates overlapping with those used in the Hendrix et al. [24] study, which defines genomic groups based on restriction enzyme banding patterns, were compared. In addition, 39 overlapping isolates included in the Hornstra et al. [19] study that used SNP signatures among MST loci were also assessed. These comparative results are outlined in Table 2 along with observed MST genotypes reported by Hornstra et al. [19]. Based on DNA restriction fingerprints, Hendrix et al. [24] differentiated C. burnetii isolates into six distinct genomic groups which is comparable to the results obtained herein. The strains assigned to 16S rRNA Groups 3 and 6 were not represented in the RE/SDS-PAGE study and vice versa. All 16 strains in common from the remaining four groups fell within their corresponding genomic groups ( Table 2). The discriminatory powers for these two typing schemes are similar with D values of 0.7151 and 0.7524 for the 16S rRNA gene-based method and the RE/ SDS-PAGE method, respectively.  The genotyping assay based on 14 SNP sites described by Hornstra et al. [19] provides discrimination between the genomic groups reported by Hendrix et al. [24] with some resolution within these groups. While SNP signatures for 39 of the 55 study isolates have previously been reported, signatures for the remaining isolates were examined by Melt-MAMA or dual-probe assays as part of this work to increase the dataset for comparison (S1 Table). The predicted genomic groups and observed MST genotypes for these additional 16 strains can be found in Table 2. To illustrate phylogenetic inference of C. burnetii based on these SNP signatures, the neighbor-joining clustering method was also used to create a tree from pairwise distances observed from 14 SNP sites within MST loci (Fig 4). While 14 polymorphic sites led to increased resolution within this tree, all but two of the 55 strains were assigned to either the correct predicted genomic group or located on the nearest branch indicating the closest evolutionary relationship. This can be illustrated by the color-coded 16S rRNA Groups indicated in Fig 4. All the isolates that clustered into the 16S rRNA Groups 2 and 5 fell within their corresponding predicted genomic group. Twenty-two of the 24 isolates of 16S rRNA Group 1 also fell within the predicted genomic group 1 with observed MST genotypes of 16 and 26. While the two remaining 16S rRNA Group 1 isolates, CM-CA1 and CM-SC1, belong to a different MST genotype within predicted genomic group 3, they are the nearest neighbors on the tree.
n/a n/a -n/a -3262 (NCBI) n/a n/a -n/a -Comparative genotyping results between the 16S rRNA gene-based method and two other typing methods. RE/SDS-PAGE genomic groups were described by Hendrix et al. [24]. Predicted genomic groups and observed MST genotypes were described by Hornstra et al. [19]. Similarly, while the 16S rRNA Group 4 isolate C. burnetii GP-AF1 populated an individual branch due to its unique MST genotype, it is also located closest in proximity to the other 16S rRNA Group 4 isolates indicated in green (Fig 4). Despite having eight fewer polymorphisms sites for typing, our 16S rRNA gene-based method did exhibit some resolution within the predicted genomic group 4 and could accurately place strains into their more discriminatory MST genotypes. For example, C. burnetii HHV-WA1, Resson and Arandale with observed MST genotypes of 1-7, 30 were assigned to 16S rRNA Group 3, whereas MST genotype 8 isolates ES-CA1, HHV-WA2, K, Pricilla and P were assigned to 16S rRNA Group 4. The difference in discriminatory power between these two typing schemes is less than 0.07, with D values of 0.7151 and 0.7837 calculated for the 16S rRNA gene-based method and the MST loci-based method, respectively. The only two strains that were not consistent by group were C. burnetii Mauriet, which was assigned to the predicted genomic group 2, and Dugway 7D77-80, which formed a single branch close to the initial bifurcation point. While the genomic group prediction for Mauriet was 2, this strain is distinctive in that it was the only one in the Hornstra et al. [19] study with observed MST genotypes 11-15, 24, 32-34. The MST loci-based SNP assays performed herein for C. burnetii Dugway 7D77-78 placed this strain in the predicted genomic group 6, which could not be discriminated by 16S rRNA gene typing. This strain, as well as Dugway 5J108-111, was isolated from rodents in Utah [19] and was not assigned a numerical MST. It is suggested that these strains could comprise a novel genotype which is unlike Dugway 5G61-63, a tick isolate that falls within 16S rRNA Group 1 and was previously shown to belong to Genomic Group 1 [42] in the Hendrix et al. study [24].

Discussion
With a broad range animal reservoir and with its potential misuse as a bioterror weapon, an unambiguous laboratory method that can genotype C. burnetii and help trace infections back to their source during an outbreak is crucial. Our analysis of the complete 16S rRNA led to the classification of all study strains into six distinct groups and allowed us to infer C. burnetii phylogeny in concordance with other typing methods. 16S rRNA genes contain nine variable regions (V1-V9) that exhibit substantial sequence diversity that is species-specific and can be useful targets for diagnostic assays. However, it has been demonstrated that no single variable region of the 16S rRNA gene can differentiate among all bacteria [43] and analyzing the entire 16S rRNA gene with the nine variable regions is important for scientific investigations. Variations in intra-genomic copies of the 16S rRNA gene must be taken into account for other bacterial biothreat pathogens [35,44], but C. burnetii contains a single gene copy so there is no potential for sequence heterogeneity. Two previous sequencing studies targeted either the 16S rRNA gene [27] or the internal transcribed 16S-23S rDNA spacer (ITS) region [28] to genotype C. burnetii. Both studies revealed very high levels of sequence homology, but the methods were not useful for epidemiological, clinical or taxonomic purposes [27,28]. The limited number of strains included in these studies likely contributed to the restricted usefulness of the typing method. Through analysis of a larger and more diverse strain collection, we demonstrated that 16S rRNA gene sequencing can be used for C. burnetii genotyping.
To avoid any potential loss of genotyping information between methods, we compared the 16S rRNA gene-based typing results to those defined by restriction enzyme banding patterns in the more historic Hendrix et al. study [24]. All newly typed isolates overlapping with those isolates used in Hendrix et al. [24] were accurately placed into these previously established genomic groups and the discriminatory power of these two typing methods was comparable. These genomic groups also correlate to previous groupings based on plasmid types [22]. Isolates in Groups 1, 2 and 3 all contain the QpH1 plasmid and QpRS is found in the genomes of Group 4 isolates. Members of Group 5 are plasmidless and Dugway isolates in Group 6 contain the QpDG plasmid [24].
We evaluated the congruity of our results and the discriminatory power of the 16S rRNA SNP-based genotyping method by comparing our inference of C. burnetii phylogeny to phylogeny predicted by the Hornstra et al. [19] MST SNP-based approach. Moreover, to expand the dataset for comparison and to build on our existing understanding of phylogeny, we performed MST genotyping for 16 C. burnetii isolates that have not been previously described. This is of particular importance in order to facilitate inter-laboratory comparisons as C. burnetii collections are often small due to the culture and containment requirements. While the calculated D values show that 14 SNPs have more discriminatory power compared to the six sites identified among the 16S rRNA gene sequences, a small number of SNPs is suitable for rooting phylogenies if there is a low substitution rate and resolution surrounding the root of the tree [45]. Indeed, 53 of the 55 newly typed strains were assigned to either the correct genomic group as predicted by Hornstra et al. [19] or located on the nearest branch indicating the closest evolutionary relationship. We observed instances of both increased and decreased resolution when comparing the genotyping schemes. Using the 16S rRNA SNP-based approach, we could discern different MST genotypes within the predicted genomic group 4, but could not resolve predicted genomic group 1 from its closest neighbor on the tree, genomic group 3. Along with typeability, reproducibility, and ease of interpretation, discriminatory power is an important criteria in evaluating bacterial genotyping methods [46]; however, higher resolution does not necessarily correspond to a more accurate inference of epidemiological relatedness and may be less suitable during analysis of many isolates over an extended period of time [47]. Thus, the genotyping method that is chosen for a particular investigation depends on the epidemiological question being addressed and whether the situation pertains to a localized outbreak or to a more large-scale study [47].
One incongruent result between the two typing approaches was observed for the French isolate C. burnetii Mauriet. While the genomic group prediction for Mauriet was 2, this strain is unique in that it was the only one in that group with observed MST genotypes 11-15, 24, 32-34 [19]. The 16S rRNA gene-based SNP method places the Mauriet isolate in a separate but neighboring group from all the other predicted genomic group 2 strains which were assigned different MST genotypes. As Mauriet was the only strain in our collection that contained an insertion polymorphism, we compared its 16S rRNA sequence with those available in Gen-Bank. C. burnetii Z3055, a strain clonal to the one responsible for the large Q fever outbreak in the Netherlands in 2007, also contained this same insertion. C. burnetii Z3055 was isolated in 1992 from an ewe placenta in Germany and has the same MST 33 genotype and VNTR profile as the strain responsible for the outbreak in the Netherlands [48]. It is possible that the French isolate Mauriet is also related to the Dutch epidemic strain as the outbreak most likely spread from Germany to the Netherlands via France [48].
While no laboratory genotyping method is considered the gold standard for C. burnetii [21], the 16S rRNA SNP-based method shows promise for a rapid and unambiguous solution that can be used to type strains into previously established genomic groups. Assays targeting these SNPs only require a real-time PCR system and can be completed with a short turnaround time; approximately 2 hours post-DNA extraction [49]. In addition, sequence-based methods are less likely to suffer the complications of inter-and intra-laboratory reproducibility [21]. The proliferation of affordable next generation sequencing platforms has contributed to a marked increase in microbial sequence data including metagenomic studies focused on full length 16S rRNA genes [50]. Unlike SNPs located in the genome-wide MST loci, the 16S rRNA SNPs identified in this study are localized to a 1500 nt region of the C. burnetii genome. 16S rRNA sequence data collected during metagenomic analyses of clinical or environmental samples could be useful for both identification and now genotyping of C. burnetii strains. In addition, 16S rRNA SNP genotyping assays could be directly applicable to these samples, in which fragmented DNA may be present, because only small regions require amplification [49,51]. Also, due to the very low mutation rate of 10 −9 to 10 −11 per base pair per generation reported for intracellular bacteria [21], the 16S rRNA SNP targets are considered stable. Future work could comprise development of timely and inexpensive assays for detection of the six 16S rRNA gene-based SNPs including; TaqMan minor-groove binding dual-probe and melt-MAMA assays or by pyrosequencing [19,52]. SNP-based genotyping could be employed as the gold-standard method for C. burnetii, but until the field comes to an agreement, harmonization across promising new techniques is crucial.

Disclaimer
The findings and conclusions in this report are those of the authors and do not necessarily represent the views of the Centers for Disease Control and Prevention or the U.S. Department of Homeland Security. The use of trade names and commercial sources is for identification purposes only and does not imply endorsement by the U.S. Public Health Services, the Department of Health and Human Services, or the U.S. Department of Homeland Security. HPM, BC, JRH, RAP, AC, AC, DH, PPS, LMW, GJK and DS declare no competing financial interests.
Supporting information S1 Table. SNP signatures among MST loci for C. burnetii strains. Signatures for these strains, not previously included in the Hornstra et al. [19] study, were examined by Melt-MAMA or dual-probe assays. (DOCX)