Risk maps for range expansion of the Lyme disease vector, Ixodes scapularis, in Canada now and with climate change
- Nicholas H Ogden1, 2Email author,
- Laurie St-Onge3,
- Ian K Barker4,
- Stéphanie Brazeau3,
- Michel Bigras-Poulin2,
- Dominique F Charron5,
- Charles M Francis6,
- Audrey Heagy7,
- L Robbin Lindsay8,
- Abdel Maarouf9,
- Pascal Michel3,
- François Milord10,
- Christopher J O'Callaghan11,
- Louise Trudel12 and
- R Alex Thompson13
© Ogden et al. 2008
Received: 22 March 2008
Accepted: 22 May 2008
Published: 22 May 2008
Lyme disease is the commonest vector-borne zoonosis in the temperate world, and an emerging infectious disease in Canada due to expansion of the geographic range of the tick vector Ixodes scapularis. Studies suggest that climate change will accelerate Lyme disease emergence by enhancing climatic suitability for I. scapularis. Risk maps will help to meet the public health challenge of Lyme disease by allowing targeting of surveillance and intervention activities.
A risk map for possible Lyme endemicity was created using a simple risk algorithm for occurrence of I. scapularis populations. The algorithm was calculated for each census sub-division in central and eastern Canada from interpolated output of a temperature-driven simulation model of I. scapularis populations and an index of tick immigration. The latter was calculated from estimates of tick dispersion distances by migratory birds and recent knowledge of the current geographic range of endemic I. scapularis populations. The index of tick immigration closely predicted passive surveillance data on I. scapularis occurrence, and the risk algorithm was a significant predictor of the occurrence of I. scapularis populations in a prospective field study. Risk maps for I. scapularis occurrence in Canada under future projected climate (in the 2020s, 2050s and 2080s) were produced using temperature output from the Canadian Coupled Global Climate Model 2 with greenhouse gas emission scenario enforcing 'A2' of the Intergovernmental Panel on Climate Change.
We have prepared risk maps for the occurrence of I. scapularis in eastern and central Canada under current and future projected climate. Validation of the risk maps provides some confidence that they provide a useful first step in predicting the occurrence of I. scapularis populations, and directing public health objectives in minimizing risk from Lyme disease. Further field studies are needed, however, to continue validation and refinement of the risk maps.
Lyme disease is an emerging infectious disease in Canada, largely due to the expansion of the geographic range of the tick vector Ixodes scapularis in eastern and central Canada . Up to 1997, the only population of I. scapularis known in Canada was that at Long Point, Ontario . However, the number of known I. scapularis populations in Canada has risen from 1 to 13  in the last decade, a period that may have witnessed the first evidence of global warming [3, 4]. Global warming is anticipated to accelerate expansion of the geographic range of I. scapularis into Canada, provided that suitable habitat and hosts occur , and may influence the emergence of tick-borne zoonoses .
If we are to limit the impact of emerging Lyme disease on human health in Canada, the appropriate public health messages and clinical information must be issued to the public and medical practitioners and targeted to the right risk locations . Early diagnosis and treatment of Lyme disease is usually simple and successful, but diagnosis and treatment of later stages of systemically disseminated Lyme disease is much more difficult and costly to the patient and health services . Furthermore, targeted surveillance for I. scapularis ticks and Borrelia burgdorferi (the agent of Lyme borreliosis) is needed to identify endemic locations, as this knowledge is important in assisting clinical diagnoses . For vector-borne diseases such as Lyme disease, the human populations at risk are at least in part defined by the geographic occurrence of the arthropod vectors, whose existence is tightly linked to climatic variables on a continental scale [9, 10], and to suitable habitat on a more local geographic scale [11, 12]. Predicting the occurrence of vectors by risk maps is, therefore, a potentially very useful tool to guide public health policy and target surveillance and intervention activities for vector-borne diseases [13, 14].
In many cases, potential geographic distributions of vectors have been predicted using statistical associations between climatic or landscape variables (or their remote-sensed proxies), which are likely associated with vector survival and/or reproduction, and observed occurrence of vectors or the diseases they transmit [11–17]. However, with a greater depth of understanding of how climate and landscape variables directly affect biological process of vector survival and abundance, these relationships alone can be used to develop risk maps, a process which can be particularly useful for predicting expansion of the geographic ranges of vectors and vector-borne diseases .
For I. scapularis, we have considerable empirical information on how climate, particularly temperature, can affect the survival of the ticks, and this information has already been synthesised in a simulation model to identify relationships between temperature and I. scapularis population survival . We also have assessed empirically how habitat in southeastern Canada may affect I. scapularis survival . Furthermore, we have information on the routes of dispersion of I. scapularis from locations where the tick is endemic at present. Recently, we identified northward migrating land birds as potentially important in expanding the range of I. scapularis for three reasons: i) a high number (circa 2%) of the estimated 3 billion birds migrating northward into and through eastern and central Canada carry an I. scapularis tick , ii) these migratory birds are capable of surmounting geographic features (lakes, sea, mountains and areas of intensive agriculture) that are obstacles to dispersal by terrestrial hosts, and iii) the long distances that migratory birds can fly during the 2–4 days that a larval or nymphal I. scapularis feeds on its host . Furthermore, some of the ticks carried by these birds are infected with tick-borne pathogens including several genotypes of B. burgdorferi and Anaplasma phagocytophilum (the agent of Human Granulocytic Anaplasmosis) . These infected ticks pose an immediate, but low-level threat to human health in Canada after dropping off the birds and moulting into the next instar [1, 23]. They could infect wildlife with these pathogens thus introducing endemic cycles of infection into newly established reproducing populations of the tick vector, which would give rise to a greater public health risk . This may, however, be a process that takes some years . Most I. scapularis ticks carried by migratory birds are nymphs, which moult into adults that feed mostly on reservoir-incompetent deer, and rarely on competent reservoir host species .
In this study, therefore, we aimed to develop risk maps for the occurrence of the Lyme disease vector I. scapularis in Canada under current climate using information on how climate and habitat may affect tick survival, and how current geographic distributions may affect tick dispersion on migratory birds. We then aimed to investigate how the geographic distribution of I. scapularis may change under possible future climate conditions.
Throughout the study, the geographic unit was the census sub-division (CSD) whose geographic area, geographic coordinates and population information were obtained as a coverage and spreadsheet from Statistics Canada . This spatial unit is increasingly used in Canadian public health because i) it is a well defined spatial unit available to users, including provincial and municipal public health organisations, ii) its spatial boundaries are relatively stable over time, and iii) accompanying data on human population acts as a denominator for surveillance activities .
Selection of risk algorithms and cut-off levels
Algorithms used to predict the occurrence of Canadian census sub-divisions containing resident I. scapularis populations.
No. of ticks at model equilibrium (T)
0.863 – 0.979
No. of ticks at model equilibrium categorised (Tc)
0.678 – 0.881
Percent forest area (F)
0.313 – 0.461
Index of larval tick immigration (range 255 km: I L )
0.713 – 0.919
Index of nymphal tick immigration (range 425 km: I N )
0.848 – 0.949
T * I N
0.869 – 0.983
Tc * I N
0.699 – 0.914
T * I L
0.743 – 0.947
Tc * I L
0.614 – 0.832
T * I N * (0.05* I L )†
0.869 – 0.983
Tc * I N * (0.05* I L )†
0.699 – 0.914
T * I N * F
0.723 – 0.919
Tc * I N * F
0.646 – 0.858
T * I N * Log10 F
0.732 – 0.933
Tc * I N * Log10 F
0.652 – 0.871
Validation of index of nymphal tick immigration
Validation of the nymphal tick immigration index.
Index of tick immigration
Newfoundland & Labrador
Prince Edward Island
There was almost complete confounding between the index of nymphal tick immigration for each CSD in southern Quebec and the latitude of the CSD, a risk factor for ticks submitted by passive surveillance that was found significant in a previous study . In a minimal multivariable model in which Ln human population was an offset, average spring normalised difference vegetation index (NDVI: a second risk factor for I. scapularis identified by Ogden et al ) and the index of nymphal tick immigration were significantly associated with the number of I. scapularis submitted by the public (coefficients = 1.845 and 0.063, SE = 0.48 and 0.01, Wald Z = 3.8 and 8.7 respectively, P < 0.001 for both). AIC for models containing latitude and the index of nymphal tick immigration were 1362 and 1366 respectively, and with both variables together in the model, AIC was 1365, i.e. the index of nymphal tick immigration and latitude were almost totally confounded, supporting the hypothesis in Ogden et al  that latitude acted as a proxy for tick dispersion on migratory birds or other animal hosts of the tick.
Risk maps under current and projected climate
In this study, we developed an algorithm for predicting the occurrence of I. scapularis populations in Canada. The algorithm predicted with reasonable accuracy the occurrence of the small number of I. scapularis populations that were known in Canada prior to this study. A warming climate was considered to cause changes in the speed of bird migration and distances of dispersion of ticks by migrating birds, and changes in the geographic range of territory with a climate warm enough for the ticks to complete their lifecycle and establish populations. With this information, the algorithm was then used as a simple model to predict the progression of I. scapularis as it spreads across Canada as the climate warms, given two scenarios for the current potential occurrence of I. scapularis populations in Canada.
There were of course many assumptions in this process. These include i) the assumptions of the I. scapularis population model , which have more effect on the gradient of the relationship between tick abundance and DD > 0°C than on the intercept with the x-axis ; ii) assumptions of the climate models and the emissions scenario; iii) assumptions that the interpolation of temperature data from meteorological stations and, in particular, output from the climate model, which have a crude 400 km spatial scale, have a negligible effect; iv) assumption that the distribution of I. scapularis-endemic areas in the USA described by Dennis et al.  represents the true distribution, and that the distribution of known I. scapularis populations in Canada represents the true distribution; v) assumption that the speed at which birds migrate is equal across Canada irrespective of water bodies or other hazards for the birds (which is possibly relevant to the Atlantic provinces but unlikely to impact dispersion elsewhere in Canada); and vi) assumption that there are no other determinants of I. scapularis establishment (e.g. host densities) that are more limiting than temperature and rates of dispersion by hosts. The latter is relevant to Prince Edward Island and Newfoundland and Labrador where there are no white-tailed deer, which are thought to be key hosts for adult I. scapularis .
Our validation of the index of nymphal tick immigration suggests that the distribution of I. scapularis-endemic sites in Dennis et al  combined with the estimate of bird dispersion ranges from these and Canadian endemic areas provides a workable model of the spatial dispersion of I. scapularis across Canada. Furthermore, the field study in southern Quebec suggested that the risk algorithm predicts the occurrence of I. scapularis surprisingly well in an area where, prior to our study, reproducing I. scapularis populations were not thought to exist. The selected risk algorithm performed better at predicting the occurrence of potential I. scapularis populations than either of its two component parts individually (model-derived relationships between temperature conditions and tick abundance, and an index for predicted numbers of immigrating ticks), i.e. there appeared to be synergy between these values. This field validation, combined with validation of the index of nymphal tick immigration, gives confidence that the maps may give us insight into the possible range of I. scapularis at present, that temperature is a likely limiting factor on I. scapularis distribution, and that climate change could drive changes to the geographic footprint of this tick and possibly Lyme disease risk in the coming decades. The prevalence of I. scapularis positive sites in both high risk and moderate risk zones were lower than the highest estimates for prevalence of positive CSDs predicted by algorithm cut-off values. This difference could be expected as absence of I. scapularis from one site per CSD does not confirm complete absence of I. scapularis from that CSD, and because at the edge of range expansion it would be expected that not all I. scapularis-suitable niches are filled by I. scapularis populations.
Clearly, the further we project into the future the distributions of I. scapularis populations, the less certain we are of them because of uncertainties in emissions scenarios, climate model outputs and rates of tick dispersion. However the risk maps as presented represent both an architecture and a starting point for decision making on policies on, and targeting of, surveillance activities. Surveillance data should then be used to continue validation of the risk maps and drive refinements due to changes in the algorithm used here, or by addition of additional information as described in the following.
Our field studies suggest that the rate of I. scapularis establishment may be somewhere between the 'fast' and 'slow' scenarios, i.e. it is likely that not all CSDs identified as having 'moderate risk' for I. scapularis populations actually contain them, but some possibly do. The rate at which I. scapularis populations establish in climatically suitable areas, given a particular rate of tick immigration needs further field study to enhance our power to predict tick establishment. Increasing knowledge of the location of I. scapularis populations in Canada will also allow us to refine the relationships between climatic suitability and tick immigration rates by more formal spatial regression analyses. Furthermore, we need to investigate to what extent land mammals such as deer  may become more important than migratory birds in dispersing I. scapularis once the ticks become more widely established in south eastern Canada. Dispersion by land mammals would introduce more biologically explicit spatial autocorrelation amongst location-specific probabilities of I. scapularis population establishment.
Habitat, as described by percentage forest cover of CSDs was not a factor that improved prediction of current I. scapularis populations even though there is a priori evidence that habitat is an important factor for I. scapularis survival, and a determinant of the density of the mammalian and avian hosts of the tick . It is likely that sufficient suitable habitat exists in all CSDs even if the percent coverage is small. Indeed, investigations not described here suggest that percentage forest cover could have a non-linear relationship with tick population occurrence: a high percentage of cover could mean a lot of incoming migratory birds and a lot of habitat suitable for tick survival, but small percentages of forest cover may mean that migratory birds are concentrated into finite areas with correspondingly high densities of bird-borne ticks. However, the estimates of forest cover are crude because a number of averaging steps were used, including averaging of remote-sensed AVHRR data into land cover maps and averaging of land cover maps per CSD. Therefore the forest cover estimates are likely to be crude indices of habitat suitability for I. scapularis, its hosts and the numbers of migratory passerines such habitats may attract.
Further studies are required to expand our current knowledge of the suitability of habitat for I. scapularis in Canada. Habitat may be a more useful predictor of the occurrence of I. scapularis at a finer spatial scale than that used here . Further studies are under way in southern Quebec and elsewhere in Canada to investigate the rates at which the agent of Lyme disease, B. burgdorferi becomes established in newly-established I. scapularis populations, to determine at what rate the risk of a tick bite becomes the risk of an infected tick bite. The rate of B. burgdorferi establishment may depend on factors that are not explicitly considered in our algorithms at present, such as reservoir host densities, although in the initial expansion zone in southeastern Canada, the densities of reservoir hosts and adult tick hosts such as deer are not considered limiting factors .
Habitat may have another effect on predictions. If off-host tick mortality is higher or lower than that set in the I. scapularis population model, then the model can slightly over- or under-estimate (respectively) tick abundance [10, 19]. The model was parameterised in woodland habitats of central/eastern Canada, and greater tick survival in woodland habitats of Manitoba is a possible explanation for the survival of the I. scapularis at the Manitoba site when the model predicted population die-out. In other studies, ROC analysis of the predictions of logistic regression models has been used to identify geographic variation in the predictive power of regression models . Here our method identified possible regional variations in the predictive power of the risk maps, although we cannot, without field studies, rule out the possibility that the microclimate at the Manitoba site is in fact suitable for tick survival as the model predicts.
For simplicity we have only used one emissions scenario and the output from one global climate model in our study. More detailed studies await i) output from the third version of the Canadian Coupled Global Climate Model (CGCM3, which is almost ready for use by the scientific community), ii) new emissions scenarios from the fourth assessment of the Intergovernmental Panel on Climate Change (IPCC), and iii) projected climate at a much finer geographic scale courtesy of regional climate models for North America . If and when I. scapularis becomes established in the more northern regions of risk identified by the risk maps, then a different unit of spatial scale than CSD will be needed as CSDs, which partly reflect population densities, become vary large in northern Canada.
We have prepared risk maps for public health use now for identifying current and possible future areas of risk for I. scapularis ticks and, by inference, possible current or future Lyme disease risk. The risk maps suggest a possible widespread expansion of I. scapularis across southeastern, and eventually south central Canada during the coming decades, while our field study suggests that this process could already have begun. We have described a process for creating an index of risk for the occurrence of arthropod vectors of public health importance in Canada using information on the biology of vector survival and dispersion. As such, this could act as a model for assessing the spatial and temporal risks from other vector-borne zoonoses in Canada that are not currently endemic but which may become so due to climate change. Such modelling and mapping studies need, however, the support of solid field and laboratory studies. The difference in the rates and ranges of I. scapularis between the two scenarios for the rates of tick establishment in suitable habitat illustrates this point. More accurate prediction of where I. scapularis occurs now (i.e. the start point for simulations) and at what rate I. scapularis will establish, and are establishing in Canadian habitats, can only be elucidated using refinements of the risk maps presented here. This will require careful and more extensive field studies than those conducted here. Such studies are therefore a priority if we are to achieve our goal of limiting the impact of Lyme disease in Canada by targeted surveillance and preventive measures.
All mapping and geospatial analyses were performed in ESRI ArcGIS 9.1 (ESRI 2004). All layers were in Canada Albers Equal Area Conic projection (datum NAD83 and ellipsoid GRS80). There were 10030 complete or partial (i.e. divided by geographic features such as rivers) CSDs east of longitude 100°W used in this study. The mean area of these CSDs, an indication of the spatial resolution of the risk maps, was 285 km2 (SD = 5750 km2), although this value was much lower, with less variation amongst CSDs, for CSDs close to the USA border (mean 85 km2, SD = 277 km2 for CSDs south of latitude 46°N in southern Quebec, Ontario and the Maritime Provinces).
Risk factors for I. scapularis occurrence used in developing a risk algorithm
We have a priori information to suggest that I. scapularis distributions in Canada are dependent on combined effects of three key variable: ambient temperature, habitat and the numbers of immigrating ticks on migratory birds [10, 19, 20]. Therefore, these variables were selected as components of a risk algorithm as described in the following.
In a previous modelling study, the maximum annual number of feeding adult female ticks at model equilibrium (as an index of overall tick abundance, and thus climatic suitability for I. scapularis) was linearly related to the cumulative annual number of degree-days > 0°C (DD > 0°C ). We therefore created an interpolated surface of DD > 0°C, using Environment Canada 1971–2000 temperature normals, for all Canadian meteorological stations that have at least 15 years data in this period. Data from 2368 meteorological stations were used and the mean distance between these stations for the whole of Canada was 24.6 km (SD = 36.8 km), although in southern Canada the mean distance was much smaller (15.4 km, SD = 12.6 km). Exact local deterministic interpolators (Radial Basis Function [RBF] and Inverse Distance Weighted [IDW]: [30, 31]) were used to interpolate the values and differences between the methods investigated. The best statistical results (mean and root mean square of predictions errors) were obtained with RBF (thin-plate spline function). However, RBF predicted values above the maximum and below the minimum measured values. In particular, RBF produced positive values for the number of ticks > 1 in 'ripples' in areas beyond the northern edge of predicted regions of climate suitability for the ticks, i.e. where all surrounding recorded values had DD > 0°C value below the threshold for tick population survival. We considered that this was a problem that was important to eliminate in our study because we wished to identify areas that are predicted to be unsuitable for I. scapularis establishment as well as those that are predicted to be suitable. IDW was the method chosen because it is a basic interpolation method that predicted values inside the minimum and maximum value range and also because the data were normally distributed. We did not account for altitude in interpolations because wide variations in altitude between meteorological stations (as occur in western Canada where they must be accounted for ) are uncommon in central and eastern Canada except in Labrador and to a lesser extent the northern end of the Appalachian chain in Quebec. In these regions, incorporation of altitude would be necessary in finer scale risk maps for targeting surveillance activities.
The maximum DD > 0°C value within each CSD was assigned as the DD > 0°C value for that CSD and then converted to the expected number of ticks at model equilibrium using the relationship in reference . The maximum value was used to avoid misclassifying CSDs as wholly unsuitable for I. scapularis establishment when they could actually contain I. scapularis populations at least at some locations. We used the relationship between DD > 0°C and the number of ticks obtained in the I. scapularis population model obtained when using temperature data from locations in southern Ontario . Below a threshold minimum value of 3063 DD > 0°C, deterministic die out of tick populations occurred in model simulations. As a result, for values of DD > 0°C at or below this threshold, the number of ticks at model equilibrium was zero.
This process of interpolating DD > 0°C, conversion to the expected number of ticks, and assigning a value for each CSD was repeated using downscaled output from the second version of the Coupled Global Climate Model (CGCM2) of the Canadian Centre for Climate Modelling and Analysis [32, 33] as previously described . We obtained output from the model that incorporated estimated forcing due to emissions of greenhouse and other gases calculated in emission scenario forcing 'A2' for global economic change for the coming decades as defined by the Intergovernmental Panel on Climate Change . Briefly, in scenario 'A2', the future world is very much as it is at present in socio-economic terms: very heterogeneous with a continuously increasing population, regionally oriented economic development, and fragmented per capita economic growth and technological change . DD > 0°C were calculated at each grid point by accumulating the daily positive mean temperature values throughout the whole year. In each case, the process was repeated for all 30 years of the three future time periods, and the climatic 'normal' of each period at each grid point was obtained by averaging 30 years of annual DD > 0°C. We therefore had interpolated surfaces of estimated I. scapularis abundance (assuming temperature was the only determinant of survival, as previously discussed [1, 5, 19]) at four time periods: the 'present day' (i.e. using recorded meteorological data for the period 1971 to 2000), and three future time slices centred around the 2020s, 2050s and 2080s.
Habitat is a crucial determinant of I. scapularis survival [19, 35]. All studies associate I. scapularis with woodland or woodland ecotone habitats [11, 12, 19, 35], but clear associations with particular woodland types and other environmental indicators within a study (e.g. soil types and drainage) conflict amongst studies [e.g. [11, 12, 19]]. In general in the USA, coniferous forest is considered unsuitable habitat for I. scapularis and protective against establishment of the tick [11, 12]. However in Nova Scotia, I. scapularis populations have established in coniferous woodland [unpublished L.R. Lindsay]. While there are likely to be common habitat factors that are fundamental to tick survival amongst all sites, it is possible that these have been overlooked in the studies to date. Habitat factors identified to date as risk factors for tick occurrence are possibly proxies for the real habitat determinants of tick survival. Because of this uncertainty, we simply assigned to each CSD the percentage forest cover in that CSD (excluding water bodies) to represent the proportion of potentially suitable habitat for I. scapularis and thus a potential determinant of I. scapularis population occurrence. Forest data based on AVHRR remote-sensed proxies for habitat were obtained from a Natural Resources Canada land cover map at a resolution of 1 km2  and percentage forest cover calculated for each CSD in ArcGIS.
Indices of tick immigration
The number of ticks immigrating into a particular CSD in Canada will be proportional to the number of I. scapularis-carrying birds that arrive in spring in that CSD, and the number of ticks that each bird carries. As described previously , a wide range of passerine species migrating into Canada each spring carry I. scapularis, and these birds fly along a wide range of routes, predominantly in a general South-North direction, but with much movement across country too . Therefore, we constructed indices of tick immigration rates for each CSD as follows.
First, a map of US counties endemic for I. scapularis was created from the data in Dennis et al  using the georeferenced data file for US counties (the spatial resolution of the data in ) in ArcGIS, which was joined to a georeferenced data file for endemic Canadian CSDs. A 'buffer' with a radius of 425 km was then created for each CSD centred on the centroid of each CSD. A distance of 425 km was chosen because this is the mean distance that a passerine will fly in 5 days (the very maximum length of time that a nymphal I. scapularis would likely remain attached to a host: ), assuming a mean 85 km flown per day . The number of US counties and Canadian CSDs endemic for I. scapularis within the 425 km radius from the centroid of each CSD was calculated and this became the index of immigration for nymphal I. scapularis for each CSD. While I. scapularis nymphs must feed for 4 days, the birds could stop in a CSD for longer periods due to inclement weather, or they could stop to establish breeding territory, so 425 km represented the maximum distance the ticks could be dispersed, but not the minimum distance they could be transported. By a similar method, an index of larval tick immigration was created using a buffer with a radius of 255 km, the maximum distance a larval tick could be transported were it to remain on a bird for three days.
The index of nymphal tick immigration was validated against data from passive surveillance for I. scapularis from 1990 to 2003 in which the same CSDs used in the current study formed the spatial unit . Two methods of validation were used. First we used the whole database of locations of I. scapularis submitted by the public, via provincial organisations to the Public Health Agency of Canada, to compare how well the index of nymphal tick immigration predicted the number of ticks submitted by the public from each CSD in negative binomial regression models in STATA version 8.0 for Windows (STATACorp, College Station, TX, USA). Province of origin was included as a fixed effect because of differences amongst provinces in the effort expended in passive surveillance, and the natural logarithm of the human population in each CSD (as estimated in the 2001 census ) was an offset in the model, because the number of submissions inevitably depends on the human population density .
Second, we investigated the index of nymphal tick immigration as a predictor of the rate of submission of I. scapularis in the existing passive surveillance programme conducted in southern Quebec. Here surveillance has been most intense and long-lived, and risk factors for I. scapularis have already been identified . These factors were mean spring NDVI (as a remotely-sensed proxy for forest cover) and latitude . This earlier study raised the hypothesis that the observed geographic distribution of I. scapularis occurrence in southern Quebec was determined by the distance birds travel from endemic areas in the USA, a south-north temperature gradient (with the variable 'latitude' possibly acting as a proxy for one or both of these factors), and the suitability of habitat for bird-borne nymphal ticks to survive the moult to become questing adults. If so, then our indices of tick immigration should also partly explain the occurrence of I. scapularis in Quebec, and in particular be a confounder of 'latitude' as an explanatory variable. We therefore re-created the negative binomial regression model used in reference  in STATA (with the natural logarithm of the human population in each CSD as an offset), and tested whether the nymphal tick immigration index partly explained I. scapularis submissions, and to what extent confounding between the index of nymphal tick immigration and latitude occurred. The fit of the different models was compared using AIC .
Selection of risk algorithms and risk map construction
While we have a priori information to suggest the dependence of I. scapularis distributions in Canada on combined effects of temperature, habitat and tick dispersion on migratory birds, we have no information on the relative magnitude of these effects. The small number of known I. scapularis populations in Canada provides little scope to explore these relationships in formal spatial logistic regression analyses, and obtain accurate estimates for the relative effects of different explanatory variables. Furthermore, here we also wished to identify locations that are not suitable for tick establishment (which is important in Lyme disease diagnosis ), which could arise if i) there were no immigrating ticks, irrespective of whether temperature or habitat were suitable, or ii) temperature was not suitable or there was no suitable habitat even though there may be many immigrating ticks. It is possible that multivariable logistic regression models could identify a location as suitable even if one explanatory variable alone predicts a very low probability for the occurrence of a tick population.
We therefore opted, as a first step in risk mapping, to use simple algorithms that combined the effects of these factors by multiplication into an index of risk for I. scapularis population establishment in Canada. A range of algorithms were chosen (Table 1), which predicted i) a risk index of zero for I. scapularis establishment in a CSD if the value for any one factor for that CSD was zero, or ii) a positive value obtained by multiplication of the different components of the algorithm. More complex non-linear relationships were not investigated due to the absence of a priori hypotheses. There was little correlation amongst the three components of these algorithms (number of ticks at model equilibrium, indices of tick immigration and percent forest cover). Pearson correlation coefficients were less than 0.21 in each pair wise comparison, with the exception that the correlation coefficients for the index of larval tick immigration and the number of ticks at model equilibrium, and the index of nymphal tick immigration were 0.58 and 0.57 respectively.
We then compared how well each algorithm predicted the occurrence of the 16 CSDs or partial CSDs in Canada in which, prior to our field validation, I. scapularis populations were known or suspected (i.e. under investigation) to have established since 2000 (one in Manitoba, three in Nova Scotia and the rest in Ontario), and 10014 presumed negative CSDs. Note that two populations (one in Nova Scotia and one in Ontario) straddle more than one CSD resulting in a greater number of I. scapularis-positive CSDs than there are known populations. All the I. scapularis populations have been identified, and their geographic scope determined, by a combination of data on passive surveillance for ticks and the occurrence of human cases, and intensive field studies using established criteria . ROC analysis in STATA version 8.0 for Windows was used for these comparisons because it explicitly generates specificity and sensitivity values for a range of cut-off values  that can be used to select cut-off values for classifying CSDs for risk mapping purposes. We did, however, compare whether AUC values in ROC analysis and AIC values for logistic regression models gave similar results in terms of algorithm selection. Logistic regression models were created in STATA version 8.0 for Windows for each algorithm. In these models I. scapularis population occurrence was the outcome, and the components of the risk algorithms were the explanatory variables. We used ROC under the assumption that, for the most part known I. scapularis populations are spatially separated [1, 20], and statistically independent of one another . Predictions were made with only I. scapularis-endemic USA counties and Canadian CSDs known to contain tick populations prior to 2000 as source locations for immigrating ticks. Only data from CSDs east of 120°W (i.e. east of the Rocky Mountains) and south of 55°N were used in this analysis to limit the number of negative values that would be certainly predicted. Throughout this evaluation process, we used algorithm values in which the number of ticks was calculated from DD > 0°C values obtained from interpolations of 1971–2000 temperature normals.
The best performing algorithm (i.e. that having the highest area under the curve in ROC analysis ) was then used to develop risk maps. ROC was used for algorithm selection rather than specifically to estimate the accuracy of each algorithm. The Youden index (J = Sensitivity + Specificity - 1 ) was calculated for each of 1138 algorithm value points on the ROC curve. The algorithm value giving two thirds of the highest Youden Index was selected as the cut-off point above which CSDs were classified as likely to contain an I. scapularis population by the end of 2019 (i.e. belonging to what we termed the 'moderate risk' group of CSDs) assuming that temperature conditions do not change during the period 2000 to 2019. This cut-off value gave a sensitivity of 94% and specificity of 66% for detecting a CSD as likely to contain an I. scapularis population according to the ROC analysis. The 10% of the CSDs in the 'moderate risk' group with the highest algorithm values were re-classified as the 'high risk' group of CSDs for I. scapularis populations. The cut-off value for the 'high risk' group had a sensitivity of > 99% and specificity of < 1% for I. scapularis population occurrence according to ROC analysis. Three further groups were: i) low risk CSDs, i.e. those with an algorithm value > 1, but below the cut-off for inclusion in the 'moderate risk' group; ii) CSDs where there is no risk of I. scapularis population establishment (i.e. the algorithm value is zero) but those where there is a risk from bird-borne adventitious ticks (i.e. the value for the index of tick immigration > 1); and iii) a 'no risk' group of CSDs where both the algorithm value and the index of tick immigration were zero (i.e. CSDs outside the 425 km range of tick dispersion by migratory birds).
Mapping risk under climate change scenarios
The selected algorithm could operate as a simple simulation model with 3 time steps. The starting point was a map of CSDs with a high index for I. scapularis population occurrence under current temperature conditions, which was combined with the map of US counties endemic for the tick . This information provided a baseline map of the possible geographic range for I. scapularis prior to climate change.
The next time step is the 2020s, for which the index of nymphal immigration was recalculated but this time using a larger radius of 458 km because bird migration speed in North America is expected to increase with higher spring temperatures . While the onset of migration of most migratory passerine species that nest in Canada may be controlled by endogenous circannual rhythms and photoperiod, the speed at which birds migrate north increases with increasing mean spring temperatures, possibly due to effects of temperature on the availability of the food the birds need for refuelling during their migration . The change in speed (from 85 to 91.6 km per day) was estimated using the relationship observed between spring (April/May) temperatures in southern Ontario and Pennsylvania, and the speed with which migratory birds complete the journey from Louisiana to southern Ontario or Pennsylvania (for a 1°C increase in temperature, birds arrive one day earlier ). The temperature rise used to predict changes in migratory bird speed was the difference in mean spring (April/May) temperature for the Chatham, Ontario meteorological station during 1971–2000 and that temperature predicted for this station for the 2020s by downscaling and interpolation of CGCM2 output under IPCC emissions scenario A2. The maximum DD > 0°C (and thus the maximum number of ticks at model equilibrium) for each CSD also changed according to the interpolated, downscaled output from CGCM2 for the 2020s as described in the previous section on temperature. CSDs were then allocated to one of the five classes (high risk, moderate risk, low risk, risk from adventitious ticks and no risk) according to the new algorithm value obtained.
In turn, these populations, combined with pre-existing positive Canadian CSDs and US counties in the baseline, pre-climate change map, provided the 'source' of ticks for I. scapularis to establish in CSDs (that were previously free of the tick) by the 2050s. Again, the index of nymphal immigration was recalculated with the radius increased from 458 to 486 km, and the DD > 0°C (and thus the maximum number of ticks at model equilibrium) for each CSD changed according to the interpolated, downscaled output from CGCM2 for the 2050s. Each CSD was then re-classified in terms of risk of I. scapularis population establishment, and the process was repeated again for the 2080s with increasing radius of bird movement (from 486 to 532 km), DD > 0°C and modeled tick numbers using CGCM2 output for the 2080s.
At present, the true geographic extent of reproducing I. scapularis populations in Canada is uncertain, at least in part because the number of newly established populations is increasing and active surveillance for these populations to date has limited scope and is incomplete. For this reason, we ran the 'model' using two different scenarios: i) the 'fast' scenario in which a CSD was considered to contain an I. scapularis population (and would thus act as 'source' for the ticks) if the CSD were in the 'moderate risk' or 'high risk' groups, and ii) the 'slow' scenario in which a CSD was considered to contain an I. scapularis population only if the CSD were in the 'high risk' group.
Field validation of risk maps under current climate conditions
Validation in the field was carried out in southern Quebec. We aimed to visit 70 sites (30 sites in 'moderate risk' CSDs and 40 sites in 'high risk' CSDs) giving 90% power to identify a difference in the possible prevalence of I. scapularis-infested sites predicted by the risk algorithm of 66% and 99% in moderate risk and high risk regions respectively with α = 0.05. However, for a number of reasons including flooding, clear cutting of trees, lack of permission and a total absence of trappable rodents, data were collected from 45 sites which gave a power of 75% to detect such a difference in prevalence. Teams of two researchers visited each site once from 5th June to 5th October 2007. Sites were selected to fulfill the following criteria: one site per CSD in two regions (Montérégie and Estrie) determined by the practicality of traveling distances, and deciduous woodlands (determined from provincial government forest cover maps ) of minimal dimensions 500 m by 150 m. Sites were visited in a haphazard fashion depending on contact with, and permission to study from, the landowners.
At each visit, rodents were trapped and host-seeking ticks were collected according to the following protocol. At each site 150 Sherman traps were placed in three parallel transects of 50 traps. The transects were placed 25–50 m apart, and within a transect, the traps were set 10 m apart. Traps were furnished with bedding and a mix of nuts, seeds and peanut butter as bait. Traps were set for two nights or one night if 15 or more rodents were captured in the first night (with α = 0.05, a sample of 15 rodents gives 80% power to identify an infested rodent if the true prevalence of infestation is 5%). All captured rodents were examined for the presence of ticks after light anaesthesia as described elsewhere . The herbage of each site was flagged to collect host-seeking ticks for 3 person-hours in the afternoon. All ticks collected from rodents and herbage were identified to species by standard keys [42, 43].
An 'index of certainty' that a reproducing population of I. scapularis ticks was present was determined for each site. Increasing values for this index gave increased confidence that ticks found at a site were from a reproducing population of ticks, rather than just bird-transported 'adventitious' ticks [20, 23]. The index was calculated for each site as follows: presence of larvae = 1 point; presence of 3 or more larvae = 1 additional point; presence of 10 or more larvae = 1 additional point; presence of nymphs = 1 point; presence of 3 or more nymphs = 1 additional point; presence of 10 or more nymphs = 1 additional point; presence of adult ticks = 1 point; presence of 3 or more adults = 1 additional point; presence of 10 or more adults = 1 additional point. Thus the index had a minimum value of 0 (no I. scapularis found at the site) and a maximum of 9 (when 10 or more ticks of each of the three instars were found at the site). All ticks found on flagging and on rodents were included in this calculation.
The index of certainty was then the outcome in ordered logistic regression models which make no a priori assumptions as to the linearity or otherwise of the index, simply assuming that 1 is 'greater' than 0, 2 'greater' than 1 etc. . Site classification (i.e. high risk or moderate risk), the value of the algorithm for the CSD in which the site occurred, and, independently, the tick immigration index and number of ticks at model equilibrium on which the algorithm for that CSD was based, were explored as predictors of the occurrence of I. scapularis. The fit of the different models in predicting the occurrence of likely I. scapularis populations was investigated using AIC. Potential confounding explanatory variables were explored including the month the site was visited (because different I. scapularis instars are active in different seasons ), the ID for the field team and the number of rodents captured at the site.
This study was funded by the Climate Change Impacts and Adaptation Programme of Natural Resources Canada. Field validation was jointly funded by Institut national de santé publique du Québec, and the Public Health Agency of Canada.
- Ogden NH, Lindsay LR, Sockett PN, Morshed M, Artsob H: The rising challenge of Lyme disease in Canada. Can Commun Dis Rep 2008, 34:1–19.PubMed
- Lindsay R, Artsob H, Barker I: Distribution of Ixodes pacificus and Ixodes scapularis re concurrent babesiosis and Lyme disease. Can Commun Dis Rep 1998, 24:121–122.PubMed
- Parmesan C, Yohe G: A globally coherent fingerprint of climate change impacts across natural systems. Nature 2003, 421:37–42.View ArticlePubMed
- Root TL, Price JT, Hall KR, Schneider SH, Rosenzweig C, Pounds JA: Fingerprints of global warming on wild animals and plants. Nature 2003, 421:57–60.View ArticlePubMed
- Ogden NH, Maarouf A, Barker IK, Bigras-Poulin M, Lindsay LR, Morshed MG, O'Callaghan CJ, Ramay F, Waltner-Toews D, Charron DF: Projections for range expansion of the Lyme disease vector Ixodes scapularis , in response to climate change. Int J Parasitol 2006, 36:63–70.View ArticlePubMed
- Kurtenbach K, Hanincova K, Tsao J, Margos G, Fish D, Ogden NH: Key processes in the evolutionary ecology of Lyme borreliosis. Nature Rev Microbiol 2006, 4:660–669.View Article
- Wormser GP, Dattwyler RJ, Shapiro ED, Halperin JJ, Steere AC, Klempner MS, Krause PJ, Bakken JS, Strle F, Stanek G, Bockenstedt L, Fish D, Dumler JS, Nadelman RB: The clinical assessment, treatment, and prevention of Lyme disease, human granulocytic anaplasmosis, and babesiosis: clinical practice guidelines by the Infectious Diseases Society of America. Clin Infect Dis 2006, 43:1089–1134.View ArticlePubMed
- Aguero-Rosenfeld ME, Wang G, Schwartz I, Wormser GP: Diagnosis of Lyme borreliosis. Clin Microbiol Rev 2005, 18:484–509.View ArticlePubMed
- Brownstein JS, Holford TR, Fish D: A climate-based model predicts the spatial distribution of the Lyme disease vector Ixodes scapularis in the United States. Environ Health Perspect 2003, 111:1152–1157.View ArticlePubMed
- Ogden NH, Bigras-Poulin M, O'Callaghan CJ, Barker IK, Lindsay LR, Maarouf A, Smoyer-Tomic KE, Waltner-Toews D, Charron D: A dynamic population model to investigate effects of climate on geographic range and seasonality of the tick Ixodes scapularis.Int J Parasitol 2005, 35:375–389.View ArticlePubMed
- Dister SW, Fish D, Bros SM, Frank DH, Wood BL: Landscape characterization of peridomestic risk for Lyme disease using satellite imagery. Am J Trop Med Hyg 1997, 57:687–692.PubMed
- Guerra M, Walker E, Jones C, Paskewitz S, Cortinas MR, Stancil A, Beck L, Bobo M, Kitron U: Predicting the risk of Lyme disease: habitat suitability for Ixodes scapularis in the north central United States. Emerg Infect Dis 2002, 8:289–297.View ArticlePubMed
- Benedict MQ, Levine RS, Hawley WA, Lounibos LP: Spread of the tiger: global risk of invasion by the mosquito Aedes albopictus.Vector Borne Zoonotic Dis 2007, 7:76–85.View ArticlePubMed
- Rakotomanana F, Randremanana RV, Rabarijaona LP, Duchemin JB, Ratovonjato J, Ariey F, Rudant JP, Jeanne I: Determining areas that require indoor insecticide spraying using Multi Criteria Evaluation, a decision-support tool for malaria vector control programmes in the Central Highlands of Madagascar. Int J Health Geogr 2007, 6:2.View ArticlePubMed
- Tran A, Poncon N, Toty C, Linard C, Guis H, Ferre JB, Lo Seen D, Roger F, de La Rocque S, Fontenille D, Baldet T: Using remote sensing to map larval and adult populations of Anopheles hyrcanus (Diptera: Culicidae) a potential malaria vector in Southern France. Int J Health Geogr 2008, in press
- Estrada-Peña A, Zatansever Z, Gargili A, Aktas M, Uzun R, Ergonul O, Jongejan F: Modeling the spatial distribution of Crimean-Congo hemorrhagic fever outbreaks in Turkey. Vector Borne Zoonotic Dis 2007, 7:667–678.View ArticlePubMed
- Brownstein JS, Rosen H, Purdy D, Miller JR, Merlino M, Mostashari F, Fish D: Spatial analysis of West Nile virus: rapid risk assessment of an introduced vector-borne zoonosis. Vector Borne Zoonotic Dis 2002, 2:157–164.View ArticlePubMed
- Tachiiri K, Klinkenberg B, Mak S, Kazmi J: Predicting outbreaks: a spatial risk assessment of West Nile virus in British Columbia. Int J Health Geogr 2006, 5:21.View ArticlePubMed
- Ogden NH, Barker IK, Beauchamp G, Brazeau S, Charron D, Maarouf A, Morshed MG, O'Callaghan CJ, Thompson RA, Waltner-Toews D, Waltner-Toews M, Lindsay LR: Investigation of ground level and remote-sensed data for habitat classification and prediction of survival of Ixodes scapularis ticks in habitats of southeastern Canada. J Med Entomol 2006, 43:403–414.View ArticlePubMed
- Ogden NH, Lindsay LR, Hanincová K, Barker IK, Bigras-Poulin M, Charron DF, Heagy A, Francis CM, O'Callaghan CJ, Schwartz I, Thompson RA: The role of migratory birds in introduction and range expansion of Ixodes scapularis ticks, and Borrelia burgdorferi and Anaplasma phagocytophilum in Canada. Appl Environ Microbiol 2008, 74:1780–1790.View ArticlePubMed
- Marra PP, Francis CM, Mulvihill RS, Moore FR: The influence of climate on the timing and rate of spring bird migration. Oecologia 2005, 142:307–315.View ArticlePubMed
- Thompson C, Spielman A, Krause PJ: Coinfecting deer-associated zoonoses: Lyme disease, babesiosis, and ehrlichiosis. Clin Infect Dis 2001, 33:676–685.View ArticlePubMed
- Ogden NH, Trudel L, Artsob H, Barker IK, Beauchamp G, Charron D, Drebot MA, Galloway TD, O'Handley R, Thompson RA, Lindsay LR:Ixodes scapularis ticks collected by passive surveillance in Canada: analysis of geographic distribution and infection with the Lyme borreliosis agent Borrelia burgdorferi.J Med Entomol 2006, 43:600–609.View ArticlePubMed
- Statistics Canada [http://www12.statcan.ca/english/census01/home/index.cfm]
- Greiner M, Pfeiffer D, Smith RD: Principles and practical application of the receiver-operating characteristic analysis for diagnostic tests. Prev Vet Med 2000, 45:23–41.View ArticlePubMed
- Dennis DT, Nekomoto TS, Victor JC, Paul WS, Piesman J: Reported distribution of Ixodes scapularis and Ixodes pacificus (Acari: Ixodidae) in the United States. J Med Entomol 1998, 35:629–638.PubMed
- Madhav NK, Brownstein JS, Tsao JI, Fish D: A dispersal model for the range expansion of blacklegged tick (Acari: Ixodidae). J Med Entomol 2004, 41:842–852.View ArticlePubMed
- Brooker S, Hay SI, Issae W, Hall A, Kihamia CM, Lwambo NJ, Wint W, Rogers DJ, Bundy DA: Predicting the distribution of urinary schistosomiasis in Tanzania using satellite sensor data. Trop Med Int Health 2001, 6:998–1007.View ArticlePubMed
- Ouranos [http://www.ouranos.ca/doc/ACFAS2006/Presentations/Frigon.pdf]
- Shepard D: A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 1968 ACM National Conference517–524.
- Duchon J: Splines minimizing rotation invariant seminorms in Sobolev spaces. Topics in Constructive Theory of Functions of Several Variables Berlin Heidelberg New York: Springer Verlag 1976, 1:85–100.
- Flato GM, Boer GJ, Lee WG, McFarlane NA, Ramsden D, Reader MC, Weaver AJ: The Canadian Centre for Climate Modelling and Analysis Global Coupled Model and its Climate. Clim Dynam 2000, 16:451–467.View Article
- Flato GM, Boer GJ: Warming Asymmetry in Climate Change Simulations. Geophys Res Lett 2001, 28:195–198.View Article
- Intergovernmental Panel on Climate Change (IPCC): Special Report on Emission Scenarios (SRES) Cambridge: Cambridge Univ. Press 2000.
- Lindsay LR, Barker IK, Surgeoner GA, McEwen SA, Gillespie TJ, Addison EM: Survival and development of the different life stages of Ixodes scapularis (Acari: Ixodidae) held within four habitats on Long Point, Ontario, Canada. J Med Entomol 1998, 35:189–199.PubMed
- Natural Resources Canada [http://ccrs.nrcan.gc.ca/optic/map_e.php]
- The Canadian Atlas of Bird-Banding [http://www.cws-scf.ec.gc.ca/publications/BBA-AOB/v1ed2/index_e.cfm]
- Akaike H: A new look at the statistical model identification. IEEE TAC 1974, AC-19:716–723.
- Brito JC, Crespo EG, Paulo OS: Modelling wildlife distributions: Logistic regression versus overlap analysis. Ecography 1999, 22:251–260.View Article
- Ministère des Ressources naturelles et de la Faune [http://www.mrnf.gouv.qc.ca/forets/connaissances/connaissances-inventaire-cartes-sief-peuplements.jsp]
- Lindsay LR, Mathison SW, Barker IK, McEwen SA, Surgeoner GA: Abundance of Ixodes scapularis (Acari: Ixodidae) larvae and nymphs in relation to host density and habitat on Long Point, Ontario. J Med Entomol 1999, 36:243–254.PubMed
- Clifford CM, Anastos G, Elbl A: The larval ixodid ticks of the eastern United States. Misc Pub Ent Soc America 1961, 2:215–244.
- Keirans JE, Hutcheson HJ, Durden LA, Klompen JSH:Ixodes ( Ixodes ) scapularis (Acari: Ixodidae): Redescription of all active stages, distribution, hosts, geographical variation, and medical and veterinary importance. J Med Entomol 1996, 33:297–318.PubMed
- Long JS, Freese J: Regression models for categorical dependent variables using Stata. College Station, Texas: Stata Press 2001.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.