Occurrence and distribution of Giardia species in wild rodents in Germany

Giardiasis is an important gastrointestinal parasitic disease in humans and other mammals caused by the protozoan Giardia duodenalis. This species complex is represented by genetically distinct groups (assemblages A-H) with varying zoonotic potential and host preferences. Wild rodents can harbor potentially zoonotic assemblages A and B, and the rodent-specific assemblage G. Other Giardia spp. found in these animals are Giardia muris and Giardia microti. For the latter, only limited information on genetic typing is available. It has been speculated that wild rodents might represent an important reservoir for parasites causing human giardiasis. The aim of this study was to investigate the occurrence and distribution of Giardia spp. and assemblage types in wild rodents from different study sites in Germany. Screening of 577 wild rodents of the genera Apodemus, Microtus and Myodes, sampled at eleven study sites in Germany, revealed a high overall Giardia prevalence. Giardia species determination at the SSU rDNA gene locus revealed that Apodemus mice, depending on species, were predominantly infected with one of two distinct G. muris sequence types. Giardia microti was the predominant parasite species found in voles of the genera Microtus and Myodes. Only a few animals were positive for potentially zoonotic G. duodenalis. Subtyping at the beta-giardin (bg) and glutamine dehydrogenase (gdh) genes strongly supported the existence of different phylogenetic subgroups of G. microti that are preferentially harbored by distinct host species. The present study highlights the preference of G. muris for Apodemus, and G. microti for Microtus and Myodes hosts and argues for a very low prevalence of zoonotic G. duodenalis assemblages in wild rodents in Germany. It also provides evidence that G. muris and G. microti subdivide into several phylogenetically distinguishable subgroups, each of which appears to be preferentially harbored by species of a particular rodent host genus. Finally, the study expands the database of sequences relevant for sequence typing of G. muris and G. microti isolates which will greatly help future analyses of these parasites’ population structure.


