ResearchOpen Access

Expansion of the Lyme Disease Vector Ixodes Scapularis in Canada Inferred from CMIP5 Climate Projections

Published:CID: 057008 by:7



A number of studies have assessed possible climate change impacts on the Lyme disease vector, Ixodes scapularis. However, most have used surface air temperature from only one climate model simulation and/or one emission scenario, representing only one possible climate future.


We quantified effects of different Representative Concentration Pathway (RCP) and climate model outputs on the projected future changes in the basic reproduction number (R0) of I. scapularis to explore uncertainties in future R0 estimates.


We used surface air temperature generated by a complete set of General Circulation Models from the Coupled Model Intercomparison Project Phase 5 (CMIP5) to hindcast historical (1971–2000), and to forecast future effects of climate change on the R0 of I. scapularis for the periods 2011–2040 and 2041–2070.


Increases in the multimodel mean R0 values estimated for both future periods, relative to 1971–2000, were statistically significant under all RCP scenarios for all of Nova Scotia, areas of New Brunswick and Quebec, Ontario south of 47°N, and Manitoba south of 52°N. When comparing RCP scenarios, only the estimated R0 mean values between RCP6.0 and RCP8.5 showed statistically significant differences for any future time period.


Our results highlight the potential for climate change to have an effect on future Lyme disease risk in Canada even if the Paris Agreement’s goal to keep global warming below 2°C is achieved, although mitigation reducing emissions from RCP8.5 levels to those of RCP6.0 or less would be expected to slow tick invasion after the 2030s.


A cause for concern to the Canadian population and public health sector is the emergence and establishment or re-establishment of climate-sensitive infectious diseases, particularly Lyme disease (Bush et al. 2014). Lyme disease is caused by the bacterium Borrelia burgdorferi and is transmitted by ixodid ticks. The main vectors are Ixodes pacifus and Ixodes scapularis in western and central/eastern North America, respectively (Eisen et al. 2016). Since the first reported human cases of Lyme disease in the early 1990s, Canada has experienced an increase in both the number of endemic regions for the disease and the number of annual diagnosed human cases (Ogden et al. 2008; PHAC 2015). Not only are these trends placing increased health stress on the Canadian population, but they are also likely resulting in added economic consequences for patients, communities, the healthcare system, and governments (Hinckley et al. 2014; Adrion et al. 2015; Zhang et al. 2006).

Tick survival depends on several factors, including climate suitability, host density, soil composition, and microhabitat characteristics (Leighton et al. 2012). Ogden et al. (2005, 2006, 2008) raised the hypothesis that climate suitability, predominantly temperature, is a key determinant of the survival of I. scapularis populations and consequently that a warming climate may be a key driver of the emergence of the tick in Canada. There are multiple ways in which climate may affect the survival of tick populations. First, climate factors such as temperature and humidity can affect host-seeking activity. Second, the development rate of ticks from one life stage to the next is highly dependent on temperature: the development rate is 0 at 0°C and becomes faster at higher temperatures (Ogden et al. 2014b).

To simulate the effect of temperature on tick population survival and seasonality, Ogden et al. (2005) developed a mechanistic I. scapularis population model comprised of 12 mutually exclusive states, each one representing a specific state in the ticks’ life cycle. The model was designed to simulate I. scapularis populations in specific geographic regions. To accomplish this, temperature data were used to vary the duration that ticks remained resident in the egg-laying and engorged-tick state and to alter questing-tick activities. Several parameters, including interstadial development rates, host-finding rates, and mortality rates, were calibrated for each development stage of the population model. Empirical validation revealed that the model was able to closely replicate the seasonal activity patterns for all three tick instars that were observed in the field. Using this model with temperature input obtained from two climate models, the effect of projected future climate on the possible geographic range of I. scapularis populations was explored under two different emissions scenarios. All the research was in agreement that northern expansion of the geographic range suitable for I. scapularis establishment may take place under future climate conditions (Ogden et al. 2006).

Wu et al. (2013) modified the mechanistic I. scapularis model developed by Ogden et al. (2005). This revised population model allowed prediction of the basic reproductive number (R0) of I. scapularis for specific geographic locations in North America. The authors used the same model framework and many data values from Ogden et al. (2005). However, this version of the model used rates of development in contrast with development duration, allowing a periodic system of delay differential equations to be developed rather than a system of delay difference equations. In turn, this allowed Wu et al. (2013) to employ the next-generation matrix mathematical technique to obtain R0 values. R0 describes the propensity for a parasite or microparasite to survive and be propagated. Conditions that allow R0>1 are those under which the tick population can persist, while conditions under which R0<1 result in die-out of the tick populations or in their failure to become established (Wu et al. 2013).

Ogden et al. (2014b) used the population model developed by Wu et al. (2013) to estimate the effect of climate change on the R0 of I. scapularis. Ogden et al. (2014b) produced maps of estimated R0 values under current and projected future climate in North America; however, these maps were obtained using surface air temperature from only one climate model simulation and one emission scenario, thus representing only one possible future. The climate model chosen for mapping in Ogden et al. (2014b) was the Canadian Regional Climate Model version 4.2.3 (CRCM4.2.3) which was constrained by both the initial and the boundary conditions of the Canadian General Circulation Model (GCM) CGCM3.1,T47 and was forced by the A2 scenario of the Intergovernmental Panel on Climate Change Special Report on Emissions Scenarios (SRES). This model was selected because it represented the annual mean changes in R0 projected by an ensemble of two GCM and three Regional Climate Model (RCM) simulations (see Supplemental Material of Ogden et al. [2014b]).

Here we explore in more detail the implications of uncertainty in future climate (associated with the variability in outputs of different models and emission scenarios) on the predicted future geographic distributions of I. scapularis. By applying the methodology detailed in Ogden et al. (2014b), we hindcast historical and forecast future effects of climate change on the R0 of I. scapularis using a complete set of GCMs from the Coupled Model Intercomparison Project Phase 5 (CMIP5) (Taylor et al. 2012). This allowed us a) to evaluate the effect of multiple Representative Concentration Pathways (RCPs) on the evolution of the R0 of I. scapularis, b) to quantify and communicate the implications of climate model variability on the projected R0 values, and c) to provide R0 estimates employing the most up-to-date climate models and emission scenarios. Our study includes results from 13 GCMs for the historical period and four RCP scenarios (65 simulations). These results detail how the full range of plausible climate futures may impact Lyme disease risk emergence and provide information on which to base adaptation strategies (Mearns 2010; Preston et al. 2013). The study also allowed us to assess how mitigation efforts may impact Lyme disease risk in Canada.


