ECO-EPIZOOTIOLOGIC STUDY OF FRANCISELLA TULARENSIS, THE AGENT OF TULAREMIA, IN QUÉBEC WILDLIFE

Tularemia is a bacterial zoonotic disease caused by Francisella tularensis. During 2003–11, 93 human cases were reported in Canada, of which over half were acquired in Québec (Institut National de Santé Publique du Québec 2011; Public Health Agency of Canada 2014). In Canada, human tularemia has been mainly associated with contact with snowshoe hares (Lepus americanus) and muskrats (Ondatra zibethicus; Wobeser et al. 2009). Two clinically predominant subspecies are recognized: type A, or Francisella tularensis tularensis, occurs in North America, and type B, Francisella tularensis holartica, occurs throughout the Northern Hemisphere (Keim et al. 2007). Type A is divided into subpopulations, type AI and type AII, characterized by different virulence and distinct geographic divisions, with type AI being more pathogenic and established in eastern North America (Staples et al. 2006; Molins et al. 2010).

Natural infections have been reported in over 200 animal species, each manifesting different degrees of susceptibility to the disease (Morner and Addison 2001). Lagomorphs and rodents are believed to be the main hosts and hematophagous arthropods, notably ticks, can play a substantial role in pathogen transmission (Jellison et al. 1961). Two main ecologic cycles have been described (Eliasson et al. 2006; Telford and Goethert 2010). The terrestrial cycle, principally associated with type A, involves essentially rabbits (Sylvilagus spp.), hares (Lepus spp.), and ticks. The aquatic cycle is typically linked to type B, in close association with aquatic environments, and involves predominantly beavers (Castor canadensis) and arvicoline rodents, such as muskrats. Exposure to the bacterium is also frequent in carnivores, including coyotes (Canis latrans). Dietary intake is potentially an important route of acquisition for coyotes because they frequently feed on hares and rodents, including snowshoe hares and muskrats (Centre d'Expertise en Analyse Environnementale du Québec 2006). Evidence suggests that coyotes are resistant to clinical disease (Gier and Ameel 1959), rendering them possible sentinels for occurrence of natural transmission cycles.

The transmission cycle of F. tularensis is complex and not clearly understood, and the primary natural reservoir species and ecologic determinants of persistence and spread remain unidentified. This study was undertaken to develop a better understanding of the ecology of the infection in wildlife. Our objectives were to 1) estimate the prevalence of F. tularensis DNA or antibody to F. tularensis in snowshoes hares, muskrats, and coyotes in four administrative regions of Québec, Canada, 2) describe the risk of infection regarding host characteristics and environmental features as they relate to the different ecologic cycles, and 3) detect spatial clusters of positive animals, which may indicate heterogeneities in environmental risk of exposure to F. tularensis.

