Phylogeography of Human and Animal Coxiella burnetii Strains: Genetic Fingerprinting of Q Fever in Belgium

Q fever is a zoonotic disease caused by the bacteria Coxiella burnetii. Domestic ruminants are the primary source for human infection, and the identification of likely contamination routes from the reservoir animals the critical point to implement control programs. This study shows that Q fever is detected in Belgium in abortion of cattle, goat and sheep at a different degree of apparent prevalence (1.93%, 9.19%, and 5.50%, respectively). In addition, and for the first time, it is detected in abortion of alpaca (Vicugna pacos), raising questions on the role of these animals as reservoirs. To determine the relationship between animal and human strains, Multiple Locus Variable-number Tandem Repeat Analysis (MLVA) (n=146), Single-Nucleotide Polymorphism (SNP) (n=92) and Whole Genome Sequencing (WGS) (n=4) methods were used to characterize samples/strains during 2009-2019. Three MLVA clusters (A, B, C) subdivided in 23 subclusters (A1-A12, B1-B8, C1-C3) and 3 SNP types (SNP1, SNP2, SNP6) were identified. The SNP2 type/MLVA cluster A was the most abundant and dispersed genotype over the entire territory, but it seemed not responsible for human cases, as it was only present in animal samples. The SNP1/MLVA B and SNP6/MLVA C clusters were mostly found in small ruminant and human samples, with the rare possibility of spillovers in cattle. SNP1/MLVA B cluster was present in all Belgian areas, while the SNP6/MLVA C cluster appeared more concentrated in the Western provinces. A broad analysis of European MLVA profiles confirmed the host-species distribution described for Belgian samples. In silico genotyping (WGS) further identified the spacer types and the genomic groups of C. burnetii Belgian strains: cattle and goat SNP2/MLVA A isolates belonged to ST61 and genomic group III, while the goat SNP1/MLVA B strain was classified as ST33 and genomic group II. In conclusion, Q fever is widespread in all Belgian domestic ruminants and in alpaca. We determined that the public health risk in Belgium is likely linked to specific genomic groups (SNP1/MLVA B and SNP6/MLVA C) mostly found in small ruminant strains. Considering the concordance between Belgian and European results, these considerations could be extended to other European countries.