General Circulation Model (GCM) Descriptions

We used a set of 13 coordinated climate models from the CMIP5 (Table 1). All models included in the CMIP5 framework must complete a set of benchmark experiments. This allows the results from each model simulation to be compared with both observations and other model simulated projections (IPCC 2013). With these available experiments, the CMIP5 allows researchers to evaluate model reliability in simulating the recent past, to provide projections of future climate change, and to gain a greater understanding of the possible causes for intermodel differences (Taylor et al. 2012).

Table 1 List of the models used in this analysis, the resolutions, and the model references.

Table 1 lists models used in the first column and its resolutions and references in the other two columns.
ModelResolution (lat°×lon°)Reference
CCSM4.00.9×1.25(100km×89km)Gent et al. (2011)
CESM1-CAM50.9×1.25(100km×89km)Meehl et al. (2013)
MIROC-ESM2.81×2.81(312km×201km)Watanabe et al. (2011)
MIROC-ESM-CHEM2.81×2.81(312km×201km)Watanabe et al. (2011)
MIROC51.4×1.4(156km×100km)Watanabe et al. (2010)
GFDL-CM32×2.5(222km×179km)Donner et al. (2011)
GFDL-ESM2G2×2.5(222km×179km)Dunne et al. (2012)
GFDL-ESM2M1.13×1.13(126km×81km)Dunne et al. (2012)
MRI-CGCM31.13×1.13(126km×81km)Yukimoto et al. (2012)
IPSL-CM5A-LR1.88×3.75(209km×268km)Dufresne et al. (2013)
IPSL-CM5A-MR1.25×2.5(139km×179km)Dufresne et al. (2013)
BCC-CSM1.12.81×2.81(312km×201km)Wu et al. (2014)
BCC-CSM1.1-M1.13×1.13(126km×81km)Wu et al. (2014)

The CMIP5 simulations include: present-day simulations, coupled-control simulations with alternate forcings, historical simulations (1850–2005), and projection simulations forced with different emission concentrations specified by several RCPs. Here we used the historical (1850–2005) and RCP simulations (2006–2100) with a daily temporal resolution.

Historical simulations cover the mid-1850s to at least 2005 (Mieville et al. 2010; Schultz et al. 2008; van der Werf et al. 2006). These simulations permit the evaluation of model performance against present and observed climate change. Additionally, historical simulations allow for the analysis of human impacts on past climate and also provide initial conditions for future scenarios of greenhouse gas (GHG) emissions and anthropogenic activities described by the RCP scenarios (Taylor et al. 2012).

The RCP climate simulations presented in the CMIP5 framework are derived from a distinct combination of economic, technological, land use, energy, demographic, and policy futures, which result from different levels of societal action that aim to prevent further changes in climate. Each of the four RCPs (RCP2.6, RCP4.5, RCP6.0, and RCP8.5) reflects a specific GHG emission trajectory and a corresponding radiative forcing starting in 2006 and culminating in 2100 (Wayne 2013). The radiative forcing levels for the four RCPs range from 2.6W/m2 to 8.5W/m2 (Masui et al. 2011; Riahi et al. 2011; Thomson et al. 2011; Van Vuuren et al. 2011). Each scenario is associated with different magnitudes of global surface temperature increase. For example, within the IPCC Fifth Assessment Report (IPCC 2013), the global mean surface temperature increases projected for 2081 to 2100, relative to 1986 to 2005, are 1.0°C (RCP2.6), 1.8°C (RCP4.5), 2.2°C (RCP6.0), and 3.7°C (RCP8.5). Historical and RCP climate simulations were used in this study to generate baseline estimates of R0 and to assess the effect of climate change on R0 of I. scapularis across a range of climate models and a range of RCP scenarios.

In this study, we used the temporal evolution of surface air temperature at daily resolution as generated by each climate model in Table 1. Because initial conditions of each GCM are sampled from an arbitrary point in a control simulation, GCM simulated temperature outputs can contain biases. To improve our projections, we used the North American Regional Reanalysis (NARR) to remove the potential biases that result from the different initial conditions employed in the simulations of each climate model. The NARR is a long-term, dynamically consistent, high-resolution, high-frequency, atmospheric and land surface hydrology data set for the North American domain (Mesinger et al. 2006). The NARR data set was generated by coupling the National Center for Environmental Prediction Eta Model and the Regional Data Assimilation System. This ensured that a range of observational data (temperature, precipitation, etc.) was incorporated in the reanalysis. At the lowest latitude, the grid resolution is roughly 0.3 degrees, or 32km×32km. All variables are provided on a Northern Hemisphere Lambert Conformal Conic grid for the following temporal frequencies: 8 times daily, once daily, and monthly (Mesinger et al. 2006).

The NARR has already been used as reference for other reanalysis projects including the California Reanalysis Downscaling at 10 km (CaRD10) (Kanamaru and Kanamitsu 2007). Other studies also support the use of the NARR for the analysis of climate processes (Fall et al. 2010). Following CMIP-WGCM-WGSIP (2011), we used the NARR surface air temperature to estimate the bias of each GCM’s surface air temperature data at each grid point. We determined the monthly biases for the period 1979 to 2005, which were then removed from the daily temperature data for each model simulation at each grid point over our spatial domain. We have followed the common practice of the modeling community: using the grid size of the most coarse-resolution model in the ensemble to interpolate the rest of the GCMs (Dunn-Sigouin and Son 2013; Kopparla et al. 2013; Willison et al. 2013). Thus, we have used the MIROC-ESM grid with 2.81°×2.81° resolution (Watanabe et al. 2011) to interpolate the CMIP5 GCMs and the NARR.

Estimating Mean R0 Values of I. scapularis

In order to estimate the temporal evolution and the spatial distribution of the basic reproductive number of I. scapularis across North America, we used Equation 1 (Ogden et al. 2014b):

R0=1.072×106DD>0°C24.65×103DD>0°C+5.556 [1]

Equation 1 describes the relationship between R0 and annual cumulative degree-days >0°C (DD>0°C). DD>0°C is defined as the annual sum of the daily temperatures in degrees Celsius above the freezing point of water. The DD>0°C were computed as the accumulation of daily temperatures >0°C for each year using the bias-corrected daily surface air temperatures at each grid point and for each model simulation. From the gridded DD>0°C data, we calculated the spatial and temporal evolution variability of R0 for each model across North America. 30-yearmultimodel mean R0 values were mapped for the time periods 1971 to 2000, 2011 to 2040, and 2041 to 2070 to facilitate comparison with the results of Ogden et al. (2014b).

