A summary of the features of all the C. burnetii strains used in this study is indicated in Table 1. The strains from the Netherlands were primarily isolated from aborted placentas of goats, sheep and from heart valves of human chronic infected patients during the Q fever outbreak period (2007–2010; Roest et al., 2011a). This study also involves strains isolated from aborted placentas of goats and cattle from France. Additionally, a few strains were isolated from Q fever infected humans from elsewhere in the world. All the strains were genotyped based on MLVA (Roest et al., 2011a). The majority of the strains used for sequencing were cultured with a minimal number of sub-passages in BGM cells (Roest et al., 2012). A few strains were cultured axenically in acidified citrate cysteine medium (ACCM-2) also maintaining a minimum number of sub-passages to prevent phase variation (Table 1; Omsland et al., 2011). Genomic DNA was isolated using the phenol-chloroform method (Tang et al., 2013) after overnight incubation with ATL lysis buffer and proteinase K (QIAGEN, Hilden, Germany). A prior DNase treatment was performed to the bacterial pellet of cell-culture cultivated strains to eliminate most of the host derived DNA.

Draft genome sequences of all C. burnetii strains (except 602) were de novo reconstructed from Illumina MiSeq paired-end 250 bp reads using SPAdes-3.6.2 (Bankevich et al., 2012). The draft genome sequence of the 602 strain was obtained by using Roche 454XL sequencing technology and was de novo assembled by Newbler 2.6. The Dutch outbreak representative strain NL3262 was fully reconstructed using a combination of sequencing technologies (Illumina MiSeq paired-end 250 bp, PacBio RS and Roche 454XL) as described previously (Kuley et al., 2016). The contigs of the draft genome sequences were scaffolded and joined into an artificial chromosome based on homology with the closest C. burnetii strain from the database using BLAT synteny tool (Table 1; Kent, 2002). Mapping the reads back on the draft genome sequences using Bowtie2 aligner (Langmead and Salzberg, 2012) confirmed its correct synteny and occasionally helped to close minor gaps. All C. burnetii (complete and draft) genome sequences were annotated by the NCBI Prokaryotic Genome Annotation Pipeline (http://www.ncbi.nlm.nih.gov/genome/annotation_prok). The (draft) genome sequences and plasmids of all sequenced C. burnetii strains have been deposited at GenBank under accession numbers indicated in Table 1. Published genomic sequences of several other C. burnetii strains used for phylogenetic studies were obtained from GenBank (Table 2).

All the sequenced strains were genotyped by MLVA using a selection of 10–12 loci described previously (Arricau-Bouvery et al., 2006; Roest et al., 2011a). As an additional verification step, genotypes of the reconstructed genome sequences as well as published genome sequences of strains from NCBI database were checked by in silico MLVA genotype analysis using isPCR (BLAT), MLVA primers and genotype references. In silico genotypes of all C. burnetii sequenced strains corresponded with known/measured MLVA genotypes. The in silico analysis of the database strains gave us information on their MLVA genotype (Table 2). Based on the number of repeats per genome loci a minimum spanning tree method was used to cluster the MLVA genotypes using Bionumerics (version 6.6) with default settings.

Reconstructed as well as database obtained C. burnetii sequences were hierarchical clustered based on their MUMi distance (Deloger et al., 2009) and displayed as a phenogram using the BioNJ algorithm (Gascuel, 1997). The underlying MUMi distance matrix was calculated from the pair-wise non-overlapping maximal unique matches (MUMs; using Nucmer version 3.22; Kurtz et al., 2004). Relative pair-wise distances (MUMi) were obtained by dividing the sum of pair-wise MUMs by the average genome size of the two compared genomes. This MUMi distance varies between 0 and 1 where 0 indicates very similar genomes and 1 is for very distant genomes. MUMi trees were visualized in SplitsTree4 (Huson and Bryant, 2006).

For SNP detection in the sequenced draft genomes and in the complete genome of strain NL3262, reads with low quality (bases with a quality score of ≤ Q20) were removed before SNP calling. High quality reads were used to map against reference NM genome (GenBank: NC_002971.3) using Bowtie2 short-read aligner (Langmead and Salzberg, 2012). SNP calling was performed using bcftools (Li et al., 2009) using default parameters. SNP sequence alignments were made for only those positions in NM reference strain where each of the sequenced strain had sequence coverage. Phylogenetic trees based on SNPs were created by Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) using standard settings (Sievers et al., 2011) and visualized in SplitsTree4 (Huson and Bryant, 2006). SNPs and insertions/deletions in CbNL01 and CbNL12 strains compared to NM were verified at nucleotide level of the mapped reads.