We conducted this study in four administrative regions in the province of Québec (Abitibi-Témiscamingue, Bas-Saint-Laurent, Chaudière-Appalaches, and Estrie; Fig. 1), Canada, where coyotes, muskrats and snowshoe hares are frequently trapped. Abitibi-Témiscamingue has an abundance of lakes and includes mostly boreal forest in the north and mixed forest in the south with some agricultural areas. Bas-Saint-Laurent and Chaudière-Appalaches feature coastal plains along the Saint Lawrence River with low-lying agricultural areas bordered by the Appalachian Mountains to the south. They contain some major rivers with a few lakes and predominantly coniferous and mixed forest in Bas-Saint-Laurent and broadleaf forest in Chaudière-Appalaches (Ministère des Ressources Naturelles et de la Faune 2010a, b, c). Estrie is a small region in the Appalachian Mountains that features 75% forest (principally mixed and broadleaf forest), 17% agricultural area, and 3% open water (Conférence Régionale des Élus de l'Estrie 2010).

The fieldwork for this cross-sectional study, involving snowshoe hares, muskrats, and coyotes, was carried out October 2012–April 2013. We calculated a target sample size of 100 carcasses per species per region to detect at least one positive animal within each subgroup with 95% confidence, assuming a prevalence of ≥3%. We recruited trappers on a volunteer basis through local trappers' associations and made every effort to cover a broad area of each region. Carcasses were collected from regular trapping activities. For coyotes, whenever possible, trappers collected a sample of blood, the left kidney with surrounding fat, and the mandible after the animal was skinned. Otherwise, the carcass was sent directly to the Université de Montréal, Québec, Canada, for processing. Muskrat carcasses were sent after skinning by trappers; snowshoe hare carcasses were sent intact. Trappers were asked to keep all specimens and carcasses frozen until shipment to the laboratory. For each carcass collected, trappers recorded the GPS coordinates or the address closest to the capture site.

We kept carcasses at approximately 4 C until completely thawed prior to processing. For each animal, a blood sample was collected from the abdominal, thoracic, or heart cavity and stored at −20 C prior to analyses. We collected samples of kidney, liver, lung, and spleen from snowshoe hares and muskrats, which were stored individually in sterile plastic bags at −80 C.

To evaluate the body condition of coyotes, we calculated the kidney fat index as the ratio of the kidney weight with the surrounding fat (including the fibrous capsule) to kidney weight after trimming surrounding fat, using a 1-g precision scale (Huot et al. 1995; Winstanley et al. 1998). The index was averaged for the two kidneys, if both were collected. For snowshoe hares, we evaluated body condition by visual evaluation of the quantity of abdominal fat surrounding the kidneys by using a scoring system from 0–2 (0: emaciated, no fat surrounding the kidney, 1: average, fat partially surrounding the kidney, and 2: very good, fat covering the whole kidney). Kidney fat provides an index of nutritional status in snowshoe hares (Hodges et al. 2006). We did not evaluate body condition for muskrats because we found no reliable method in the literature.

We classified coyotes as juveniles (<1 yr), subadults (1–2 yr), or adults (>2 yr) according to cementum age analysis of the lower canine at Matson's Laboratory LLC, Milltown, Montana, USA (Jean et al. 1986). We did not estimate the age of muskrats and snowshoe hares because no positive results were obtained and no reliable method was found in the literature, respectively.

We thawed and then centrifuged blood samples (1,500 × G for 5 min). The supernatant was placed in sterile microtubes and shipped on ice to the National Microbiology Laboratory in Winnipeg, Manitoba, Canada. Upon receipt, sera were recentrifuged, and the supernatants were tested by a microagglutination test, as described by Sato et al. (1990). We considered antibody titers ≥64 indicative of past infection and referred to those animals as positive (see Discussion section for the rationale of this cutoff). Additionally, we tested coyotes with titers ≥16 for antibodies against Brucella abortus by similar agglutination tests to rule out a potential cross-reaction between F. tularensis and B. abortus (Behan and Klein 1982).

We tested snowshoe hares and muskrats for F. tularensis by using real-time PCR on a pool of kidney, liver, lung, and spleen samples for each animal. In addition, based on serologic results, we performed supplementary PCR on individual organs of the two antibody-positive snowshoe hares and 30 hares collected nearest to these animals. PCR was also performed on one kidney of each coyote (n=46) with antibody titers ≥16. Each of these coyotes was paired with the nearest antibody-negative coyote whose kidney was also tested for F. tularensis by PCR.

We extracted DNA from approximately 0.1 cm³ of each sample by using DNeasy Blood and Tissue Kits (Qiagen, Hilden, Germany) following the manufacturer instructions. The DNA was subject to an in-house developed real-time PCR assay targeting fopA, and any positive results were confirmed via PCR for two unique signature sequences. Subsequently, multiple locus variable number tandem repeat analysis (MLVA) was conducted for positive samples by using primers and probes, as described by Vogler et al. (2009a, b). We analyzed the MLVA electrophoretic data, as described by Antonation et al. (2015).

A sufficient number of positive results to perform risk factor analysis were obtained only for coyotes. We used ArcGIS version 10.2 (Esri, Redlands, California, USA) to extract environmental features in three defined circular buffers with radii of 2, 4, and 6 km around each coyote capture site corresponding to three scales of possible home ranges (respectively, 12.5 km², 50 km², 113 km²; Messier and Barrette 1982; Crete et al. 2001). We used 1) vegetation covers from the Ecoforestery Information System (Ministère des Forêts, de la Faune et des Parcs of Québec, Québec, Canada), 2) wetland covers from Ducks Unlimited (Ducks Unlimited Canada, Stonewall, Manitoba, Canada), and 3) water bodies maps from CanMap Water (DMTI Spatial Inc., Markham, Ontario, Canada) to define five main land cover types: forest, shrubland-grassland, agriculture, wetland, and waterside. For the latter, we calculated an external buffer of 10 m around each inland water body to represent access to the water for the coyote. We calculated the proportion of each land cover type within three spatial scales, using the area of the buffer, excluding water bodies.