Because of the quadratic nature of Equation 1, R0 intersects the degree-day domain axis at two points. This results in two threshold temperature conditions for tick persistence, one at 2843DD>0°C and the other at 1795DD>0°C. We show the spatial distribution only for the R0>1 values that were reflective of temperature conditions equal to or above 2843DD>0°C. Lyme disease risk occurs west of the Rocky Mountains, but here B. burgdorferi is transmitted by a different tick species, Ixodes pacificus (Ogden et al. 2009). Thus, we masked out all geographic locations west of the Rocky Mountains.

We also excluded from our analysis any area with an elevation greater than 500 meters based on the criteria of Diuk-Wasser et al. (2012). This work suggests that the presence of infected nymphs is closely associated with lower elevation, low vapor pressures deficit, and low seasonal temperature extremes. No infected nymphs were collected above 510 meters. This upper-elevation limit may well be due to temperature conditions being too cold at higher altitudes, although other habitat characteristics may also be unsuitable for ticks at higher altitudes. For this reason, and because the resolution of climate model output was too coarse to identify colder temperatures of the relatively few locations over 500 m in southeastern Canada, these areas were masked out.

The relationship between R0 and DD>0°C is sensitive to changes in factors such as host densities, suitable host ranges, and host-finding success (Wu et al. 2013). The estimated temperature threshold for tick establishment is largely affected by off-host tick mortality rates and by variations in temperature during the warmest months of the year (Wu et al. 2013). By applying this function to all regions pertinent to this study, we are assuming that these important factors for I. scapularis establishment are constant across both space and time.

The Kolmogorov-Smirnov test at the 95% confidence interval (Wang et al. 2003) was used to determine the statistical significance of the R0 multimodel mean increases that were estimated for each of the future time periods, 2011 to 2040 and 2041 to 2070, relative to the historical period 1971 to 2000. This was conducted for each RCP simulation. Also, statistically significant differences in the estimated R0 values, obtained from successive RCP scenarios (RCP2.6 versus RCP4.5; RCP4.5 versus RCP6.0 and RCP6.0 versus RCP8.5), were evaluated using the Kolmogorov-Smirnov test at the 95% confidence interval for the periods 2011 to 2040 and 2041 to 2070.

Additionally, we focused on three regions (northern Ontario, southern Ontario, and Nova Scotia) to illustrate the evolution of the multimodel R0 mean values, the progression of climate model uncertainty, and the variations among RCPs over the period 1860 to 2099. Ontario extends from 42°N to 57°N latitude and 75°W to 95°W longitude (Baldwin et al. 2001). Over this large area, climate exhibits large spatial variability and thus a wide range of R0 values across the province. Ontario’s 10 major cities exist below 50°N, making southern Ontario home to more than 85% of Ontario’s total population (Government of Ontario 2015). I. scapularis populations have emerged in southern Ontario and Nova Scotia but are thought to be absent from the region of northern Ontario studied here.


R0 Values of I. scapularis in North America

For the historical period of 1971 to 2000, the multimodel mean yielded R0 values ≥ 3 limited to areas in Ontario below 45°N, which is slightly higher than, but close to, the R0 values obtained using observed climate data (Ogden et al. 2014b). These lower-latitude areas of Ontario were identified as having the highest R0 values in eastern and central Canada (Figure 2a). For the historical time period, the estimated R0 values from all GCMs were in good agreement, displaying a 95% confidence interval of <1 (Figure 2d). Estimates of the geographic range where R0>1 are consistent with the observed geographic distribution of blacklegged ticks for the period 2003 to 2013 (Bouchard et al. 2015; Ogden et al. 2014a) (Figure 1). The surveillance data from Bouchard et al. (2015) and Ogden et al. (2014a) showed that few tick populations were found in the warmest parts of Ontario, which may be attributed to the relative lack of woodland habitat (Clow et al. 2016). Additionally, the predicted evolution of R0 estimates for southern Ontario and Nova Scotia were compared with the evolution in the number of Canadian Census Subdivisions (CSD) with probable established I. scapularis populations over the years between 1991 and 2008, estimated from surveillance data by Leighton et al. (2012). In the two selected regions, tick populations started to emerge, after a period of sustained warming, increasing R0 toward or over 2 for several consecutive years (2000 to 2005 in Ontario and 1998 to 2003 in Nova Scotia) (Figures S1 and S2).

Grid map of Northings (y-axis) and Eastings (x-axis) showing the regions where at least one tick and where no tick was found according to the surveillance conducted 2003–2012. The second grid map of Northings (y-axis) and Eastings (x-axis) plots the R0 estimates for the period 1971–2000.

Figure 1. Top: Data from active tick surveillance conducted between both 2003–2012 (Bouchard et al. 2015) and 2008–2013 (Ogden et al. 2014a). Locations where at least one I. scapularis tick has been found are marked by red triangles. White triangles mark locations where no I. scapularis ticks were collected. Bottom: Historical multimodel mean R0 estimates for the period 1971–2000. The R0 bar provides a color gradient showing the colors, which correspond to specific R0 values. The numbers around both the top and bottom maps depict the northings and eastings associated with the domain.

Consistent with the findings of Ogden et al. (2014b), future projected changes in climate (IPCC 2013) were suggested to have an effect on the DD>0°C and consequently on R0 of I. scapularis; increases in the multimodel mean R0 values estimated for the period 2011 to 2040 relative to 1971 to 2000 were statistically significant under all four RCP scenarios (Figures 2b and 3b and Figures S3b and S4b). Changes in R0 were also statistically significant under all four RCP emission scenarios between the historical period 1971 to 2000 and the period 2041 to 2070 (Figures 2c and 3c and Figures S4c and S5c).

Six grid maps plotting Northings (y-axis) and Eastings (x-axis). (a), (b) and (c) plot the R0 values for 1971–2000, RCP45 2011–2040, and RCP45 2041–2070, respectively. (d), (e) and (f) plot the 95% confidence interval for 1971–2000, RCP45 2011–2040, and RCP45 2041–2070, respectively.