Alignment of the complete genome sequence of strain NL3262 with the sequences of the reference NM, and strains which were related to the outbreak CbNL01 genotype (Z3055, NL-Limburg) (D'Amato et al., 2014; Hammerl et al., 2015) was performed by using Mauve 2.4.0 aligner software with default parameters (Darling et al., 2004). Protein orthologs analysis of C. burnetii amino acid sequences obtained from GenBank files from NCBI were performed using Proteinortho software (Version 5) with default settings (Lechner et al., 2011). Orthologs and unique proteins in sequenced strains were identified from the orthologous list generated by Proteinortho. Additionally, ortholog analysis of sequenced strains was conducted using the proteome of NM strain as a reference. Pseudo-genes were excluded from the ortholog analysis. Orthologs from different strains were plotted by means of Edwards-Venn diagrams using Venerable package (Version 2.0). For a more detailed genomic sequence comparison, Artemis Comparison Tool (ACT) (Carver et al., 2005) was used to visualize whole genome MUMmer alignments. Detailed analyses were performed on virulence related genes such as variants (mutations/deletions) in LPS encoding region, which were performed by comparing C. burnetii genome sequences to LPS encoding region of NM-I (GenBank: AF387640). Also, variants found in the CbNL01 and CbNL12 strains were subsequently pairwise checked against the reference NM strain using ACT and ClustalW (Thompson et al., 1994). We also performed a COG (Cluster of Orthologous Groups) enrichment-analysis on all non-synonymous genes in CbNL01 and CbNL12 strains and analyzed the proportions of different proteins for each COG category.

In the present study we report the new genome sequences of 11 C. burnetii strains. The reconstruction of the genome sequences of the NL3262 and NLhu3345937 strains was described before (Kuley et al., 2016). Table 1 is a summary of the main features of the sequenced genomes. Except for strain NL3262, all the other strains are reported as draft genome sequences. The assembled contigs of these sequenced strains are joined into artificial chromosomes based on scaffolding against the closest genome sequence present in the database. In C. burnetii genome sequences, the insertion sequence (IS) elements are dispersed all over the chromosome and are not found on the plasmid (Seshadri et al., 2003). These elements are repetitive in nature and as a consequence, the majority of contig-breaks in the draft genome assemblies mapped to the positions of an IS sequence, when analyzed against the NM reference genome sequence (GenBank: NC_002971.3). As expected, all C. burnetii sequenced genomes possess a single circular chromosome of around 2 Mb and a QpH1 plasmid of 37 kb size (except the Schperling and Scurry strains). The bacterial genomes of all sequenced strains are predicted to contain between 1800 and 2062 coding sequences (CDSs). Around 35–39% of these CDSs were annotated as encoding hypothetical products based on the NCBI Prokaryotic Genome Annotation Pipeline. In silico MLVA typing of the genome sequences corresponded with similar profiles of strains as measured by lab-based MLVA measurements (Table 2).

Whole (complete and draft) genome-based phylogenetic analysis was used to infer the relationship between sequenced C. burnetii strains of different genotypes. The strains were initially isolated from different hosts such as goats, sheep, cattle and human patients (both acute and chronically infected; Table 1). The analysis also included whole genome sequences of C. burnetii strains that have become available in the NCBI database (Table 2). Genome sequences of some of the strains from the database, present as draft sequences (several contigs), were scaffolded similarly to our reconstructed genomes and joined into artificial chromosomes based on the reference NM strain. The phylogenetic analysis of these genome sequences were based on MUMi distance. This approach takes into account the number of maximum unique and exact matches (MUMs) of a given minimal length shared by the two genome sequences being compared (Kurtz et al., 2004; Deloger et al., 2009). The phylogenetic analysis was conducted for genomes excluding the plasmids, as presence or absence of plasmids did not result in any difference in hierarchical clustering of C. burnetii strains.

Based on the clustering analysis, the strains could be divided into 4 main clades (indicated as 1a, 1b, 2 and 3 in Figure 1). The majority of the C. burnetii sequenced strains segregated into two clusters [cluster 1 (1a and 1b) and cluster 2, Figure 1]. These two clusters correspond to the two major genotypes (CbNL01 and CbNL12) of strains isolated from The Netherlands, respectively. Strain NL3262 and the database strain NL-Limburg clustered closely together (Cluster 1a, Figure 1). The genome of strain NL3262 is a complete reconstructed genome derived and curated using different sequencing technologies (Kuley et al., 2016). The NL-Limburg genome is a draft genome generated with a single sequencing technology (PacBio) and its assembly constitutes 4 contigs representing the chromosome (Hammerl et al., 2015). The contigs of NL-Limburg cover the complete genome of NL3262 without any gaps (as visualized by ACT; data not shown). The ends of the contigs in NL-Limburg contain overlapping sequences of around 50,000 bp. Due to these additional sequences (which are probably the result of assembly error) the genome of strain NL-Limburg (represented as ^** in Figure 1) is clustered at some distance away from strain NL3262 (MUMi distance 0.025). Upon removal of these repetitive sequences the MUMi distance between strains NL-Limburg (indicated in cluster 1a without ^**) and NL3262 decreases to 0.01, indicating nearly identical genomes with differences based on 8 SNPs. Draft genome sequences of the database strains (Cb109 and EV-Cb_C13) belonging to the CbNL01 genotype (Figure 2 and Table 2) clustered with the CbNL01 strains (602, CbCVIC1, 3345937, 42785537) originating from The Netherlands in cluster 1b (Figure 1). Strains NL3262 and NL-Limburg belong to the CbNL01 genotype strains as well (Figure 2 and Table 2; Hammerl et al., 2015; Kuley et al., 2016). The major differences between the complete and draft genome sequences of the strains belonging to the clusters 1a and 1b are the genomic locations of highly repetitive transposon sequences. The genome sequences of strains NL3262 and NL-Limburg contain 121 and 128 IS elements respectively scattered around the chromosome. These IS elements are missing in draft genome sequences, most probably due to their presence at contig-break regions, as indicated above. Thus, due to these structural differences, the draft genome sequences cluster very closely to each other (Cluster 1b) relative to the complete genome sequences (Cluster 1a). The complete sequence of NL3262 and the nearly complete genome sequence of NL-Limburg (belonging to genotype CbNL01) differ with a MUMi distance of around 0.05–0.06 from the draft genome sequences in cluster 1b which also belong to CbNL01 genotype.

Cluster 2 includes the draft genome sequences of strains Henzerling, Heizberg and the sequences of the CbNL12 genotype group of strains 601, 701CbB1, 2574, and 18430 (Figure 1). Apart from these, cluster 2 contains the database strains belonging to the CbNL12 genotype (draft genome sequences: CbB1, Cb_B18, and EV-Cb_BK10) as well as the reference NM like genotype strains (NM, CbRSA331, Cb175, and Cb185) (Figure 2 and Table 2). The complete genome sequence of strain Z3055 was placed at an equidistance (MUMi distance 0.02) from clusters 1 (CbNL01) and 2 (CbNL12) (Figure 1). An in silico PCR analysis of strain Z3055 showed few MLVA markers similar to CbNL01 genotype but the majority indicative for the NM genotype. This resulted in the presence of Z3055 in-between the nodes of the NM and CbNL01 genotype strains in an MST analysis (Figure 2). Based on the MST and phylogenome analysis, we assessed the genome sequences of strains Z3055 and NM to be very closely related. Cluster 3 included strain Scurry and the database strain CbuG_Q212, both belonging to the same genotype (Figures 1, 2). These strains clustered separately as they contain plasmid-homologous regions integrated in their genomes unlike other strains (Figure 1 and Table 1; Willems et al., 1997). From the phylogenetic analysis, the Schperling strain was found to be the most distant strain as compared to the other de novo sequenced strains used in this study (MUMi distance 0.05; Figure 1). The genome sequences of strains Q321, CbuK_Q154, and Goat Q177 were most similar to the strain Schperling with a MUMi distance of around 0.025–0.03. However, except for strain Q321, these strains belonged to different genotypes (Figure 2). Hence, a phylogenetic clustering of the strains based on an MLVA classification rather than on the host-origin background indicated a high-level of genome sequence similarity of genotype-specific C. burnetii strains.

A phylogenetic tree was also reconstructed using the identified SNPs among the sequenced strains with respect to the reference strain NM using bcftools (Figure 3; Li et al., 2009). The SNP tree was constructed using Clustal Omega (Sievers et al., 2011). At SNP level, strains belonging to similar genotypes clustered together with clade 1 containing the CbNL01 genotype strains; clade 2 containing the strains Henzerling and Heizberg and clade 3 containing the CbNL12 genotype strains. This indicates that the presence of similar SNPs in the genomes of the same genotype strains relative to the reference strain NM. Strains Schperling and Scurry clustered separately, indicating that these strains contain SNPs distinct from the SNPs present in other sequenced strains compared to the reference strain NM. The data corresponds to the presence of distinct genotypes of Schperling and Scurry compared to other sequenced strains (Figure 2). This SNP based clustering of genotype-specific strains further confirms the highly similar genome sequences of strains belonging to the same genotype.

A comparison of the genome sequences of C. burnetii shows that all sequenced strains are highly conserved in terms of size, number of coding sequences and in nucleotide composition (Table 1). The degree of conservation is illustrated by the comparison of the coding capacity of these genomes: 1447 predicted proteins are orthologs in all sequenced 13 C. burnetii strains in this study (with a total CDS count ranging from 1800 to 2062) among which 1409 of the predicted proteins are orthologs with reference NM. Remaining 38 proteins not shared with NM are majorly hypothetical protein products mostly occurring due to modifications (such as point mutations, small deletions/ insertions of <50 bp) in existing NM pseudo-genes (Supplementary Table 1). The 1409 orthologs among all strains constitute the coding core genome. All sequenced C. burnetii strains contain the QpH1 plasmids, except for strain Schperling containing QpRS plasmid and the plasmidless strain Scurry. The QpH1 plasmids from sequenced strains were all very similar to the reference NM QpH1 plasmid (Accession No. NC_002118.1) with differences on an average of 58 point mutations. The draft plasmid sequence of Schperling was similar to the published sequence of plasmid QpRS (Accession No. Y15898.1) with differences in 63 SNP and a gap of 1810 bp (data not shown).

Ortholog analysis of the proteome of the sequenced C. burnetii strains showed the presence of a few unique proteins annotated as hypothetical protein in CbNL01 and CbNL12 strains (data not shown). Further analysis of genes encoding these unique proteins showed that these genes were present at the beginning or end of contigs in assembled draft genome sequences. As contig ends of draft genome sequences are highly repetitive, genes annotated in these regions were not considered as unique genes, because these genes were truncated transposase genes. Taking this factor into account and other manual annotation corrections, we found no unique complete genes encoded in CbNL01, CbNL12, Henzerling, and Heizberg strains. Instead of unique genes for each strain we did find orthologs that were genotype-specific. These orthologs were shared between strains of the same genotype only and were annotated as hypothetical proteins. When we compared these genotype-specific genes at the nucleotide level to the genome of the reference strain NM, mutations or partial ORF deletion in existing NM pseudo-genes were observed for which altered functionality or regulation cannot be directly inferred (Supplementary Table 2).

By comparing the proteome of the 13 sequenced strains (Table 1) 50 and 98 unique protein encoding genes were identified in strains Scurry and Schperling, respectively. Out of these, 3 and 17 genes are exclusively found in Scurry and Schperling, whereas the other genes were orthologs with the closely related strains CbuG_Q212 and CbuK_Q154, respectively (Figure 1). These unique genes in Scurry and Schperling strains were annotated as hypothetical proteins and result from point mutations in existing pseudo-genes of CbuG_Q212 and CbuG_Q154 (data not shown).

A comparison of the coding capacity of the genomes of strains belonging to the same genotype showed 98, 99.4 and 100% orthologs in predicted proteins of CbNL01, CbNL12, and the human strains (same genotype Henzerling and Heizberg strains), respectively (Figure 4). Fifty-one unique genes (non-orthologs) of NL3262 encode for transposase associated proteins (Figure 4A). These genes are not properly annotated in 602, CbCVIC1, NLhu3345937, and 42785537, because they fall within contig-break regions of the draft genome sequences. The high level of similarities of the proteomes between these strains indicates that the gene content is highly conserved within strains of the same genotype. As all the 13 C. burnetii strains included in this study were initially isolated from different hosts, a comparison was made to assess the percentage of proteins shared by strains derived from each host-species. These comparisons resulted in the presence of 94.6, 100, 87.9, and 87% orthologs in predicted proteins of goat, cattle, human-acute, and human-chronic isolated strains respectively (Supplementary Figure 1). The lower percentage of orthologs shared by strains obtained from various host species is due to the presence of different genotypes within each host-species, except for the cattle-derived strains, as both cattle strains were of same CbNL12 genotype. For example; the strains isolated from goat include both CbNL01 (NL3262, 602, CbCVIC1) and CbNL12 (601) genotypes (Supplementary Figure 1B). Further, the strains isolated from acute and chronic human Q fever patients contain more divergent strains (Schperling and Scurry respectively). Because of this, the protein orthologs of acute and chronic human strains are comparatively less (~87%, Supplementary Figures 1C,D). Taken together, the percentage of orthologs in strains obtained from similar host-species is less conserved than the strains obtained from similar genotype. This clearly indicates the presence of genotype-specific functional characteristics of C. burnetii strains.

To assess the genome architecture of strain NL3262 and to study its putative genome rearrangements, the complete genome sequence of strain NL3262 was aligned with the genome sequence of strains NL-Limburg, NM, Z3055, Cbuk_Q154, CbuG_Q212 and Dugway using the Mauve genome aligner (Figure 5). Among the different C. burnetii strains, the CbNL01 genotype strains NL3262 and NL-Limburg encode the highest number of transposase genes (106-112 IS110 family transposase, 3 ISAs1 family transposase and 12-13 remnant transposase genes, respectively). Other completely sequenced C. burnetii genomes generally contain only 30-60 transposase genes (Beare et al., 2009; D'Amato et al., 2014). Compared to the genome sequence of NL3262, chromosomal rearrangements in strains NL-Limburg, NM, Z3055, Cbuk_Q154, CbuG_Q212, and Dugway strains have resulted in 2, 19, 21, 32, 21, and 23 Locally Collinear Blocks (LCBs) of the same gene content and gene order respectively (Figures 5A–F). NM and Z3055 strains share almost all LCBs indicating a similar genome structure (Supplementary Table 3). Cumulatively, the LCBs of NM, Z3055, Cbuk_Q154, CbuG_Q212, and Dugway represent 116 genomic breakpoints relative to the genome of strain NL3262 (Figures 5B–F). In the genome sequence of strain NL3262, 102 (88%) of the breakpoints occurred at intact or a remnant transposase sequences. These transposase mostly belonged to IS110 family and very rarely to IS30 and ISAs1 family (Supplementary Table 3). The presence of transposase genes at almost all genome breakpoints suggests an important role of homologous recombination in order to establish the observed genome rearrangements of the Dutch outbreak strains. Unlike other strains, the gene order in strains NL3262 and NL-Limburg are nearly identical, yielding only 2 LCBs. The absence of rearrangements between NL3262 and NL-Limburg, indicates a similar genome structure of these CbNL01 genotype Dutch outbreak strains isolated from goat and human patient (Figure 5A).

Genome sequencing of a number of Dutch outbreak strains belonging to the CbNL01 MLVA genotype resulted in 132, 135, and 123 contigs for strains NLhu3345937, 42785537, and CbCVIC1, respectively. In contrast, sequencing of the CbNL12 strains 601, 2574, 18430, and 701CbB1 resulted in 43, 40, 66, and 42 contigs, respectively. As described above, the majority of the contig-break points contained repeat sequences and corresponded mostly to transposase genes when mapped against the NM reference genome. Based on this observation, we anticipate that the observed number of contigs roughly corresponds to the number of transposase genes encoded by the strain. Based on this assumption, the CbNL01 and CbNL12 strains contain around 130 and 47 contigs, respectively. The assumed high numbers of transposons in the draft genome sequences of CbNL01 strains also matches with the number of annotated transposase genes in strains NL3262 (121 intact and remnant transposons) and NL-Limburg strains (128 intact and remnant transposons) belonging to the same genotype group. Taken together, the largest number of transposons are encoded by the CbNL01 Dutch outbreak strains as these strains contain two to three fold more transposase genes compared to other strains of C. burnetii sequenced so far (Beare et al., 2009; D'Amato et al., 2014).

To further asses the differences in the genome sequences of the Dutch strains from different genotypic groups, we performed a detailed NM reference-centric variant analysis. For this analysis we used the complete genome sequences of the NL3262 strain (Kuley et al., 2016) and the draft genome sequence of the CbNL01 and CbNL12 strains. We found an average of 2514 SNPs between the CbNL01 strains and NM. Among these, 417 SNPs were present in intergenic regions, 1041 were synonymous SNPs in coding regions and 1056 were non-synonymous SNP which corresponded to 702 mutated genes. Between the CbNL12 strains and NM an average of 2078 SNPs were found, with 340 SNPs in intergenic regions, 889 SNPs in coding regions (synonymous SNP) and 849 non-synonymous SNPs which corresponded to 581 mutated genes. With respect to NM, these mutated genes in CbNL01 and CbNL12 genotype strains majorly encoded hypothetical proteins, membrane proteins, transporter proteins, DNA repair proteins, translation-related proteins and a few virulence-associated proteins (explained in detail below). Compared to the genome of NM, the genomes of the CbNL01 strains showed a deletion of 3,600 bp (containing ankyrin-repeat containing protein; CBU_0072), a deletion of 2,300 bp (containing several hypothetical proteins; CBU_0877-CBU_0880) and a larger deletion of 9,700 bp (containing the peptidoglycan catabolism proteins; CBU_1101-CBU_1112). A deletion of around 2,250 bp (containing several hypothetical proteins; CBU_0016-CBU_0019) was observed in the CbNL01 and CbNL12 genome sequences, relative to the NM genome sequence (Supplementary Figure 2).

With respect to mutations in membrane protein encoding genes, a total of 64 and 48 genes were mutated in the CbNL01 and CbNL12 strains out of the 104 membrane protein encoding genes in the NM strain. Thirty eight of these genes contained the same mutation in both genotype strains. Additional mutations in 26 hypothetical membrane proteins were specific for all the CbNL01 strains. Although, the exact functions of these hypothetical protein products are not known, changes in these cell surface proteins could potentially contribute to an altered antigenic profile of the CbNL01 strains. Large numbers of mutations were also seen in transporter genes, where a total of 41 and 25 genes were mutated in the CbNL01 and CbNL12 strains (from a total of 75 transporters in NM). Among these, 20 genes contained identical mutations in both genotype strains, whereas 21 genes contained mutations specific for the CbNL01 strains. These transporter genes belonged mostly to the ABC (ATP-binding cassette transporters) and MFS (major facilitator superfamily) family transporters. Further COG analysis of the mutated transporter genes showed that these genes could be primarily involved in transportation of amino acid, carbohydrate and ions (COG categories E, G, and P respectively). Moreover, we also identified mutations in several translation-related proteins in the outbreak CbNL01 strains. Around 29 of these mutated genes were also found in the CbNL12 strains, whereas 15 mutated genes were specifically found in the CbNL01 strains (from a total of 121 translation genes in NM; Supplementary Figure 2).

Recognized virulence factors of C. burnetii are genes involved in primary defense mechanisms against oxidative stress. These include several DNA repair genes involved in base excision repair, nucleotide excision repair, mismatch repair, addAB mediated recombinational repair systems, and oxidative stress enzymes (Mertens and Samuel, 2012). Our variant analysis showed high number of mutations in the DNA repair genes; wherein 13 mutated genes were shared in both genotype strains and mutations in 7 genes were specific to the CbNL01 strains (from a total of 26 DNA repair genes in NM). On the contrary, we did not find any mutations in genes encoding for oxidative stress enzymes among the various strains indicating that these genes are highly conserved. The secondary defense mechanisms of C. burnetii includes the manipulation of host cell processes (Mertens and Samuel, 2012). A Dot/Icm type IV secretion system (T4SS) is predicted to be encoded by the genome of C. burnetii homologous to the L. pneumophila Dot/Icm system. This is an important transfer system associated with the delivery of effector proteins from the bacterial cell to the host cytosol mediating PV formation and other cellular events required for bacterial maintenance and survival (Seshadri et al., 2003; van Schaik et al., 2013; Kuley et al., 2015a). Among the T4SS genes, 4 genes were mutated in both genotype strains and an additional 4 genes were mutated specifically in the CbNL01 strains. Furthermore, 12 effector proteins characterized by the presence of eukaryotic like protein domains were also mutated in the CbNL01 and CbNL12 strains relative to NM. These genes are predicted to be virulence-related as the encoded proteins contain eukaryotic-like domains which may mimic host proteins and module host responses required for successful bacterial infection (Supplementary Figure 2). Taken together, compared to the reference strain NM, genotype-specific mutations of genes were present in the Dutch strains. Furthermore, a relatively larger number of non-synonymous mutations were present in critical genes of the CbNL01 strains compared to the CbNL12 strains.

Repeated in vitro passages of C. burnetii strains have been shown to induce antigenic variation and loss of virulence characteristics due to the transition into a truncated LPS structure. This truncation of LPS is associated with chromosomal deletions of O-antigen coding genes, located in a 38 kb region in the genome resulting in antigenic phase variation from phase I to phase II (Hoover et al., 2002; Denison et al., 2007; Kuley et al., 2015b). The genome sequences established in this study, were compared with the LPS encoding region (CBU0676–CBU0706) of NM-phase I (GenBank: AF387640) to visualize deletions in the LPS encoding region. The human strains Henzerling, Heizberg, Schperling, and Scurry contained no deletions in this region (Table 3). However, a deletion of 201 bp (NM-phase I LPS encoding region: 14604-14805) was consistently observed in this region in all sequenced (601, 2574, 701CbB1, 18430) and database strains (CbB1, CbB18, EV-Cb_BK10) of the CbNL12 genotype (Table 3). This deletion corresponded to a portion of CBU_0686 gene (old gene name: JB153-5), a paralog predicted to encode a 2-Oxoacid dehydrogenase (Hoover et al., 2002). As a paralogous gene, deletion of a portion of this gene might not affect the LPS synthesis as the function might be compensated by the other gene in the genome. Similarly, a deletion of 354 bp (NM-phase I LPS encoding region: 20476-20830) was consistently observed in this region in all sequenced (NL3262, 602, CbCVIC1, NLhu3345937, 42785537) and database strains (NL-Limburg, Cb109, EV-Cb_C13) of the CbNL01 genotype. This deletion corresponded to a portion of CBU_0691 gene (old gene name: JB153-10, Table 3) and is predicted to be involved in virenose synthesis which is one of the unusual sugars only present in the C. burnetii LPS-phase I structure (Toman et al., 1998). Due to the deletion, CBU_0691 gene is frame-shifted and might consequently not encode for enzymes involved in virenose synthesis. We hypothesize that this deletion of 355 bp may induce the start of phase shifting in CbNL01 strains. Since the genome sequences of 602 (CbNL01 genotype) was obtained from in vitro grown cell cultures, we looked back into our previous functional analysis data of the 602 strain cultured in cells maintaining a low passage number. Our transcriptome measurements showed that the CBU_0691 gene was actively transcribed in cells (Kuley et al., 2015b). Although the deletion of a portion of this gene results in frame shifting of CBU_0691, it is still possible that gene CBU_0691 is transcribed but that no protein product is formed. In addition, a complete gene deletion (corresponding to CBU_0677 to CBU_0682) and a deletion of part of the genes (corresponding to CBU_0676 and CBU_0683) were observed in the NLhu3345937 strain of the CbNL01 genotype. The NLhu3345937 strain was passaged in cells for only 11 times and already showed deletion of several O-antigen genes, indicating a high rate of phase shifting in this strain under in vitro culture conditions (Table 3). Apart from these deletions, non-synonymous mutations in CBU_0688 and CBU_0698 were present in all CbNL01 strains relative to NM. Additionally, mutations in 12 O-antigen genes were observed in all CbNL01 and CbNL12 strains, compared to NM (Supplementary Figure 2). Taken together, we observed genotype-specific deletions and mutations in the LPS encoding region at a higher rate in CbNL01 strains than in CbNL12 strains.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.