Q fever is a zoonotic disease caused by the bacteria Coxiella burnetii. Domestic ruminants are the primary source for human infection, and the identification of likely contamination routes from the reservoir animals the critical point to implement control programs. This study shows that Q fever is detected in Belgium in abortion of cattle, goat and sheep at a different degree of apparent prevalence (1. 93%, 9.19%, and 5.50%, respectively). In addition, and for the first time, it is detected in abortion of alpaca (Vicugna pacos), raising questions on the role of these animals as reservoirs. To determine the relationship between animal and human strains, Multiple Locus Variable-number Tandem Repeat Analysis (MLVA) (n=146), Single-Nucleotide Polymorphism (SNP) (n=92) and Whole Genome Sequencing (WGS) (n=4) methods were used to characterize samples/strains during 2009-2019. Three MLVA clusters (A, B, C) subdivided in 23 subclusters (A1-A12, B1-B8, C1-C3) and 3 SNP types (SNP1, SNP2, SNP6) were identified. The SNP2 type/ MLVA cluster A was the most abundant and dispersed genotype over the entire territory, but it seemed not responsible for human cases, as it was only present in animal samples. The SNP1/MLVA B and SNP6/MLVA C clusters were mostly found in small ruminant and human samples, with the rare possibility of spillovers in cattle. SNP1/MLVA B cluster was present in all Belgian areas, while the SNP6/MLVA C cluster appeared more concentrated in the Western provinces. A broad analysis of European MLVA profiles confirmed the hostspecies distribution described for Belgian samples. In silico genotyping (WGS) further identified the spacer types and the genomic groups of C. burnetii Belgian strains: cattle and goat SNP2/MLVA A isolates belonged to ST61 and genomic group III, while the goat SNP1/MLVA B strain was classified as ST33 and genomic group II. In conclusion, Q fever is widespread in all Belgian domestic ruminants and in alpaca. We determined that the public health risk in Belgium is likely linked to specific genomic groups (SNP1/MLVA B and SNP6/MLVA C) mostly found in small ruminant strains. Considering the concordance INTRODUCTION Q fever is a zoonosis caused by Coxiella burnetii, an intracellular gram-negative bacterium. Phylogenetic studies based on 16S rRNA classified C. burnetii within the Legionellas order of the g-Proteobacteria (Stein et al., 1993). Although C. burnetii is able to infect a wide range of host species (Cutler et al., 2007), the primary animal reservoirs relevant for human infection are domestic ruminants (cattle, goat, sheep). In infected ruminants, the main symptoms are reproduction disorders along with abortion, stillbirth and infertility (To et al., 1998;Woldehiwet, 2004;Agerholm, 2013). Excretion of the microorganisms occurs massively during abortion, with the release of billions of bacteria in the placenta and birth fluids. The bacteria are also secreted via milk, feces and vaginal discharges after parturition for a variable time depending on the ruminant species (Berri et al., 2001;Arricau Bouvery et al., 2003;Guatteo et al., 2006;Rodolakis et al., 2007;Rousset et al., 2009). An additional source of human infection, especially in remote areas, might be represented by exposure to wild animals through tick bites or direct contact with wildlife (Jado et al., 2012;Gonzaĺez-Barrio et al., 2016a).
Humans are easily infected through inhalation of contaminated aerosols (Maurin and Raoult, 1999). In most cases the pathology presents no or flu-like symptoms that may clear up in a few weeks, while in other cases (1%-5% of chronic manifestations) the infection can lead to serious life-threatening issues. Therefore, it is crucial to track and mitigate transmission routes from primary animal infection sources to humans. Current epidemiological data indicate human outbreaks related largely to goat and sheep infections (Roest et al., 2011;Georgiev et al., 2013).
Given the tight connection between ruminants' infection and the risk of human exposure, it is important to establish the relationship between animal and human strains, to be able to identify the likely risk for public health. Molecular typing allows to identify strain-host interconnections and to determine transmission dynamics, useful for the development of targeted intervention or prevention strategies (Riley and Blanton, 2018;Riley, 2019). Nevertheless, C. burnetii genotyping is challenging due to the difficulties in obtaining animal, but especially human samples suitable for the genetic characterization. Considering the complexity of Q fever direct diagnosis and the low chronicity rate, human samples positive for C. burnetii are extremely limited, except during outbreak episodes. In addition, the inefficient and difficult culturing of C. burnetii strains from clinical or animal samples impedes the direct detection of these pathogens by other methods than PCR.
For the molecular characterization of C. burnetii, a harmonized reference method is yet to be designed. Initially, C. burnetii isolates were classified in six distinct genomic groups (GG I-VI) via DNA fingerprinting (Hendrix et al., 1991). Nowadays, this technique is substituted by more discriminant and reproducible methods, even though the Hendrix GG classification has never been abandoned. In fact, GG are also predicted from recent results obtained with newer methods (i.e. MST, ParSNP, CanSNP, MLVA, microarray) (Beare et al., 2006;Hornstra et al., 2011;Karlsson et al., 2014;Piñero et al., 2015;Hemsley et al., 2019). The GG I includes strains from tick, cattle and humans, while group II and IV are found in humans and small ruminants. GG III, V and VI are the most species-specific as they comprehend isolates from big ruminants (GG III), humans (GG V) and rodents (GG VI) (Hendrix et al., 1991;Hemsley et al., 2019). In 2006, two novel genetic groups, GG VII and GG VIII, were identified (Beare et al., 2006). The GG VII consisted of wild animal strains, and the GG VIII included strains that originated from small ruminant and human samples (Jado et al., 2012;Gonzaĺez-Barrio et al., 2016b).
Following the Dutch outbreak, a single nucleotide polymorphism (SNP) typing scheme was developed with the advantage of being a rapid method, providing however an intermediate discriminatory power. A panel of 10 SNPs was selected: 7 located in single copy genes and 3 in the multicopy gene IS1111 (Huijsmans et al., 2011). Successively, other SNP schemes were developed based on the identification of polymorphisms present in MST loci or in intragenic conserved regions (Hornstra et al., 2011;Karlsson et al., 2014).
In Belgium, the prevalence of C. burnetii remained fairly unknown until 2009, when the Dutch outbreak pushed the Belgian authorities to initiate active surveillance of Q fever in small ruminants and cattle farms. From December 2009 to March 2013, goat farms located in 6 out of 10 Belgian provinces were positive for C. burnetii (Boarbi et al., 2014). MLVA and SNP genotypic types of strains in goat (Mori et al., 2013;Boarbi et al., 2014;Dal Pozzo et al., 2015) and SNP genotypes in cattle samples (Dal Pozzo et al., 2015) provided an initial fragmented viewpoint of Q fever in animals. However, there was a complete lack of genetic data in sheep and other animal species as well as an absolute absence of genotypic information from clinical human cases. To assess the likely animal sources important for human contamination in Belgium and to advance a One Health viewpoint of Q fever in the country, this work investigated genetic characteristics and spatial segregation of C. burnetii strains based on MLVA, SNP and whole genome sequence (WGS) data arising from strains cultured in axenic, cell culture and egg environments. Genotyping results and clinical presentation of human cases are presented and discussed, together with the large set of data obtained in different animal species. We showed that Q fever is widespread in Belgium and we identified C. burnetii as an abortive agent in alpaca, which has never been detected before.