We assessed the availability of suitable habitat for snowshoe hare in each buffer by using the habitat suitability model for this species in Québec (Forêt Modèle Bas-St-Laurent et Université du Québec à Rimouski 2003). This model uses information from the Ecoforestry Information System database to classify forest stands in terms of their capacity to offer forage and cover to snowshoe hares (high=0.75, medium=0.5, low=0.25, and none=0). We also calculated an area-weighted average of habitat suitability for each buffer to give a global index of habitat suitability for hares. We calculated the proportion of high-level habitat suitability and the global index of habitat suitability for each spatial scale, as described for the land cover variables. For the GPS position of each coyote, we obtained the mean temperature and precipitation from an interpolated grid (from 2002 to 2010) from the National Land and Water Information Service of Canada (McKenney et al. 2011).

We estimated prevalence of PCR-positive and antibody-positive animals with 95% exact confidence limits for each species. We used exact logistic regression models, performed in SAS software version 9.3 (SAS Institute Inc., Cary, North Carolina, USA), to model the antibody status of coyotes according to animal level and environmental characteristics (Tables 1, 2). The exact logistic regression is recommended for analyzing sparse datasets, such as the one generated from this study (Derr 2000). Based on our proposed causal Web diagram (Fig. 2), two models were considered that included variables from individual characteristics, water-related land, climatic conditions, as well as vegetation cover for model “I” and potential hare habitat for model “II.” If the linearity assumption was not respected based on a graphic method (Dohoo et al. 2003), we categorized continuous variables using the 25th and 75th percentiles (from the 4-km buffer) as cutoffs. We considered all explanatory variables with P values ≤0.20 from univariate analyses eligible for inclusion in the multivariable logistic regression analysis (Tables 1, 2), with the exception of potential confounders, which we included in the full models, regardless of the P value. We calculated pairwise Pearson correlations among explanatory variables selected from univariate analyses on the continuous scale. If a strong correlation (|r|>0.7) was detected, we included only one variable in the full models based on the strongest biologic plausibility. We obtained the final models using a backward selection procedure with P>0.15 (exact conditional analysis) as the rejection criterion, allowing for liberal inclusion of variables. We kept potential confounding factors in the final models if their removal would change the coefficient of other predictors of interest by more than 30% or if their P value was ≤0.15.

When comparing the three spatial scales, the best performing scale was identified for model I and model II by comparing two performance metrics: 1) the area under the curve (AUC) of the receiver-operator characteristic curve and 2) the model classification accuracy when using the sensitivity-specificity equivalence point as a cutoff value for classification (Caraguel et al. 2011). We calculated classification statistics by using the ROCR package version 1.0-5 in R version 3.1.2 (R Development Core Team 2014). For the variables kept in models I and II comprising the highest AUC, the impact of the potential imperfect sensitivity of the diagnostic tests on the odds ratio estimates was assessed in deterministic sensitivity analyses performed in two-way contingency tables, keeping specificity at 1.00 or 0.99 (Fox et al. 2014).