Figure 2. Multimodel mean R0 values for I. scapularis (maps a−c) and the climate model variability (95% confidence interval) (maps d−f) for the period 1971–2000 using the historical simulation and for the periods 2011–2040 and 2041–2070 using RCP4.5 simulations. R0>1 values were mapped for areas east of the Rocky Mountains and for elevations below 500 m. The R0 bar provides a color gradient showing the colors, which correspond to specific R0 values. The numbers around the maps depict the Northings and Eastings associated with the domain. Black stippling in maps b and c shows the spatial distribution of statistically significant changes in R0 between each future time period (2011–2040 and 2041–2070) and the historical time period (1971–2000). Statistical significance was calculated using the Kolmogorov-Smirnov test.

Six grid maps plotting Northings (y-axis) and Eastings (x-axis). (a), (b) and (c) plot the R0 values for 1971–2000, RCP85 2011–2040, and RCP85 2041–2070, respectively. (d), (e) and (f) plot the 95% confidence interval for 1971–2000, RCP85 2011–2040, and RCP85 2041–2070, respectively.

Figure 3. Multimodel mean R0 values for I. scapularis (maps a−c) and the climate model variability (95% confidence interval) (maps d−f) for the period 1971–2000 using the historical simulations and for the periods 2011–2040 and 2041–2070 using RCP8.5 simulations. R0>1 values were mapped for areas east of the Rocky Mountains and for elevations below 500 m. The R0 bar provides a color gradient showing the colors, which correspond to specific R0 values. The numbers around the maps depict the Northings and Eastings associated with the domain. Black stippling in maps b and c shows the spatial distribution of statistically significant changes in R0 between each future time period (2011–2040 and 2041–2070) and the historical time period 1971–2000. Statistical significance was calculated using the Kolmogorov-Smirnov test.

The northern expansion of the ticks’ geographic range as well as the increases in the multimodel mean R0 values were of a similar magnitude across all RCP scenarios for the period 2011 to 2040. During this time period, there were no statistically significant differences between the R0 multimodel mean estimates under RCP2.6 versus RCP4.5 or under RCP4.5 versus RCP6.0. However, statistically significant differences were found between RCP6.0 and RCP8.5 (Figure S5). This is consistent with the IPCC WGI report (IPCC 2013), in which climate model outputs using all RCP scenarios project similar temperature increases for the first two decades after 2005. Similarly, the climate model variability was of a similar magnitude across the four RCP scenarios. Over this time period, the range of the 95% confidence intervals began to display a trend, expanding both with decreasing latitudes and also in accordance with increasing R0 estimates (Figures 2e and 3e and Figures S3e and S4e).

The possible geographic range of ticks projected for 2041 to 2070 covers most of Atlantic Canada, Manitoba, and areas south of 55°N in Quebec for all RCP scenarios. Under RCP4.5, RCP6.0, and RCP8.5, the ticks’ geographic range is projected to expand into areas of Newfoundland and Labrador. The highest R0 values and the greatest range expansion were obtained using the RCP8.5, followed closely by the RCP4.5. The lowest R0 values and the smallest range expansion were found using the RCP2.6. Because the maps presented here show data up to 2070, they do not capture R0 estimates for the latter 30 years of the century when DD>0°C and therefore R0 increases become greater under RCP6.0 compared with those under RCP4.5 (Figure S6). This being said, the differences in the estimated R0 mean values for RCP2.6 versus RCP4.5 and for RCP4.5 versus RCP6.0 were small; no statistically significant differences were identified for Canada (Figure S5). Differences between the R0 estimates for RCP6.0 and RCP8.5 were again statistically significant in the major part of the geographic domain (Figure S5).

Although the range of the 95% confidence intervals for the multimodel ensemble mean increased with time, our results suggest that, even when the lowest bound of the climate model variability is considered, the northern expansion of the ticks’ latitudinal limit continues to be notable for all CMIP5 RCP simulations. For example, the hindcasted climate (1971 to 2000) at 47°N, 80°W was not suitable for tick establishment, resulting in R0 below 1. However, by 2099, the mean R0 under RCP2.6 and RCP8.5 was estimated at approximately 4.0±2.0 and 6.0±3.0, respectively (Figure S3c and Figure 3c). This demonstrates that the lowest value within the 95% confidence interval for both the RCP2.6 and RCP8.5 simulations results in a change in the R0 from below 1 to above 1 in at least some locations in southern Canada.

R0 values in Nova Scotia and Southern Ontario

The predicted evolution of R0 values for the coming century for Nova Scotia and southern Ontario (south of 50°N) are presented in Figure 4. The results for northern Ontario (north of 50°N) are presented in Figure 5.

Two shaded line graphs show R0 values (y-axis) for the years 1914, 1940, 1967, 1993, 2010, 2046, 2073, and 2099 (x-axis) for Southern Ontario and Nova Scotia. The simulation lines indicate historical period, RCP26, RCP45, RCP60, and RCP80.

Figure 4. Annual estimates of the R0 multimodel mean and the climate model variability (95% confidence interval) for Ontario south of 50°N (left) and Nova Scotia (right).

Shaded line graph shows R0 values (y-axis) for the years 1914, 1940, 1967, 1993, 2010, 2046, 2073, and 2099 (x-axis) for Northern Ontario. The simulation lines indicate historical period, RCP26, RCP45, RCP60, and RCP80.

Figure 5. Annual estimates of the R0 multimodel mean and the model variability (95% confidence interval) for Ontario north of 50°N.

For the historical period in southern Ontario and Nova Scotia, estimated R0 multimodel mean values experienced some annual fluctuation, but they remained relatively stable up to the early 1970s. R0 multimodel mean values then began to increase with surface air temperature warming in the 1970s (IPCC 2013). By 2005, the annual multimodel mean R0 values for Nova Scotia and southern Ontario were estimated at approximately 2.2±1.0 and 1.8±0.72, respectively.

Between 2005 and 2040, the evolution of the R0 values was of similar magnitude for all four RCP simulations. All RCP simulations displayed R0 mean values >1 for Nova Scotia and southern Ontario by 2040. Only after 2040 did the annual R0 multimodel means and the climate model variability become increasingly dependent on the specified RCP simulation. By 2099, estimated R0 mean values were the greatest under RCP8.5 simulations and the lowest under RCP2.6 simulations. In Nova Scotia, the annual R0 mean values were projected to range from 1.4 to 5.2 (under RCP2.6) and from 3.1 to 16 (under RCP8.5). In southern Ontario, the annual R0 mean values were projected to range from 0.92 to 5.0 (under RCP2.6) and from 2.0 to 16 (under RCP8.5). By 2099, R0 values > 1 fell within the 95% confidence interval for all four RCP simulations in Nova Scotia and for all simulations with the exception of RCP2.6 in southern Ontario (Figure 4).