Sample Collection
The prevalence analysis included samples from the abortion official program (total analyzed samples N=40,457 from cattle, N=349 from goats, N=1,310 from sheep) that is based on the mandatory/voluntary report of abortions that are successively analyzed for the presence of C. burnetii. Along this program, collected samples consists of placenta, cotyledon, vaginal swabs, fetal stomach content, fetal spleen, liver or any other relevant samples. The prevalence study incorporated also alpaca samples (N=33) retrieved from routine passive surveillance of farms in reference settings and tested for the presence of C. burnetii.
The genotyping study included Q fever positive samples from the abortion official program (N=103 from cattle, N=2 from goats, N=11 from sheep), bulk tank milks (BTM) from the BTM monitoring program (N=22 from goats, N=4 from sheep), alpaca sample from routine passive surveillance (N=1) and biological human samples (N=5 cardiac valves, aortic aneurism and blood) analyzed in the context of human National Reference Centre activities.
Culturing of C. burnetii C. burnetii strains were isolated from infected animal samples using two different strategies, including (1) in vivo isolation in OF1 mice (Charles Rivers, Wilmington, MA) and amplification in embryonated chicken eggs as previously described (Mori et al., 2017); (2) in vitro propagation in Vero cells followed by axenic cultivation in Acidified Citrate Cysteine Medium-2 (ACCM-2) as previously reported (Omsland et al., 2008;Omsland et al., 2009). Both procedures were adopted as strategy (1) was the only one available when strain isolation started. In addition, not all strains were able to grow in ACCM-2 (Kersh et al., 2016), therefore strategy (1) must be used for certain isolates.

DNA Extraction and Diagnostic Real-Time PCR
DNA was extracted either directly from field samples for MLVA and SNP genotyping, or from C. burnetii culture (embryonated egg or axenic culture) for whole genome sequencing (WGS) studies. DNA was extracted from clinical human and veterinary samples with the MagMax ™ Isolation Kit (Applied Biosystems; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. DNA was extracted directly from 200 µl of BTM, or homogenized abortion material or human samples (max. 1 gr in 1 ml MilliQ water). Next, 1/10 th of the eluted DNA was analyzed by real-time PCR targeting the IS1111 repetitive element for the detection of C. burnetii as previously described (Klee et al., 2006;Mori et al., 2013). Results of the real-time PCR assays, performed with the 7500 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) or Light cycler ® 480 Instrument II (Roche Molecular Systems, Inc., US), were expressed as cycle threshold (Ct) values. Samples displaying a Ct-value < 40 were considered to be positive.
DNA derived from embryonated eggs or axenic culture used for the WGS was obtained with the silica-based column method of the DNeasy Blood & Tissue Kit (Qiagen ® , Hilden, Germany) from 200 µl sample following the manufacturer's instructions. DNA was analyzed for the presence of C. burnetii by real-time PCR as described above. The total DNA quality and concentration was assessed by Nanodrop 1000 (Isogen Life Science, De Meern, Netherlands) measurements.

MLVA and SNP Genotyping
DNA extracted directly from human and animal samples was analyzed for MLVA and SNP typing. We refer to these data as in vitro data. DNA positive samples having a Ct-value ≤ 30 were used for MLVA typing performed for 13 markers (MS3, MS12, MS21, MS22, MS30, MS36, MS23, MS24, MS27, MS28, MS31, MS33, MS34) (Arricau- Bouvery and Rodolakis, 2005;Svraka et al., 2006;Tilburg et al., 2012). Amplicons were analyzed by capillary electrophoresis on a CEQ 8000 Genetic Analysis System (Beckman Coulter, Indianapolis, IN, USA). The fragment sizes were defined by comparison to a 600 bp internal standard (Beckman Coulter, Indianapolis, IN, USA). Only MS12, whose unit size (126 bp) is not compatible with capillary electrophoresis, was run on 1% agarose gel electrophoresis and the PCR product size was determined according to a 100 bp DNA-ladder (Invitrogen; Thermo Fisher Scientific, Inc.). The number of repeats in each marker was established by extrapolating the sizes of the obtained amplicons relative to those obtained from the Nine Mile strain (RSA 493), used as a reference control strain. For MS33, readjustments of marker units were additionally operated following WGS analysis. The SNP typing was performed on C. burnetii positive DNA samples by real time PCR as previously reported (Huijsmans et al., 2011).

Genotype Clustering and Geographical Distribution
MLVA and SNP data were analyzed using BioNumerics version 6.6 software (Applied Maths, Belgium). Minimum spanning Tomaiuolo et al. Phylogeography of Belgian C. burnetii Strains trees were built on categorical data to assess the genetic relationship between C. burnetii profiles. MLVA analyses were performed on Belgian data (this study) and on public data, available from the MLVA online database 1 and published in the literature. The geographical distribution of Belgian profiles was determined using the online LocalFocus Go tool 2 .