We assessed the presence of spatial clusters of antibody-positive coyotes by using the Kulldorff spatial scan test in SatScan version 9.2 (Kulldorff 1997). We used the Bernoulli model with a maximum cluster size of up to 50% of the population at risk, with no geographic overlap in clusters permitted. The significance (α<0.15) of clusters was determined through 999 Monte Carlo replicates. Cluster analyses were conducted separately for Abitibi-Témiscamingue given that this region is remote from the three others.

We collected 1,141 carcasses in the four regions during the trapping season, October 2012–April 2013 (Fig. 1), which included 385 coyotes (190 female, 190 male, and five unknown), 345 snowshoe hares (162 female, 177 male, and six unknown), and 411 muskrats (199 female, 201 male, and 11 unknown). The body condition of snowshoe hares was 151 below average, 115 average, and 65 very good (14 missing). Individual characteristics of coyotes are presented in Table 1. Missing data were due to lack of assessment, sometimes because of poor conservation or missing parts of carcasses.

Anti-F. tularensis antibodies were detected only in 11 coyotes and two hares (Fig. 1; Table 3; see Supplementary Table S1 for titers). In contrast, antibodies against B. abortus were not detected (titer <20) in the 46 samples of coyotes that had F. tularensis antibody titers ≥16.

All tested organ pools of muskrats (n=411) and snowshoe hares (n=345) were PCR negative, as were the kidneys of all 46 coyotes with antibody titers ≥16 and of 46 location-matched antibody-negative coyotes (median distance of 0.24 km [range=0–57 km] between each location-matched coyotes). However, DNA of F. tularensis was detected in spleen or lung tissues processed individually (i.e., not pooled) from the two antibody-positive snowshoe hares collected in Bas-Saint-Laurent. All individual organs of the other 30 hares located closest to these two animals (median=17 km, range = 0–28 km) were PCR negative. The MLVA patterns of the two positive hares clustered with type AI isolates within the National Microbiology Laboratory in-house database. The two patterns differed by only one of the 11 loci used, locus M03, a marker well-known to mutate rapidly (diversity D=0.95; data not shown).

In the risk factor analysis for presence of antibodies in coyotes, all continuous variables did not respect the linearity assumption and were therefore categorized. The age, forest, agriculture, wetland, high-level habitat suitability for hares, and global index of habitat suitability for hares variables were selected from univariate analyses for at least one spatial scale (Tables 1, 2).

Collinearity was detected between forest and agriculture variables, as well as the high-level habitat suitability and global index of habitat suitability variables. Forest and high-level habitat suitability were selected for inclusion in the full models. Higher odds of detecting antibodies were observed in subadult coyotes compared with juveniles in models I and II (Tables 4, 5), in areas with a higher percentage of forest in model I (Table 4) and in areas with higher suitability for hares in model II (Table 5). Temperature was retained as a potential confounder only for the model II using a radius of 4 km (Table 5).

Very similar results from the AUC and the sensitivity-specificity approach were obtained in models I and II, making it difficult to distinguish the best performing scale of analysis (Tables 4, 5). Furthermore, we could not perform the sensitivity analyses for misclassification on model I because of the presence of zero cells. Imperfect sensitivity of the diagnostic test for the model II comprising the highest AUC (radius of 4 km) was associated with underestimation of the odds ratios for age and high-level habitat suitability by <5% for the tested range of value (Table 6).

No spatial clusters of antibody-positive coyotes were detected (P>0.54; Fig. 1). Cluster analysis was not conducted for Abitibi-Témiscamingue because only one coyote was antibody positive.