For northern Ontario, historical mean R0 values were below 1 and remained so until the 2040s for all RCPs. From the 2040s onward, R0 mean values increased up to 4 using RCP8.5. However, for the other RCPs, R0 remained close to 1 by 2099 (Figure 5).


Changes in R0 that are projected to occur under future climate conditions will likely influence public health by altering the distribution and increasing the magnitude of Lyme disease risk across Canada. Climate change may lead to a larger area of Canada experiencing climate suitable for tick establishment, inducing a change from R0<1 to R0>1. Climate change may also affect the speed at which ticks’ range expands and may potentially increase tick abundance (Leighton et al. 2012; Ogden et al. 2014b; Stafford et al. 1998).

Our study attempts to build upon the findings of Ogden et al. (2014b) by exploring the potential impacts that the full range and most up-to-date climate models and emission scenarios may have on future projected R0 estimates for I. scapularis. To accomplish this, R0 projections were estimated using a set of thirteen CMIP5 climate models and all four RCP scenarios. CMIP5 GCMs take into account more detailed and numerous climate processes performed at higher resolutions than the CMIP3 models employed by Ogden et al. (2014b). Such improvements enhance the projections of climate change and therefore the estimates of R0. Additionally, by employing a large ensemble of GCMs and more complete and numerous future scenarios, our work provides a more comprehensive estimate of the possible impact of climate change on Lyme disease expansion over Canada.

The RCP8.5 is the most similar emission scenario in terms of the trajectory of GHG emissions over the 21st century to the SRES A2 scenario, which was used in Ogden et al. (2014b). Although both scenarios are similar, the RCP8.5 scenario leads to a wider range expansion and larger increases in the R0 by 2099. These differences may be caused by improvements in the RCP8.5, such as updates of base-year land-cover statistics, updates in the resource inventory of the Agricultural Modeling Framework, and aggregation of the pasture and natural grasslands in the land-use projections (Riahi et al. 2011). Our results validate the findings of Ogden et al. (2006) by suggesting that there is a general consensus amongst the climate models and emission scenarios that the geographic range suitable for tick establishment will expand northward and R0 increases will occur under future climate conditions. Our results also demonstrate that climate change may affect the R0 of I. scapularis despite the existing uncertainties in climate change projections. For example, Figure 4 indicates R0 multimodel mean values were estimated to increase under all four RCP scenarios (Figure 4). R0 values >1 fell within the 95% confidence interval for all four RCP simulations in Nova Scotia and for all simulations with the exception of RCP2.6 in southern Ontario. Relative to the historical period, 1971 to 2000, R0 multimodel mean increases estimated for the periods 2011 to 2040 and 2041 to 2070 were statistically significant under each RCP scenario for the entire province of Nova Scotia and Ontario south of approximately 47°N (Figures 3 and 4 and Figures S3 and S4).

Climate change projections made by RCP climate simulations were referenced at the Paris Climate Conference. Mitigation associated with RCP2.6 is the most likely of the four RCPs to reach the governments’ goal to limit warming to 2°C or 1.5°C, relative to preindustrial levels (IPCC 2013). Our results suggest that increases in R0 values and expansions of the ticks’ range in Canada are statistically significant even under RCP2.6, and, thus, Lyme disease risk will likely emerge in Canada even if the Paris Agreement’s goal is achieved. Although the R0 increases projected under all four RCP scenarios for the periods 2011 to 2040 and 2041 to 2070, relative to the historical period 1971 to 2000, were found to be statistically significant, the difference in the magnitude of these increases was statistically significant only between RCP8.5 and RCP6.0 (Figure S5). This suggests that mitigation that results in reductions in both emissions and climate warming towards the three lower RCP scenarios may result in public health benefits starting in the 2030s; However, these benefits would not be realized if current emissions continue, as represented by the RCP8.5 scenario (Edenhofer 2014; IPCC 2013). Nonetheless, this being said, a large percentage of the Canadian population will be inhabiting or pursuing outdoor activities in southern regions predicted to be climatically suitable for I. scapularis populations in the early decades of this century (Leighton et al. 2012). Thus, mitigation may not have a marked impact on Lyme disease risks for the Canadian population in these regions. For those living further north, mitigation efforts that reduce emissions from those simulated in RCP8.5 to those in RCP6.0 may have much greater impact on whether or not these populations become at risk from Lyme disease (Figure 5).

Due to the use of GCM simulations, our results provide simple estimates of R0 for I. scapularis. Although it is widely accepted that GCMs are adequate tools for simulating large-scale climate processes, these simulations likely result in large uncertainty in their predictive ability on a local scale (Semenov and Stratonovitch 2010; IPCC 2013). The horizontal resolution of the climate models used in this study ranged from 100km×89km to 312km×20km. These climate models do not account for the small-scale climate processes or terrain and land-use changes that may affect tick populations. Furthermore, the physics and parameterizations of climate processes vary among the CMIP5 GCMs used in this study. For example, although all CMIP5 modeling groups are supplied with the same specified land-use changes, the dynamic vegetation differs among each modeling group. This can lead to variable representations of land surface characteristics in each GCM.

Further uncertainties in our projections may stem from Equation 1. By applying this equation, we are implicitly assuming that factors that play an important role in tick population establishment (i.e., host densities, mortality rates, leaf litter, precipitation extremes, etc.) are spatially and temporally homogeneous. A number of parameters in the I. scapularis model used for this study have assumed values, and the influence of these values on model outcome has been explored in local and global sensitivity analyses (Ogden et al. 2005; Ogden et al. 2014b). Empirical studies have, however, demonstrated that, in southeastern Canada, the threshold temperature conditions for R01 for I. scapularis, obtained by model simulations using these values, are accurate (Gabriele-Rivet et al. 2015). Clearly though, two factors in addition to suitable temperature conditions must be present for I. scapularis populations to persist: woodland habitat to provide habitat for hosts and refuges for the ticks from extremes of temperature and host (particularly deer) densities above thresholds for tick population persistence. The patchy occurrence of these factors is likely responsible for the patchy emergence of I. scapularis at a local level in Canada (Gabriele-Rivet et al. 2015; Clow et al. 2016); however, at a national level, empirical data suggest that temperature conditions have been the main determinant of the rates and trajectories of I. scapularis spread, even when accounting for alternative factors (Leighton et al. 2012). There is no expectation that these other factors will be limiting in the future (Ogden et al. 2014b), but this cannot be ruled out. Further study is needed to monitor their effects. Our conclusions are based on the assumption that temperature is indeed the main limiting factor on northward expansion of the tick populations and will remain so in the future.