WGS Sequencing and In Silico Genotyping
Only DNA extracted from embryonated eggs or axenic culture was used for WGS, since the DNA quantity deriving directly from field samples was insufficient for sequencing. WGS was performed on Illumina MiSeq platform using standard protocols (Garcia-Graells et al., 2020). For sequences obtained from embryonated eggs, read quality and filtering was calculated by ShortRead 1.16.331 package (Morgan et al., 2009) (Bioconductor 3 ) and FastX 0.0.1333 (FastX-toolkit 4 ), while adapters were trimmed with cutadapt 1.2 (Martin, 2011). PhiX and Gallus gallus 4.0 (galGal4) contaminants were removed using bowtie 2.0.0-beta5 (78.65% and 70.05% of reads were removed for the bovine CbBEB1 and the caprine CbBEC1 genomes, respectively). De novo assembly of processed pair reads was executed into CLC Genomics Workbench version 6.0.2. For sequences obtained from axenic cultures, read quality control and trimming were performed with FastQC (Galaxy Version 0.72) and Trimmomatic (Galaxy Version 0.38.0). De novo assembly was realized using SPAdes (Galaxy Version 3.12.0). Contigs were subsequently analyzed for in silico genotyping. We refer to these data as in silico data. In silico MLVA typing was achieved for the 13 markers described above using Primersearch tool (Galaxy Version 5.0.0.1), which indicated the amplicon size corresponding to specific primer pairs of each marker. In silico SNP typing was performed via CLC Sequence Viewer version 8.0 using the probe/primer sequences described in Huijsmans et al. (2011) for the Dutch SNP panel, in Karlsson et al. (2014) for the canonical SNP (CanSNP) panel of intragenic region, in Hornstra et al. (2011) for the classification based on the CanSNP panel of MST loci. For the in silico MST, 10 primer pairs (Cox2, Cox5, Cox18, Cox20, Cox22, Cox 37, Cox51, Cox56, Cox57, Cox61) described in Glazunova et al. (2005) were used to identify the corresponding spacer sequences via CLC Sequence Viewer version 8.0. Spacer types were then assigned by blasting the obtained sequences against all C. burnetii spacers collected in the C. burnetii MST database 5 .

Statistical Analyses
Prevalence data and 95% confidence interval (C.I.) were calculated using EPITOOLS 6 . The Hunter-Gaston discriminatory Index (HGDI) was calculated using Comparing Partitions tool 7 to estimate the discriminatory power of individual VNTR loci (Hunter and Gaston, 1988). To determine differences between and within clusters, the analysis of similarity (ANOSIM) (Clarke, 1993) was carried out using Primer-e Software Version 7.

Human Cases of Q Fever
In this study, we present five PCR-positive human cases associated with chronic (patient 1 to 4) and acute (patient 5) Q fever. Clinical presentations consisted of endocarditis, vascular infections and fever (the characteristics are summarized in Table 1). Three cases (patients 1 to 3) were submitted for diagnosis because of a blood culture negative endocarditis, while in one (patient 4) symptoms were related to a persistent abdominal aortic aneurysm. Endocarditis from patient 2 was discovered during a valve replacement following an aortic valve stenosis. In one case (patient 5), positive PCR result were detected in blood. The patient presented asthenia, irregular and nocturnal episodes of fever and pancytopenia. Positive serology was detected in four out of five patients (patient 1,2,4,5; data not shown). The most likely region of infection of all patients was determined based on medical request information and/or practitioner/patient interviews by local public health services. All cases where the likely contamination occurred in Belgium were associated with a SNP1 genotypic type, and MLVA clusters B5, B6 and B7. This genotype is closely related to the genotype of the strain responsible for the Dutch epidemic. The genotype of patient 4 belonged to SNP6/MLVA cluster C2 differing in repeat units size in MS22, MS36 and MS33 from similar strains found in Belgium.

Animal Cases of Q Fever
Samples from the abortion official program provided data to estimate apparent prevalence of Q fever cases linked with abortion in cattle between 2011 and 2015 (after 2015, Q fever was no more systematically tested in all cattle abortive samples) and in sheep and goats between 2010 and 2019 ( Table 2). The apparent prevalence of C. burnetii in bovine abortion samples was 1.93% (95% CI: 1.80%-2.07%), with an annual incidence ranging from 1.26% to 2.84% (N on average=8091). The highest incidence was observed in 2014. Beyond these years, C. burnetii positive rate in bovine could not be estimated due to the interruption of the program and a drastic decrease in sampling (<100 per year). C. burnetii apparent prevalence in goats was 9.19% (95% CI:6.59%-12.69%), with a range of annual incidence between 0%-17.65% (N on average=35) and a peak in the year 2019. Prevalence in sheep was 5.50% (95% CI: 4.39%-6.87%), reaching similar annual incidences as in goats 0-18.34% and a peak in 2017. Interestingly, among the passive surveillance samples, two cases of abortion positive to C. burnetii were detected in alpaca (Vicugna pacos). Differential diagnosis excluded Toxoplasma gondii, Staphylococcus sp., Candida lambica, Brucella sp., Campylobacter sp., Chlamydia sp., fungi which for all was negative. One of the two alpaca samples could be included in this genotyping study.
In addition, this study included goat and sheep BTM positive cases, gathered during the BTM surveillance program (2009ongoing). In this case, C. burnetii positive rate has been described elsewhere (Boarbi et al., 2014).