The low antibody prevalence in snowshoe hares was consistent with antibody prevalences (0% to 1.6%) reported in various provinces of Canada (Hoff et al. 1970; Zarnke and Yuill 1981; Akerman and Embil 1982; Cayouette 1993). Higher antibody prevalences in snowshoe hares were described in two studies in Ontario, one of these following an outbreak in muskrats (Mac Lulich 1937; Ditchfield et al. 1960). The low antibody prevalence frequently seen in snowshoe hares could be due to infrequent exposure to F. turalensis or high susceptibility to disease. Many authors have suggested that most snowshoe hares, unlike cottontail rabbits, are unaffected by the infection (Green 1938, 1939; Jellison and Parker 1945), but significant mortality in snowshoe hares occurs occasionally during tularemia outbreaks and in experimentally infected animals (Jellison et al. 1961; Miller 1974). Our detection of F. tularensis DNA in organs from two antibody-positive snowshoe hares suggests that infected snowshoe hares can develop antibodies and survive infection even with the most virulent type A1 strain.

No evidence of past or active infection to F. tularensis was found in muskrats, as in previous studies conducted in endemic areas in Saskatchewan, Canada (F.A. Leighton pers. comm.), and in Germany (von Keyserlingk et al. 2011). However, a prevalence of 5% by culture and 41% by agglutination tests was reported in muskrats in Vermont following a large outbreak of tularemia in humans (Young et al. 1969). Infection in muskrats usually results in sudden death; affected populations show high mortality and, during outbreaks, infection is detectable in dead animals (McDermid 1946; Parker et al. 1951; Fyvie et al. 1959; Ditchfield et al. 1960; World Health Organization 2007). This could explain the absence of antibody-positive muskrats in our study, but it is also possible that muskrats were not exposed to F. tularensis if the bacterium was not present in an aquatic cycle in their habitat.

We detected anti-F. tularensis antibodies in 2.9% of coyotes. In the US, reported antibody prevalences range from 0% to 32% or more (Thorpe et al. 1965; Trainer and Knowlton 1968; Gier et al. 1978; Gese et al. 1997; Arjo et al. 2003; Gese et al. 2004; Bischof and Rogers 2005; Chronert 2007). We did not detect infection by PCR in kidneys of any of the 92 coyotes tested, including the 46 coyotes with antibody titers ≥16, which, to our knowledge, has not been investigated previously in free-ranging coyotes. The detection of antibodies coupled with the absence of evidence of active infection supports a commonly held belief that coyotes are resistant to tularemia. In fact, apart from very young pups, experimentally infected coyotes recover from infection and, additionally, are not likely to disseminate the bacteria in the environment (Stagg et al. 1956; Lundgren et al. 1957).

Subadult coyotes were more likely to be antibody-positive than juveniles, which has not been seen previously (Gese et al. 1997; Arjo et al. 2003; Gese et al. 2004). These results may be the reflection of longer exposure in subadults compared with younger coyotes. Moreover, in experimental infections, very young coyotes occasionally die from F. tularensis infection (Parker and Francis 1926; Marchette 1960), which could result in a disproportionate removal of exposed young vs. subadult coyotes. Sex was not associated with antibody detection in coyotes, as previously reported (Gese et al. 2004); neither was body condition, suggesting no or limited lasting clinical consequences of infection.

The results from model I for all spatial scales suggest that forests have ecologic characteristics that support F. tularensis transmission. Wooded areas are suitable for species of ticks in Québec that can be vectors of tularemia, such as the American dog tick (Dermacentor variabilis; Bishopp and Trembley 1945; Hopla 1974; Artsob et al. 1984). Alternatively, model II suggests, for all spatial scales, that snowshoe hares might play a significant role in the infection of coyotes, perhaps due to their prey-predator relationship. The coyote that was collected closest to the two antibody-positive hares (16 km away) was also antibody positive, supporting the possibility of transmission between these two species. Rabbit ticks (Haemaphysalis leporispalustris), a potential vector of F. tularensis, frequently infest hares and rabbits (Bishopp and Trembley 1945; Telford and Goethert 2010). According to a study conducted in Québec, all snowshoe hares examined in the summer of 1992 were heavily and exclusively infested with these ticks, suggesting a high potential for their involvement in F. tularensis transmission (Cayouette 1993). Thus, snowshoe hares associated with rabbit ticks may constitute a potential polyspecies reservoir of F. tularensis, as previously proposed (Hopla 1974; Cayouette 1993).