The difference between I. scapularis establishment and B. burgdorferi transmission risk should be taken into account; the ticks themselves do not threaten public health in the absence of B. burgdorferi and other tick-borne pathogens. To date, the establishment of B. burgdorferi in Canada has occurred following that of I. scapularis; however, the former lags several years behind the initial establishment of tick populations (Ogden et al. 2010). A 5-year gap likely occurs between tick and B. burgdorferi establishment in Eastern Canada while a shorter, 3-year gap occurs in Central Canada (Ogden et al. 2013). Also, tick phenology and the fitness of different B. burgdorferi strains may be affected by climate change (Ogden et al 2008; Ogden et al. 2013; Moore et al. 2014; Monaghan et al 2015), but this analysis is beyond the scope of this paper.

The results presented here provide the most detailed analysis to date of how climate change may impact the future distribution of I. scapularis at the expanding northern edge of the ticks’ range. R0 estimates were obtained using an ensemble of 13 CMIP5 GCMs and all the available RCP scenarios. These results suggest that increases in Lyme disease risk will occur in many parts of southern Canada even under the most optimistic future climate conditions, which emphasizes the need for continuous development and implementation of public health adaptation strategies.


This study has shown that all CMIP5 RCP simulations project statistically significant changes in the R0 of I. scapularis over time, as a result of temperature increases associated with climate change. The results suggest that increases in R0 multimodel mean values and northward range expansion of I. scapularis may occur even under the most optimistic future climate change scenario projected by the RCP2.6 simulations, where GHG concentrations are assumed to be greatly reduced over time as result of stringent policies (Wayne 2013). However, mitigation may have a greater impact on preventing Lyme disease risk emergence in more northern regions of Canada. This information reinforces the need for continuous planning and implementation of location-specific robust Lyme disease adaptation strategies to ensure that the health of the Canadian population is protected from the range of possibilities by which climate change may affect I. scapularis and consequently alter Lyme disease risk in Canada.