Molecular Phylogeny of Belgian C. burnetii Strains
MLVA and SNP analysis were conducted to provide a comprehensive viewpoint of C. burnetii strains responsible for animal and human infections in Belgium. All human samples were characterized at the genetic level because they were extremely rare and difficult to obtain, by contrast only positive animal samples with a Ct value ≤ 30 were selected for the genetic characterization. Above this threshold, the genetic analysis could be unsuccessful due to a low bacterial DNA load. A total of 146 C. burnetii positive samples were included for MLVA typing with 13 markers (MLVA13) ( Table S1). Only samples with a complete MLVA13 profile, or missing maximum two out of 13 markers, were used to build up the minimum spanning tree (n=114) and to evaluate genetic relationships between strains. The Unweighted Pair Group Method with Arithmatic Mean (UPGMA)-generated minimum spanning tree graph ( Figure  1A), including references strains (n=8), depicted a central group including most of the cattle samples, some goat and sheep samples and giving origin to two smaller divisions represented by other goat, sheep samples and all human strains. Interestingly, within the central group no human samples were present. To identify specific phylogenic patterns, the minimum spanning tree built on the MLVA13 data was collapsed to include strains with similarity above 90% (only one marker difference). This resulted in a more compact tree divided in three clusters (A, B, C) and different sub-clusters (A1-A12, B1-B8, C1-C3) validated by the ANOSIM results (R = 0.851, p = 0.1%).
The subcluster A1 included most of the cluster A samples (91%), indicating the presence of relatively homogeneous profiles within cluster A. By contrast, the collapsing did not change sample distribution within clusters B and C, confirming their heterogeneity and division in separate subclusters.
In parallel, a total of 93 samples were characterized by SNP typing, revealing the presence of 3 major genotypes over the different years in Belgium: SNP1, SNP2, and SNP6 ( Figure 2A and Table S2). The Nine Mile reference strain was confirmed as SNP3 type. The SNP2 profile prevailed upon Belgian small and big ruminant samples and defined the alpaca genotype. No human strain belonged to this genomic group, just as in the MLVA cluster A. By contrast, the SNP1 and SNP6 types were mainly detected in goat, sheep and human samples with a single case of abortion in cattle associated with an SNP1 profile, and one milk cattle sample presenting a SNP6 profile. To compare the MLVA and SNP typing, 52 samples characterized by both methods were used to construct two minimum spanning trees ( Figure 2): one using SNP data ( Figure 2A) and the other one using MLVA13 data ( Figure 2B). Figure 2B showed that the SNP divisions were preserved among the MLVA classification, with however increased discriminant power conferred by the

Geographical Distribution of Belgian C. burnetii Strains
The geographical distribution of C. burnetii strains circulating in Belgium was assessed by both MLVA and SNP data. Plotting MLVA clusters and subclusters on the Belgian map linked phylogeny with geographical distribution of each strain with a high degree of discriminatory power. The geographical localization of the MLVA profiles displayed that cluster A was dispersed over the territory, although it was particularly concentrated in the West Flanders (2) and East Flanders (3) provinces (56%) ( Figure 3A). Cluster B was randomly present in all Belgian areas, while cluster C appeared more concentrated in the Western provinces (75%). The geographical distribution of the Belgian C. burnetii strains based on the SNP classification confirms that the SNP6 (cluster C) genomic group is more localized in the Hainaut province (7) (89%), situated in the Western part of the country ( Figure 3B).

Phylogeny of European C. burnetii Strains
Considering that live animals are traded every day in Europe, we examined the relationship between Belgian and other European (EU) C. burnetii strains. A database containing 955 MLVA profiles from 13 European countries was created assembling  data from the present study, the MLVA online database and published articles (Table S5). Due to the lack of uniformity in the choice of markers, the comparison of all European profiles was challenging. For this reason, together with a large MLVA profiling (MLVA12, Figure S1), an MLVA7 analysis was performed to include more countries, and giving a broader view on the EU C. burnetii phylogeny (Figure 4). In the two analyses, complete profiles or those missing maximum 2 marker values, were used to build minimum spanning trees grouped by species ( Figure 4A and Figure S1A) and by countries ( Figure 4B and Figure S1B). Most of the European cattle strains were grouped together in a central section, including Belgian caprine strains. Some French human samples from the 1990s were present as well. Strains from caprine, human and other animals (bull, rodent, lamb) were split into two lateral groups ( Figure 4A). The division in three main clusters described for the Belgian strains could be extended to France, Italy, Spain, Portugal and Hungary as they presented samples in all three clusters. By contrast, Dutch, Polish and Croatian samples were limited to the central and one lateral group, indicating the possibility of two major genomic groups ( Figure 4B). The CbNL01 strain, responsible for the Dutch outbreak, was present in Dutch, Belgian and French samples.