Background
Giardiasis caused by the flagellated protozoan Giardia duodenalis (syn. G. lamblia and G. intestinalis) is one of the most frequent gastrointestinal parasitic diseases in humans worldwide [1,2]. Giardia duodenalis is considered as a species complex of genetically different but morphologically identical organisms with varying zoonotic potential and host preferences. Humans and a wide range of mammal species including livestock, pets and wildlife are susceptible to these parasites [3]. This species complex consists of at least eight genetically distinct groups (referred to as assemblages A to H) based on sequence analyses of a few genomic markers such as glutamate dehydrogenase (gdh), beta-giardin (bg) and small-subunit rRNA (SSU rDNA) genes [4][5][6][7][8][9]. It has been suggested to consider part of these assemblages as distinct Giardia species due to their genetic distance and differing host preferences [3]. Giardia duodenalis assemblages A and B show the broadest host range, which includes humans and, therefore, they are considered potentially zoonotic [3,8,10]. Assemblages C-H are mainly found in certain mammal taxa: C and D in canids, E in ruminants, F in felids, G in rodents, and H in marine mammals [3,11]. The overall impact of zoonotic transmission on the epidemiology of human infections remains unclear [3].
Besides G. duodenalis, to date two more Giardia species, Giardia microti and Giardia muris, are known to infect mammals. Both species are found in small rodents with supposedly different host preferences. Giardia microti is thought to be a parasite mainly of rodents belonging to the family Cricetidae such as voles and muskrats and G. muris mainly of rodents belonging to the family Muridae [3,[12][13][14][15]. Thus, wild rodents can be infected with different Giardia species, including G. microti, G. muris and various assemblages of G. duodenalis (assemblages A, B and G) [3,[15][16][17][18]. To evaluate the epidemiology and the potential zoonotic risk, molecular surveys for Giardia spp. in various wild rodents and lagomorphs, including the North American beaver (Castor canadensis), muskrat (Ondatra zibethicus), brown rat (Rattus norvegicus), house mouse (Mus musculus) and prairie vole (Microtus ochrogaster), have been attempted to determine the Giardia species or assemblage types. However, systematic molecular studies are still comparatively rare, in particular if one considers the broad range of existing rodent host species [8,9,13,15,[17][18][19][20]. In fact, earlier studies suggested a link of waterborne human giardiasis outbreaks to a source in wildlife, in particular beavers, that led to the classification of Giardia as a zoonotic agent [21,22]. However, the distribution of different Giardia species in rodents of various genera and their geographical distribution based on molecular studies is not sufficiently clarified [21].
Rodents can carry a multitude of pathogens, including important zoonotic viruses, bacteria and parasites [23][24][25] and as infections with Giardia spp. are common in wild rodents it has been discussed whether they could play an important role as a reservoir for these potentially zoonotic parasites as well [17,19,24,26].
The goal of this study was to investigate the occurrence and distribution of Giardia spp., including Giardia species allocation, in wild rodents from different study sites in Germany.
Subsequently DNA was extracted and analyzed by previously described real time-PCR (qPCR) [27] and seminested PCR assays (SSU-PCR) [4,28] targeting the SSU rDNA gene locus. The respective number of samples yielding specific PCR products was 471 of 569 (83%, 95% CI: 79-86%) for the qPCR and 371 of 569 (65%, 95% CI: 61-69%) for the SSU-PCR assay. Of note, in all three rodent genera a higher number of infected animals were detected by qPCR compared to IFA analyses, possibly due to higher sensitivity of the former assay as previously suggested [27,29] (Table 1). The qPCR and SSU-PCR tests overall confirmed the higher Giardia prevalence in Microtus and Myodes compared to Apodemus (Table 1 and data not shown).
Since samples of Apodemus spp. were less frequently positive than those of Microtus and Myodes, we estimated parasite abundance to test whether a lower abundance could in part explain this finding. For this purpose we performed a semi-quantitative analysis of the cyst numbers in fecal samples using the IFA data and of parasite DNA abundance in feces using the threshold cycle (Ct) value of the qPCR results. The latter served as a proxy for parasite numbers. DNA samples of Apodemus feces showed significantly higher Ct-values when amplifying Giardia SSU rDNA compared to samples from Myodes (delta Ct~2; P < 0.001) and Microtus (delta Ct~4; P < 0.001) ( Table 2). This implies on average 4-16 times less Giardia DNA content in Apodemus samples compared to those of Myodes and Microtus (Table 2) and, most likely, reflects the lower abundance of parasite cysts in Apodemus animals confirmed by semi-quantitative analysis of IFA data ( Table 2). The qPCR Ct-values for Myodes and Microtus were also  SSU PCR (P ≤ 0.0001). Comparison between groups was done by using Fisher's exact test followed by multiple testing correction (Bonferroni-Holm procedure) and P-values are presented in Additional file 2: Table S1 significantly different, but such a difference was not noted at the resolution of the semi-quantitative IFA test (Table 2). Thus, Giardia burden is likely to differ among rodents of different genera and decreasing when comparing Microtus and Myodes vs Apodemus.
Giardia muris and G. microti predominate in wild rodents while zoonotic Giardia spp. were rarely detected To determine the Giardia species in wild rodents the sequences of 371 (of n = 571 investigated) PCR products of the semi-nested SSU PCR (see above) were determined, analyzed and compared to reference sequences (Table 3). Overall, 358 sequences could be typed and 355 of these sequences could be assigned to Giardia spp. and three sequences to Octomitus intestinalis, a sister lineage of Giardia spp. [30]. The analysis revealed that 97% of Microtus and 96% of Myodes were infected with G. microti. Apodemus spp. were predominantly infected with G. muris (80%) but a sizeable fraction contained G. microti (18%) (see Table 3 for details). There were apparent differences in the G. muris/G. microti proportion between A. agrarius and A. flavicollis, but absolute numbers of samples from these two species were too low to allow meaningful conclusions (Table 3).
Only five samples (of n = 358, 1.4%; Table 3) were found positive for G. duodenalis, and respective sequences could be assigned to assemblage A (in case of one Apodemus, two Myodes samples) and B (in two samples from Myodes). Attempts to further subtype these samples at different genomic loci (tpi, bg and gdh) were not successful (data not shown).
Of the three O. intestinalis-positive samples (two from Microtus, one from Myodes) two were positive in the IFA analysis and all samples were also positive in the Giardia-specific qPCR analysis. Further subtyping at bg locus (see below) from the two IFA positive samples revealed G. microti as a result indicating possible G. microti/O. intestinalis co-infections in these animals.
Sequence analysis of SSU rDNA, bg and gdh genes revealed high variation in G. microti and G. muris Sequence analysis was performed first on 317 SSU rDNA sequences (277 G. microti, 5 G. duodenalis, 32 G. muris and 3 O. intestinalis) for which a complete sequence from both amplicon DNA strands was available (see Additional file 3: Figure S2 for dendrogram of Bayesian phylogenetic analysis of all SSU rDNA samples, including corresponding data of host species, locality, season, habitat and year of collection). Unique sequences were deduced from this data set and analyzed for possible phylogenetic relatedness which indicated high divergence within G. muris and G. microti sequences (Fig. 1). This prompted us to also type samples at the gdh and bg gene loci by established methods [5,7] (Fig. 2).
Within the 32 G. muris SSU rDNA sequences, 10 unique sequence types were identified with four types represented by ≥ 3 identical sequences. Distance analysis revealed that 31 sequences differed by 0-10 single nucleotide polymorphisms (SNPs) to each other over a fragment length of 245-247 base pairs (bp). In contrast, one sample (ssu406, from M. glareolus) differed by [19][20][21][22][23][24][25][26] SNPs to these 31 sequences and this sequence was identical to the G. muris sequence AF113895 that was used for reference. The separation of the 31  sequences deduced here from current G. muris Gen-Bank entries AF113895 and X65063 was also strongly supported by phylogenetic analyses that revealed a separation into two newly identified sequence clusters. Each of these sequence clusters was preferentially detected in different Apodemus hosts as determined by using phylogenetic trait analysis (Fig. 1b). Analysis of bg sequences identified five G. muris sequence types for this typing fragment. These sequences were also separated from bg GenBank sequence EF455599 used as reference further supporting the existence of phylogenetically distinguishable subgroups within the G. muris population (Fig. 2). No gdh sequences could be amplified of G. muris DNA containing samples and, notably, also no database entry exists for this gene. Giardia microti SSU rDNA differed by 0-17 SNPs to each other over a fragment length of 250 bp. Within the G. microti samples, 106 unique sequences were identified, identical sequences, respectively. The phylogenetic analysis is compatible with a genetic substructure of the population as previously suggested [13] although, expectedly, support of nodes was low due to the short sequence fragment used for the analysis (Fig. 1a). However, analysis of 29 G. microti gdh (unique sequence types deduced from a respective set of 49 total G. microti sequences) and 59 G. microti bg (deduced from a total set of 118 G. microti sequences of this locus) sequences strongly supports the existence of different phylogenetic subgroups in G. microti (Fig. 2). Testing for a significant association with host distribution (Microtus vs Myodes) of the phylogenetic subgroups using phylogenetic trait analysis indicated that subgroups are preferentially harbored by one of the two host genera (Fig. 2).

Discussion
Rodents investigated in this study belonged to mice of the genus Apodemus (A. agrarius, A. flavicollis or A. sylvaticus) or to voles of the genus Microtus (M. agrestis or M. arvalis) or Myodes (M. glareolus). The overall occurrence of Giardia spp. in these animals was very high. Apodemus spp. were mainly infected with G. muris, whereas G. microti was predominantly found in Microtus spp. and M. glareolus. Our data argue that G. muris and G. microti subdivide into several phylogenetically distinguishable subgroups, each of which appears to be preferentially harbored by species of a particular rodent host genus. Only a low proportion of samples (1.4%) contained G. duodenalis assemblages A or B. This implies a very low potential risk for transmission of zoonotic Giardia types from wild rodents in Germany. In contrast, several previous studies estimated a higher potential risk for zoonotic transmission from Giardia-infected wild rodents [17,19,24,26,[31][32][33][34]. However, in only a few of these reports molecular analysis was performed to identify the underlying Giardia species. These studies concerned very different geographical regions and covered a variety of different rodent species: for example, G. duodenalis assemblages G and B in mice (Mus spp.) and rats (Rattus spp.) [17,20], G. duodenalis assemblage B and G. microti in muskrats (O. zibethicus) [9,13], and G. microti in prairie voles (M. ochrogaster) [13]. A study from Sweden found G. muris and G. duodenalis assemblage G in M. musculus and R. norvegicus, and G. microti in the one A. flavicollis analyzed, but no evidence for zoonotic G. duodenalis genotypes [15]. In contrast, beavers (Castor spp.) [9,19] and pet chinchillas (Chinchilla lanigera) in particular [35], were found to harbor G. duodenalis assemblage B, indicating a potentially higher risk for zoonotic transmission associated with these host species. In summary, these and the present data highlight the need for local molecular typing analyses to estimate the actual zoonotic risk for Giardia infections emanating from a particular rodent host.
Previous studies from Poland also demonstrated high Giardia spp. prevalences in A. flavicollis, M. glareolus and M. arvalis [36,37] and, in agreement with our results, showed lower parasite cyst loads in samples of Apodemus vs those from Microtus and Myodes. Systematic typing was not performed in these investigations, but four of the Microtus and Myodes samples were exemplarily typed and similar to our study revealed G. microti [38]. In this context it should be noted that fecal excretion of Giardia spp. is often not uniform. Therefore, quantitative data on individual time points should be interpreted carefully. However, it is assumed that such effects are averaged out when analyzing aggregated data as that presented here where the means of different groups were compared. Future studies will be necessary to investigate whether host-parasite sequence type relationships as reported here are a consistent finding also in other rodent populations.
Many studies of Giardia spp. in wild rodents used microscopy-based detection of parasite cysts and subsequently attempted to type Giardia by sequencing of PCR amplicons of one or a few marker genes. Often this approach was reported to be inefficient which suggests shortcomings in typing protocols. For example, a recent report on 12 rodent samples (ten R. norvegicus, two M. musculus) showed that G. duodenalis assemblage G could be typed at all four tested genomic loci (bg, tpi, gdh, SSU rDNA), but G. muris and G. microti-like samples only at the SSU rDNA locus [18]. Primer mismatches that may negatively affect PCR efficacy or low target DNA amounts that may reduce assay sensitivity could be possible reasons. In agreement with these data we also observed lower typing efficiency at bg and gdh loci in comparison to the SSU rDNA locus. We also recognized mismatches of published primers [4,9,13,28] in the G. muris SSU rDNA reference sequences and, therefore, modified primer sets for our semi-nested SSU rDNA PCR assay. Modification resulted in higher positivity rates in particular for G. muris containing samples. In samples of Microtus spp. and Myodes spp., species predominantly infected with G. microti, the effect was much less pronounced. Thus to complement IFA data, we consider the semi-nested PCR at the SSU rDNA locus as the most reliable approach to detect the corresponding Giardia species. This approach will even detect non-Giardia species represented by O. intestinalis. Alternative PCR methods to discriminate Giardia spp. and G. duodenalis assemblages such as the recently developed tests that amplify 5.8S rDNA or internal transcribed spacers (ITS) of ribosomal gene loci [38] may also be suitable but have not been evaluated in this study. PCR and sequencing of bg and gdh loci was more appropriate though for subtyping and description of population structure due to higher sequence heterogeneity. Because only two bg and no gdh sequences were deposited in public databases for G. muris and none of both for G. microti our work significantly enlarges the database which will greatly facilitate future comparative studies.
Our analysis not only confirmed that G. microti is more closely related to the G. duodenalis complex than to G. muris [13] but allowed also to distinguish genetic subgroups and provides evidence that these appear to preferentially infect different hosts. The existence of such subgroups was previously suggested for G. microti based on a very limited number of ribosomal RNA gene sequences [13] and is now supported by other markers through our analysis. This allows the interpretation that possibly various assemblages exist within G. muris and G. microti populations comparable to what has been described for G. duodenalis. Our study further indicates that G. muris also can be divided in sub-types that are preferentially associated with different hosts, in this case illustrated by A. agrarius and A. flavicollis. Most of the G. muris SSU rDNA sequences newly described were highly divergent from the two published sequences used as references. Both of the latter were derived from samples of non-European rodents [39,40]. The concept of non-uniformity of G. muris is further supported by bg sequences, in spite of the low PCR efficiency at this locus. Further studies are warranted to determine the broader context and significance of phylogenetic relationships (e.g. host range, localities) within G. muris and G microti. This will require improved typing tools.
We also identified by semi-nested PCR at the SSU rDNA locus three sequences from O. intestinalis, a sister lineage from Giardia spp. known to infect wild rodents and together with Giardia spp. form the subfamily Giardiinae [30]. The life-cycle of O. intestinalis is not well characterized and, in particular, information about the morphology and antigenicity of cysts is scarce. All three samples (two Microtus, one Myodes) were also positive in the qPCR analysis and two samples in the IFA analysis. The bg sequence typing was also successful in the latter two samples and revealed G. microti. We therefore consider it likely that co-infections of G. microti/O. intestinalis, with a dominance of O. intestinalis, occurred in the respective hosts and that sequence similarity at the SSU rDNA locus of Octomitus spp. to Giardia spp. resulted in detection of the former. This is also corroborated by the high sequence similarity of the primer sequences used in the current study with published reference and primer sequences of the SSU rDNA sequence of O. intestinalis [30]. Co-infections may also occur with G. microti and G. muris which would impact our interpretations. However, we consider it unlikely that co-infections with different Giardia species are common in our sampled rodent population because, if true, mixed sequences should have been observed frequently.

Conclusions
In the present study, G. muris was found to be present mainly in Apodemus, and G. microti mainly in Microtus and Myodes hosts. Furthermore, our findings argue for a very low prevalence of zoonotic G. duodenalis assemblages in wild rodents in Germany. Evidence is provided that G. muris and G. microti subdivide into phylogenetically distinct subgroups that each prefers species of a particular rodent host genus (in the case of G. microti) or family (in the case of G. muris). The study also expands the database of sequences relevant for sequence typing of G. muris and G. microti isolates. In future studies this will greatly help to analyze the population structure of these parasites in more detail.
Fecal samples (1-2 pellets per animal) were collected from 577 individual animals and preserved in 1 ml H 2 O supplemented with a cocktail of antibiotics in order to inhibit bacterial growth. Samples were shipped at room temperature and stored at 4°C until further processing (storage time between 1 and 10 months). Rodent species were identified mainly by field inspection, in some cases molecular typing was performed on fecal samples using cytochrome b as a target gene as described earlier [42].

IFA test
All samples were examined using the MERIFLUOR® Cryptosporidium/Giardia test (Meridian Bioscience, Luckenwalde, Germany) following the manufacturer's protocol and results were qualitatively assigned as being positive or negative for Giardia. In addition, samples were used for a semi-quantitative assignment of cyst numbers in the microscopic sample, which approximates 10 μl of volume: +, 1-10 cysts; ++, 11-50 cysts; +++, > 50 cysts.

DNA isolation
Genomic DNA was extracted from samples using one of the two commercial kits following the protocols of the manufacturers: QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany) or Maxwell® 16 FFPE Plus LEV DNA Purification Kit (Promega, Mannheim, Germany). The final elution volume was 70 μl. Both kits have been tested and performed similarly in subsequent PCR applications. The amount of nucleic acid in samples was photometrically estimated at OD 260 .
qPCR All samples were tested using a TaqMan-based Giardia specific qPCR as earlier described [27,43] with minor modifications. Briefly, a 62 bp fragment of SSU rDNA was amplified in a total PCR volume of 25 μl Titration experiments with DNA from G. duodenalis trophozoites (WB6, ATCC 50803) revealed a detection limit (analytical sensitivity) in our technical setting of approximately one trophozoite (equals approximately four genome equivalents) in the PCR. Analytical specificity was calculated as 98% (calculated from 49 control PCRs of which one was slightly positive with a Ct-value of 38.5). Inhibition was not found to be an issue in the PCR analyses (Only one of 35 tested samples showed a slight inhibition with a delta-Ct value > 2 in a test set with internal amplification control).
We also used the Ct-value of the qPCR analysis as a proxy for parasite numbers. We assumed that target DNA integrity is similar in the samples' fecal matrix of all three rodent genera and that PCR efficacies were equal, although we were not able to formally test these assay parameters.

Typing PCR and sequencing
A fragment of the SSU rRNA gene was amplified using a combination of previously described and modified PCR protocols [4,28]. A nested PCR protocol was used with the initial primer pair RH11-derivates (equal mix of 5'-CAT CCG GTC GAT CCT GCC-3' and 5'-CAT CCG GTT GAT CCT GCC-3') and Gia2150c Fragments of the bg and gdh genes were amplified according to previously described nested-PCR protocols [5,7]. Briefly, PCRs were run in a total volume of 50 μl and included 1-2 μl target DNA, 200 μM dNTPs, 1× PCR MangoTaq buffer containing 3 mM MgCl 2 (Bioline), 2.5 U of MangoTaq polymerase (Bioline), and 200 nM of each primer. The reactions were performed for 35 cycles using following conditions: 1st PCR (94°C for 45 s, 65°C (for bg) and 56°C (for gdh) for 30 s and 72°C for 45 s) and nested-PCR (94°C for 45 s, 65°C (for bg) and 56°C (for gdh) for 30 s and 72°C for 45 s). Final extension was done for 7 min at 72°C.
Some of the obtained sequences were analyzed to identify the most similar sequence deposited in the GenBank public database using the built-in Basic Local Alignment Search Tool (BLAST) algorithms (http:// blast.ncbi.nlm.nih.gov/Blast.cgi). Sequence names are deduced from a unique running sample number preceded by an abbreviation of the sequenced genes (e.g. ssu349 = SSU rDNA sequence of rodent isolate number 349). Nucleotide sequences generated in this study have been deposited in the GenBank database under accession numbers KY114167-KY114486 (SSU rDNA) and MG676938-MG677109 (bg and gdh). Accession numbers of each sequence are listed in Additional file 4: Table S2.

Phylogenetic analysis
Alignments of completely sequenced PCR products (without primer sequences) and respective regions of the reference sequences were produced using MUSCLE [44] integrated subroutine of Geneious version 9.1. (Biomatters Ltd.). Sequences containing polymorphic positions, a wellknown phenomenon in G. duodenalis [45], were included in the analysis when the sequencing data of both strands was available. Bayesian analysis of sequence alignments was performed by using the BEAUti and BEAST software packages [46]. PhyML analysis was done with ATGC online tool PhyML 3.0 [47] using the Smart Model Selection tool to select the best substitution model and subsequent PhyML analysis using the best substitution model with Subtree-Pruning-Regrafting (SPR) tree searching and bootstrap performance of 100 [48]. Trees were annotated using the iTOL online tool [49,50]. Phylogenetic trait analysis was done using the program BaTS [51].

Statistical analysis
Data were organized in a spread sheet and subsequently imported into the statistics software package STATA 14. 1 (StataCorp, College Station, USA). Proportions were calculated and analyzed for binomial exact 95% confidence intervals (95% CI). To test for any difference between proportions of groups the Fisher-Freeman-Halton test was used [52]. For comparison between groups, the two-sided Fisher's exact test was used and a multiple correction using the Bonferroni-Holm procedure was applied and adjusted P-values are reported [53]. A P-value of ≤ 0.05 was considered statistically significant. Due to the low number of cases we have chosen not to use multivariate analyses using logistic regression models. For qPCR data (Ct values) the median and min/max values of groups were calculated and non-parametric Kruskal-Wallis test was performed followed by Dunn's test of multiple comparisons to assess statistical significance using the software package GraphPad Prism 7.03 (a P-value of ≤ 0.05 was considered statistically significant).