This work was supported by grants from the Natural Sciences and Engineering Research Council of Canada Discovery Grant (NSERC-DG) program, a NSERC-CREATE award Training Program in Climate Sciences, the Atlantic Computational Excellence Network (ACEnet), and the Canadian Foundation for Innovation to H. Beltrami (NSERC DG 140576948). H. Beltrami holds a Canada Research Chair. M.M., A.G.G. and F.J.C.V. are funded by the NSERC-CREATE Training Program in Climate Sciences based at St. Francis Xavier University. M.M. is also funded by the Nova Scotia Graduate Scholarship. The data for this paper are available at the CMIP5 data center at


  • Adrion ER, Aucott J, Lemke KW, Weiner JP. 2015. Health care costs, utilization and patterns of care following Lyme disease. PLoS ONE 10:e0116767, PMID: 25650808, doi:10.1371/journal.pone.0116767. Crossref, MedlineGoogle Scholar
  • Baldwin DJB, Desloges JR, Band LE. 2001. Physical geography of Ontario. Ecology of a Managed Terrestrial Landscape: Patterns and Processes of Forest Landscapes in Ontario, Perera AH, Euler DL, Thompson ID, eds. Vancouver: UBC Press, 12–29. Google Scholar
  • Bouchard C, Leonard E, Koffi JK, Pelcat Y, Peregrine A, Chilton N, et al. 2015. The increasing risk of Lyme disease in Canada. Can Vet J 56(7):693–699, PMID: 26130829. MedlineGoogle Scholar
  • Buis A. 2011. Climate change may bring big ecosystem changes. [accessed 1 December 2015]. Google Scholar
  • Bush EJ, Loder JW, Mortsch LD, Cohen SJ. 2014. An overview of Canada's changing climate. In: Canada in a Changing Climate: Sector Perspectives on Impacts and Adaptation, Warren FJ, Lemmen DS, eds. Ottawa: Government of Canada, 23–64. Google Scholar
  • Clow K, Ogden NH, Lindsay LR, Michel P, Pearl D, Jardine C. 2016. Distribution of ticks and the risk of Lyme disease and other tick-borne pathogens of public health significance in Ontario, Canada. Vector Borne Zoonotic Dis 16(4):215–222, PMID: 26870937, doi:10.1089/vbz.2015.1890. Crossref, MedlineGoogle Scholar
  • CMIP-WGCM-WGSIP Decadal Climate Prediction Panel. 2011. Data and bias correction for decadal climate predictions. International CLIVAR Office, CLIVAR Publication Series, 150. [accessed 12 May 2017]. Google Scholar
  • Diuk-Wasser MA, Hoen AG, Cislo P, Brinkerhoff R, Hamer SA, Rowland M, et al.2012. Human risk of infection with Borrelia burgdorferi, the Lyme disease agent, in eastern United States. Am J Trop Med Hyg 86:320–327, PMID: 22302869, doi:10.4269/ajtmh.2012.11-0395. Crossref, MedlineGoogle Scholar
  • Donner LJ, Wyman BL, Hemler RS, Horowitz LW, Ming Y, Zhao M, et al.2011. The dynamical core, physical parameterizations, and basic simulation characteristics of the atmospheric component AM3 of the GFDL global coupled model CM3. J Climate 24:3484–3519, doi:10.1175/2011JCLI3955.1. CrossrefGoogle Scholar
  • Dufresne JL, Foujols M-A, Denvil S, Caubel A, Marti O, Aumont O, et al.2013. Climate change projections using the IPSL-CM5 earth system model: From CMIP3 to CMIP5. Clim Dyn 40(9):2123–2165, doi:10.1007/s00382-012-1636-1. CrossrefGoogle Scholar
  • Dunne JP, John JG, Adcroft AJ, Griffies SM, Hallberg RW, Shevliakova E, et al.2012. GFDL's ESM2 global coupled climate-carbon earth system models. Part I: Physical formulation and baseline simulation characteristics. J Climate 25:6646–6665, doi:10.1175/JCLI-D-11-00560.1. CrossrefGoogle Scholar
  • Dunn-Sigouin E, Son S. 2013. Northern hemisphere blocking frequency and duration in the CMIP5 models. J Geophys Res Atmos 118:1179–1188, doi:10.1002/jgrd.50143. CrossrefGoogle Scholar
  • Edenhofer O. 2014. Future Pathway for Adaptation, Mitigation and Sustainable Development. COP 20. 2 December 2014, Lima, Peru. Google Scholar
  • Eisen RJ, Eisen L, Beard CB. 2016. County-scale distribution of Ixodes scapularis & Ixodes pacificus (acari: Ixodidae) in the continental United States. J Med Entomol 53:1–38, doi:10.1093/jme/tjv237. Crossref, MedlineGoogle Scholar
  • Fall S, Niyogi D, Gluhovsky A, Pielke RA, Kalnay E, Rochon G. 2010. Impacts of land use land cover on temperature trends over the continental United States: assessment using the North American Regional Reanalysis. Int J Climatol 30:1980–1993, doi:10.1002/joc.1996. CrossrefGoogle Scholar
  • Gabriele-Rivet V, Arsenault J, Badcock J, Cheng A, Edsall J, Goltz J, et al.2015. Different ecological niches for ticks of public health significance in Canada. PLoS One 10(7):e0131282, PMID: 26131550, doi:10.1371/journal.pone.0131282. Crossref, MedlineGoogle Scholar
  • Gent PR, Danabasoglu G, Donner LJ, Holland MM, Hunke EC, Jayne SR, et al.2011. The community climate system model version 4. J Climate 24:4973–4991, doi:10.1175/2011JCLI4083.1. CrossrefGoogle Scholar
  • Government of Ontario. 2015. About Ontario. [accessed 5 October 2015]. Google Scholar
  • Hinckley AF, Connally NP, Meek JI, Johnson BJ, Kemperman MM, Feldman KA, et al.2014. Lyme disease testing by large commercial laboratories in the United States. Clin Infect Dis 59(5):676–681, PMID: 24879782, doi:10.1093/cid/ciu397. Crossref, MedlineGoogle Scholar
  • IPCC, 2013: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [Stocker T.F., Qin D., Plattner G.-K., Tignor M., Allen S.K., Boschung J., Nauels A., Xia Y., Bex V. and Midgley P.M. (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp. Google Scholar
  • Kanamaru H, Kanamitsu M. 2007. Fifty-seven-year California Reanalysis Downscaling at 10 km (CaRD10). Part II: comparison with North American Regional Reanalysis. J Climate 20:5572–5592, doi:10.1175/2007JCLI1522.1. CrossrefGoogle Scholar
  • Kopparla P, Fischer EM, Hannay C, Knutti R. 2013. Improved simulation of extreme precipitation in a high-resolution atmosphere model. Geophys Res Lett 40(21):5803–5808, doi:10.1002/2013GL057866. CrossrefGoogle Scholar
  • Leighton PA, Koffi JK, Pelcat Y, Lindsay LR, Ogden NH. 2012. Predicting the speed of tick invasion: an empirical model of range expansion for the Lyme disease vector Ixodes scapularis in Canada. J Appl Ecol 49(2):457–464, doi:10.1111/j.1365-2664.2012.02112.x. CrossrefGoogle Scholar
  • Masui T, Matsumoto K, Hijioka Y, Kinoshita T, Nozawa T, Ishiwatari S, et al.2011. An emission pathway for stabilization at 6 Wm-2 radiative forcing. Clim Change 109:59–76, doi:10.1007/s10584-011-0150-5. CrossrefGoogle Scholar
  • Mearns LO. 2010. The drama of uncertainty. Clim Change 100:77–85, doi:10.1007/s10584-010-9841-6. CrossrefGoogle Scholar
  • Meehl GA, Washington WM, Arblaster JM, Hu A, Teng H, Kay JE, et al.2013. Climate change projections in CESM1 (CAM5) compared to CCSM4. J Climate 26:6287–6308, doi:10.1175/JCLI-D-12-00572.1. CrossrefGoogle Scholar
  • Mesinger F, DiMeso G, Kalnay E, Mitchell K, Shafran PC, Ebisuzaki W, et al.2006. North American Regional Reanalysis. Bull Amer Meteor Soc 87:343–360, doi:10.1175/BAMS-87-3-343. CrossrefGoogle Scholar
  • Mieville A, Granier C, Liousse C, Guillaume B, Mouillot F, Lamarque J-F, et al.2010. Emissions of gases and particles from biomass burning during the 20th century using satellite data and an historical reconstruction. Atmos Environ 44(11):1469–1477, doi:10.1016/j.atmosenv.2010.01.011. CrossrefGoogle Scholar
  • Monaghan AJ, Moore SM, Sampson KM, Beard CB, Eisen RJ. 2015. Climate change influences on the annual onset of Lyme disease in the United States. Ticks Tick Borne Dis 6:615–622, PMID: 26025268, doi:10.1016/j.ttbdis.2015.05.005. Crossref, MedlineGoogle Scholar
  • Moore SM, Eisen RJ, Monaghan A, Mead P. 2014. Meteorological influences on the seasonality of Lyme disease in the United States. Am J Trop Med Hyg 90:486–496, PMID: 24470565, doi:10.4269/ajtmh.13-0180. Crossref, MedlineGoogle Scholar
  • Ogden NH, Bigras-Poulin M, O’Callaghan CJ, Barker IK, Lindsay LR, Maarouf A, et al.2005. A dynamic population model to investigate effects of climate on geographic range and seasonality of the tick Ixodes scapularis. Int J Parasitol 35(4):375–389, PMID: 15777914, doi:10.1016/j.ijpara.2004.12.013. Crossref, MedlineGoogle Scholar
  • Ogden NH, Bouchard C, Kurtenbach K, Margos G, Lindsay LR, Trudel L, et al.2010. Active and passive surveillance and phylogenetic analysis of Borrelia burgdorferi elucidate the process of Lyme disease risk emergence in Canada. Environ Health Perspect 118(7):909–914, PMID: 20421192, doi:10.1289/ehp.0901766. LinkGoogle Scholar
  • Ogden NH, Koffi JK, Pelcat Y, Lindsay LR. 2014a. Environmental risk from Lyme disease in central and eastern Canada: a summary of recent surveillance information. Can Commun Dis Rep 40:74–82. Crossref, MedlineGoogle Scholar
  • Ogden NH, Lindsay LR, Morshed M, Sockett PN, Artsob H. 2009. The emergence of Lyme disease in Canada. Can Med Assoc J 180:1221–1224, PMID: 19506281, doi:10.1503/cmaj.080148. Crossref, MedlineGoogle Scholar
  • Ogden NH, Lindsay LR, Leighton PA. 2013. Predicting the rate of invasion of the agent of Lyme disease Borrelia burgdorferi. J Appl Ecol 50:510–518, doi:10.1111/1365-2664.12050. CrossrefGoogle Scholar
  • Ogden NH, Maarouf A, Barker IK, Bigras-Poulin M, Lindsay LR, Morshed MG, et al.2006. Climate change and the potential for range expansion of the Lyme disease vector Ixodes scapularis in Canada. Int J Parasitol 36(1):63–70, PMID: 16229849, doi:10.1016/j.ijpara.2005.08.016. Crossref, MedlineGoogle Scholar
  • Ogden NH, Radojevic M, Wu X, Duvvuri VR, Leighton PA, Wu J. 2014b. Estimated effects of projected climate change on the basic reproductive number of the Lyme disease vector Ixodes scapularis. Environ Health Perspect 122:631–638, doi:10.1289/ehp.1307799. LinkGoogle Scholar
  • Ogden NH, St-Onge L, Barker IK, Brazeau S, Bigras-Poulin M, Charron DF, et al.2008. Risk maps for range expansion of the Lyme disease vector, Ixodes scapularis, in Canada now and with climate change. Int J Health Geogr 7:24. PMID: 18498647, doi:10.1186/1476-072X-7-24. Crossref, MedlineGoogle Scholar
  • PHAC. 2015. Lyme disease and other tick-borne diseases: information for healthcare professionals. [accessed 4 September 2015]. Google Scholar
  • Preston BL, Dow K, Berkhout F. 2013. The climate adaptation frontier. Sustainability 5:1011–1035, doi:10.3390/su5031011. CrossrefGoogle Scholar
  • Riahi K, Rao S, Krey V, Cho C, Chirkov V, Fischer G, et al.2011. RCP 8.5—a scenario of comparatively high greenhouse gas emissions. Clim Change 109:33–57, doi:10.1007/s10584-011-0149-y. CrossrefGoogle Scholar
  • Schultz MG, Heil A, Hoelzemann JJ, Spessa A, Thonicke K, Goldammer JG, et al.2008. Global wildland fire emissions from 1960 to 2000. Global Biogeochem Cycles 22(2), doi:10.1029/2007GB003031. CrossrefGoogle Scholar
  • Semenov MA, Stratonovitch P. 2010. Use of multi-model ensembles from global climate models for assessment of climate change impacts. Clim Res 41:1, doi:10.3354/cr00836. CrossrefGoogle Scholar
  • Stafford KC, 3rd, Cartter ML, Magnarelli LA, Ertel SH, Mshar PA. 1998. Temporal correlations between tick abundance and prevalence of ticks infected with Borrelia burgdorferi and increasing incidence of Lyme disease. J Clin Microbiol 36(5):1240–1244. MedlineGoogle Scholar
  • Taylor KE, Stouffer RJ, Meehl GA. 2012. An overview of CMIP5 and the experiment design. Bull Amer Meteor Soc 93:485–498, doi:10.1175/BAMS-D-11-00094.1. CrossrefGoogle Scholar
  • Thomson AM, Calvin KV, Smith SJ, Kyle GP, Volke A, Patel P, et al.2011. RCP4. 5: a pathway for stabilization of radiative forcing by 2100. Clim Change 109:77–94, doi:10.1007/s10584-011-0151-4. CrossrefGoogle Scholar
  • van der Werf GR, Randerson JT, Giglio L, Collatz GJ, Kasibhatla PS, Arellano AF. 2006. Interannual variability in global biomass burning emissions from 1997 to 2004. Atmos Chem Phys 6:3423–3441, doi:10.5194/acp-6-3423-2006. CrossrefGoogle Scholar
  • Van Vuuren DP, Edmonds J, Kainuma M, Riahi K, Thomson A, Hibbard K, et al.2011. The representative concentration pathways: an overview. Clim Change 109:5–31, doi:10.1007/s10584-011-0148-z. CrossrefGoogle Scholar
  • Wang J, Tsang WW, Marsaglia G. 2003. Evaluating Kolmogorov's distribution. J Stat Soft 8:, doi:10.18637/jss.v008.i18. CrossrefGoogle Scholar
  • Watanabe M, Suzuki T, O’ishi R, Komuro Y, Watanabe S, Emori S, et al.2010. Improved climate simulation by MIROC5: mean states, variability, and climate sensitivity. J Climate 23:6312–6335, doi:10.1175/2010JCLI3679.1. CrossrefGoogle Scholar
  • Watanabe S, Hajima T, Sudo K, Nagashima T, Takemura T, Okajima H, et al.2011. MIROC-ESM 2010: model description and basic results of CMIP5-20c3m experiments. Geosci Model Dev 4:845–872, doi:10.5194/gmd-4-845-2011. CrossrefGoogle Scholar
  • Wayne G. 2013. [accessed 2 April 2017]. Google Scholar
  • Willison J, Robinson WA, Lackmann GM. 2013. The importance of resolving mesoscale latent heating in the North Atlantic storm track. J Atmos Sci 70:2234–2250, doi:10.1175/JAS-D-12-0226.1. CrossrefGoogle Scholar
  • Wu T, Song L, Li W, Wang Z, Zhang H, Xin X, et al.2014. An overview of BCC climate system model development and application for climate change studies. Acta Meteorol Sin 28(1):34–56, doi:10.1007/s13351-014-3041-7. CrossrefGoogle Scholar
  • Wu X, Duvvuri VR, Lou Y, Ogden NH, Pelcat Y, Wu J. 2013. Developing a temperature-driven map of the basic reproductive number of the emerging tick vector of Lyme disease Ixodes scapularis in Canada. J Theor Biol 319:50–61, doi:10.1016/j.jtbi.2012.11.014. Crossref, MedlineGoogle Scholar
  • Yukimoto S, Hosaka M, Tanaka T. 2012. A new global climate model of the meteorological research institute: MRI-CGCM3: model description and basic performance. J Meteor Soc Japan 90A:23–64, doi:10.2151/jmsj. CrossrefGoogle Scholar
  • Zhang X, Meltzer MI, Peña CA, Hopkins AB, Wroth L, Fix AD. 2006. Economic impact of Lyme disease. Emerging Infect Dis 12:653–660, doi:10.3201/eid1204.050602. Crossref, MedlineGoogle Scholar