WGS of Belgian C. burnetii Strains
In order to have a deeper and resolute view about C. burnetii genetic characterization, the genome of four Belgian isolates was sequenced ( Table 3). Considering the complexity to sequence DNA directly from field samples, strains were first cultured, using the egg amplification strategy and/or the axenic culture isolation method, and successively sequenced. Isolates were selected with the aim of getting a representative bovine and/or caprine genome for each Belgian genetic group, although strain sequencing strongly depended on the success of culturing, which is extremely difficult for C. burnetii (Kersh et al., 2016). Two bovine (CbBEB1, CbBEB2) and two caprine (CbBEC1, CbBEC2) isolates were successfully cultured, sequenced and further analyzed for the in silico genotyping. All sequences presented an intact lipopolysaccharide (LPS) DNA sequence which characterizes virulent Phase I strains (data not shown), contrary to avirulent phase II bacteria which possess truncated LPS.
In silico SNP typing confirmed in vitro SNP results, with the CbBEB1, CbBEB2 and CbBEC2 strains belonging to the SNP2 type, while the CbBEC1 strain belonged to the SNP1 type. The in silico MLVA profiles of the sequenced strains were congruous with the in vitro results, except for three values (corresponding to MS24 marker for CbBEB1 and MS34 marker for CbBEB2 and CbBEC2), located in the microsatellite panel ( Table 3).
The sequence of these isolates allowed in silico MST-typing and to establish for the first time the spacer types (ST) corresponding to the Belgian C. burnetii strains. Cattle and goat SNP2 isolates belonged to ST61, while the goat SNP1 strain was classified as ST33.
Sequences were also grouped using the CanSNP panel in MST loci which confirmed that the caprine SNP1 strain was close to  the human and goat reference strains, while the SNP2 isolates were classified along with other ruminant reference strains (Table 4). Additionally, by comparison with published phylogeny (Hornstra et al., 2011;Hemsley et al., 2019), the SNP2 isolates were classified as genomic group III, while the SNP1 isolate was assigned to the genomic group II as described in Hendrix et al. (1991). The CanSNP panel in the intragenic region did not allow a separation between SNP2 and SNP1 isolates, therefore it did not allow to discriminate between genomic group III and II (Table S6).