We detected no spatial clusters of antibody-positive coyotes. This is in contrast with the general belief that tularemia persists in natural foci (Pavlovsky 1966), as supported by studies in which clusters were detected (Goethert and Telford 2009; Svensson et al. 2009). However, we cannot rule out the possibility that each positive coyote was linked to a separate small natural focus that occurred in the past or that the geographic location, where the coyote was captured, was remote from the location of exposure to the infection, because coyotes can roam long distances.

Our study has some limitations inherent to the cross-sectional design, including the inability to determine the timing of bacterial infection. Second, lack of statistical power due to low antibody prevalence likely limited our chances of detecting significant associations and spatial clusters. Also, an opportunity sampling was used in this study, which might have resulted in the selection of nonrepresentative zones of sampling by trappers. Likewise, within trapping areas, the potential higher probability of trapping young animals and the potential exclusion of moribund animals by trapping methods may have underestimated the prevalence of antibody and prevalence of infection, respectively. Furthermore, blood samples were partially diluted with body fluids, and organs were sometimes partially autolyzed, which may have reduced the sensitivity of detection tests. Pooling of organs for PCR may have decreased the probability of detecting F. tularensis DNA, as was observed for our two antibody-positive hares, possibly because all organs were not equally infected. However, this potential lack of test sensitivity gave conservative estimates for the risk factors according to simulations and therefore was mostly a concern in prevalence estimation. We used an antibody titer cutoff of ≥64, which is lower than the ≥128 used to diagnose an infection in humans. However, while no cutoff has been validated in wildlife species, F. tularensis DNA was detected in a hare with an antibody titer of 64, which supported use of 64 rather than a higher threshold for detecting serologic responses in early infection, as well as in recovered animals. This was also supported by the likely absence of cross-reactions with the anti-Brucella antibody. Finally, this study was essentially exploratory, justifying our choice of using a liberal statistical significance (P≤0.15) for our analysis. Therefore, the associations found in our study must be seen as preliminary and need to be validated by further studies.

Our results collectively suggest that a terrestrial cycle is present in the area studied: positive results in snowshoe hares and coyotes, identification of the strain type AI, and forest or high level of hare habitat suitability as risk factors. This is also supported by a recent report of type AI in 78% of humans and 100% of hare isolates from 1998–2011 within Québec (Antonation et al. 2015). Furthermore, we found no evidence of the aquatic cycle, with all muskrats found negative, and none of the water-related variables being statistically relevant in our models. Despite our results, the existence of the aquatic cycle in Québec cannot be ruled out. In fact, type B was also isolated in humans in areas of high muskrat density in Québec (Antonation et al. 2015). Further studies should be conducted to evaluate the possible existence of the aquatic cycle, as well as the role of snowshoe hares and ticks in the terrestrial cycle of F. tularensis in Québec. Moreover, the use of trapped coyotes for the surveillance of F. tularensis circulation in the environment could be further evaluated.

Funding was provided by the Public Health Agency of Canada, the Ministère des Forêts, de la Faune et des Parcs, and the Research Group on Epidemiology of Zoonoses and Public Health. This work was also made possible through graduate student research awards from the Fonds de recherche du Québec–Nature et technologies, Lucie Besner grant, Université de Montréal (Faculty of Veterinary Medicine and Faculty of Graduate Studies) attributed to V.G.-R. We thank Pierre Canac-Marquis, Denise Bélanger, Benoit Lévesque, Thérèse Lanthier, Andres Ivan Ramirez, Josée Guillemette, Jessie Longpré, Sonia Bernier, Huguette Mallet, the Matson's Laboratory LLC, the Centre Québécois sur la santé des animaux sauvages, the National Land and Water Information Service of Canada, and all trappers for their important contribution to fieldwork.

Supplementary material for this article is online at http://dx.doi.org/10.7589/2015-04-096.