DISCUSSION Q Fever in Human and Animal Cases
Q fever is a worldwide zoonosis caused by infection with C. burnetii. Domestic ruminants are the principal sources of human contamination (de Bruin et al., 2012), and there are suggestions that wild animal might occasionally play a role too (Jado et al., 2012;Gonzaĺez-Barrio et al., 2016a). Consequently, it is important to understand the transmission dynamics within a specific epidemiological context and to tackle the likely routes  (Table S2) and the host species. from reservoir animals critical for the public health. In this study, Q fever was investigated and detected in all domestic ruminants and in humans. In all, but one human case the disease was chronic, with the major clinical manifestation associated with endocarditis, which generally represents 60%-70% of all Q fever chronic cases (Brouqui et al., 1993;Maurin and Raoult, 1999). One case was due to vascular infection, that together with osteomyelitis, hepatitis, encephalitis and other rarer symptoms, represents the minor manifestations of the chronic form of the disease (Hatchette and Marrie, 2001;Manchal et al., 2020). The last case was a febrile patient with a positive PCR result in blood, indicating a probable recent involvement with the disease. The patient presented also pancytopenia, a rare manifestation of acute Q fever.
In animals, C. burnetii was detected in abortion of domestic ruminant species at different degrees of apparent prevalence. In bovine samples, C. burnetii annual incidence remained stable over the years (from 1.26% to 2.84%); by contrast, in goat and sheep abortion products, it was more susceptible to fluctuations (from 0% to~18%). Of note, the sample size of the bovine samples was bigger than the caprine/ovine ones (thousand samples for bovine, hundred samples for caprine/ovine), highlighting the importance of large sample sizes for a more accurate estimation of C. burnetii prevalence. Differences in   In silico genotyping of isolates sequenced after isolation by in vitro propagation in Vero cells followed by axenic cultivation in ACCM-2. b In vitro genotyping of isolates directly from positive animal samples. c In silico genotyping of isolates sequenced after in vivo isolation in OF1 mouse strain and amplification in embryonated chicken eggs. d Coding convention according to Arricau-Bouvery et al. (2006).
NA, Not available due to unsuccessful typing.  Tick  USA  1935  A  T  T  T  C  G  G  T  G  T  G  G  G  T  /  16  I  Henzerling a  Human  blood   Italy  1945  G  G  T  T  C  G  G  T  G  T  G  C  G  T  Br.II.001  18, 25  II   CbBEC1  Goat milk  Belgium  2010  G  G  T  T  C  A  G  T  G  T  G  C/G  G  T  Br.II.007  33  II b  CB2 a  Human  blood   France  /  G  G  T  T  C  A  G  T  G  T  G  C  G  T  Br.II.007 11-15,  24,  32-34   II   CbBEB1 Cow abortion material monitoring programs, sampling strategy and testing methods impedes a clear comparison of prevalence data among different countries and the estimation of Q fever true prevalence in Europe and worldwide (Guatteo et al., 2011;The European Union One Health 2018Zoonoses Report, 2019. In this study, we found C. burnetii in abortion material of alpaca. To the best of our knowledge, this is the first report of association of abortion with a positive PCR result on Q fever in alpaca. Infection of camelids by C. burnetii represents an important issue in several countries (i.e. Egypt, Algeria, Iran, Kenia, Saudi Arabia), where considerable levels of seroprevalence (Benaissa et al., 2017;Browne et al., 2017;Khalafalla et al., 2017;Mobarez et al., 2017;Klemmer et al., 2018) and some abortion cases in camels were described (El-Deeb et al., 2019), but none has previously been observed in alpaca. Investigation on the role of C. burnetii in abortion in these animals and their potential role as reservoir for human contamination deserve certainly attention in the future.

Genotyping and Geographical Distribution of Strains
C. burnetii molecular typing defined the genotypes circulating in Belgium and highlighted the relationship between animal and human strains, assessing the risk for human infection. MLVA and SNP typing methods were used to have a comprehensive viewpoint of Q fever in animals and humans through samples obtained from 2010 to 2019.
Overall, three MLVA clusters (A, B, C) subdivided in 23 subclusters (A1-A12, B1-B8, C1-C3) and three SNP types (SNP1, SNP2, SNP6) were identified. The high discriminative power conferred by the MLVA method was crucial for the identification of subclusters and accomplished an in-depth analysis. Although the SNP2 type/MLVA cluster A was confirmed to be the most abundant genotype in Belgium, it seemed not responsible for the human cases considered in this study, as it was only present in ruminant samples. However, interspecies transmission and maintenance of this genotype seemed to be high as it was found in samples from cattle, goat, sheep and alpaca. The SNP1/MLVA B and SNP6/MLVA C clusters were mostly found in small ruminant and human samples, with the (extremely rare) possibility in cattle. Similar incidental finding was observed also during the Q fever outbreak in the Netherlands (Huijsmans et al., 2011;Roest et al., 2013). The heterogeneous population structure in clusters B and C suggests that these strains could bear a higher capacity to evolve, possibly linked with the improved ability to infect more host species, including humans. Increased genome plasticity and mutations in membrane proteins and predicted virulence-associated genes were indeed identified in the CbNL01 strain, identical to our SNP1/MLVA cluster B strains (Kuley et al., 2017). Based on the overall results, we suggest a genomic group-host specificity, pointing at SNP1 and SNP6 genomic groups as responsible for human infections in Belgium. Both genotypes may cause not only chronic, as demonstrated here, but also acute Q fever (Huijsmans et al., 2011). In Dal Pozzo et al. (2015 a novel genotype SNP5 was reported in Belgium, which was not confirmed in the present study. Noteworthy, Dal Pozzo et al. (2015) reported only a single sample with this genotype, which differs in a single marker from SNP1 type. Detection in a larger set of samples and/or isolation of the strain is warranted to confirm the circulation of this genotype in Belgium.
Regarding the MLVA data, as others, we observed a lack of data for MS23 and MS33 markers, especially in small ruminant samples by both in vitro and in silico MLVA typing. This suggests that also in Belgian strains an increased variability or mutations occurs in these regions, hampering their detection by the two methods (PCR and Illumina sequencing). The instability of these VNTR markers was previously proposed (Arricau- Bouvery et al., 2006;Svraka et al., 2006) and the unfeasibility to establish the number MS23 and/or MS33 repeats by in vitro or in silico MLVA was confirmed in other studies Ceglie et al., 2015;Di Domenico et al., 2018). The presence of the recognition site for the insertion element IS1111 in front of the repeat units of MS23 and MS33 makes these markers easy targets for IS1111 insertion (Sidi-Boumedine et al., 2015).
To acquire a broader view on the MLVA types of C. burnetii strains present in Europe, we created adatabase containing 955 MLVA profiles from 13 European countries grouping published data or those present in available databases. The minimum spanning trees built on these MLVA profiles revealed the presence of 3 main groups: a central one including most of the cattle samples, and two lateral groups containing small ruminants and human samples, reflecting the distribution described for Belgium and in Jouliéet al. (2017). The Nine Mile (tick) strain was located at the extremity of the trees, indicating an important phylogenetic distance with the European field strains. Considering that this strain is used for the production of the Coxevac ® vaccine employed for ruminant prevention, it raises the question of whether it is the best option to protect ruminant herds against Q fever.
In the European population structure, Belgian MLVA profiles are close to French and Dutch genotypes. Indeed, the SNP1, SNP2, and SNP6 types were also found in France and in the Netherlands (Huijsmans et al., 2011). In the Netherlands, the SNP1 and SNP6 groups were associated with human and goat samples (cattle spillover was present as well), while the SNP2 type only included animal strains (Huijsmans et al., 2011), supporting our hypothesis on the likely sources for the public health risk. As also indicated in Mori et al. (2013), the MLVA profile of SNP1 goat samples from 2010 were identical to the CbNL01 strain responsible for the Dutch outbreak. Three of these samples came from Antwerp, a region bordering the Netherlands, and one from the Namur province, situated at the Franco-Belgian border. The origin of the Namur sample is uncertain, although the owner could probably have exchanged animals with the bordering countries or within the country. Internal animal trade could be also responsible for strain transfer inside the Belgian territory. Of note, there has been an evolution of this genotype over time. In 2012, one sample from the Namur area presented an identical CbNL01 MLVA profile except for a Tomaiuolo et al. Phylogeography of Belgian C. burnetii Strains single marker (MS36). Over the years, positive samples differed from the Dutch strain in more than one MLVA marker, keeping yet the SNP1 type features, and spreading all over the territory. It is possible that, as during the Dutch outbreak , the expansion of the C. burnetii SNP1 genotype in Belgium is occurring as well. The SNP2 group was also distributed throughout the country, even if Dal Pozzo et al.
(2015) detected SNP2 bovine samples (2011)(2012)(2013)(2014) mostly in the southern provinces of Belgium. This was probably due to the localization of the farms visited during that study. Instead, the SNP6 type was concentrated in the Hainaut and Namur provinces, two regions bordering France, confirming previous studies (Mori et al., 2013;Dal Pozzo et al., 2015).

In Silico Genotyping
This study presented additionally WGS of four C. burnetii isolates originated from Belgium. Three genomes (the CbBEB1, CbBEB2, CbBEC2) were obtained from SNP2 animal samples (two cattle and one goat, respectively), while an additional one (the CbBEC1) corresponded to a goat SNP1 sample. The in silico genotyping of these isolates was used to complete our genomic study. For the first time, the spacer types corresponding to C. burnetii Belgian strains were identified by in silico MST. Cattle and goat MLVA A/SNP2 isolates belonged to ST61. This is a new genotype only recently detected in Polish dairy cattle herds (Szymańska-Czerwińska et al., 2019). The goat MLVA B/SNP1 strain was ST33. This ST was also found in human and small ruminant samples (rare in cattle) from France, German and the Netherlands, and it was recently described in goat isolates from the UK (Chochlakis et al., 2018;Hemsley et al., 2019). Samples were additionally classified using two in silico CanSNP typing methods. The CanSNP in intragenic region was based on the identification of SNPs present in conserved region of the core genome (Karlsson et al., 2014). It was created as an alternative to the CanSNP in intergenic MST regions (Hornstra et al., 2011), considered as highly variable regions. In this study, the CanSNP in intragenic region method was unable to discriminate between SNP1 and SNP2 samples. By contrast, we could classify our samples using the CanSNP in MST regions. This classification confirmed that the caprine SNP1 strain was close to human and goat reference strains, while the SNP2 isolates were classified along with other ruminant reference samples. In addition, it was possible to predict the genomic groups (Hendrix et al., 1991) by comparison with the published phylogeny (Hornstra et al., 2011;Hemsley et al., 2019). The SNP1 strain (or ST33) belonged to the genomic group IIb, which included human, goat and sheep isolates (Hemsley et al., 2019). The SNP2 strain (or ST61) was part of the genomic group III, mainly composed by big ruminant samples (Hemsley et al., 2019). European isolates were predominant in the genomic group IIb and III (Hemsley et al., 2019). The genomic group III included also ST20, a genotype not found in Belgium, but present in several European countries (UK, France, Germany, Netherland, Italy, Spain, Hungary) (Glazunova et al., 2005;Astobiza et al., 2012;Reichel et al., 2012;Tilburg et al., 2012;Di Domenico et al., 2014;Sulyok et al., 2014). The ST20 was principally detected in cattle, rarely in small ruminants, and involved in farm animal outbreak (Reichel et al., 2012). During the Dutch outbreak, this strain type was only associated with animal samples and not with human infection Roest et al., 2013). Nevertheless, a few ST20 human cases were detected in France in the 1990s (Glazunova et al., 2005). Interestingly, the minimum spanning tree that we created assembling MLVA European data, presented a central group including most of the cattle strains and a few French human samples from 1990s, suggesting that the association of humans with this genotype could be ascribed to a particular situation during that time frame.
The in silico genotyping analysis assessed the reliability of the in vitro MLVA and SNP results confirming that the SNP1 type is related to human and small ruminant infections, contrary to the SNP2 type that is mostly linked to animal infections.
In conclusion, the public health risk in Belgium appears to be linked to specific genomic groups and not to strains found predominantly in cattle. These conclusions are supported by studies in other European countries (Jado et al., 2012), although a limited number of human samples were applied here. The fact that human cases are more genetically linked to small ruminants strains could be attributed to the different manifestation of Q fever in the ruminant species (abortion storms-and therefore the release of enormous amount of infected material in the environment-in goats and sheep, sporadic cases in cattle) or to the intrinsic genetics of the strains. In either case, interspecies transmission of caprine/ovine genotypes can incidentally occur in cattle (i.e. SNP1 found in a rare case of bovine abortion). In these situations, exposure to cattle might be a risk for human health, therefore regular genotyping surveillance in the bovine sector should be promoted. These considerations could be extended to other European countries, especially those bordering Belgium, because of the concordance between Belgian and other European data as discussed above.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, SAMN16591806, SAMN16591807, SAMN16591808, SAMN16591832, SAMN16591833.

ETHICS STATEMENT
The approbation of the Bioethics Committee was not required in this study because data accessed here had not been collected for research purposes but as part of the routine data collection for epidemiological surveillance in accordance with article 18 of the law of 08/12/1992 of the Belgian Government regarding the protection of the privacy of the individual when dealing with personal data.

FUNDING
The NRC activity is supported by the Belgian Ministry of Social Affairs through a fund of the Health Insurance System. This work was partly supported by the project RF 10/6228 of the Federal Public Service of Health, Food Chain Safety and Environment. The Federal Agency for the Safety of the Food Chain (FASFC) funded the C. burnetii detection during the official surveillance programs in animals.