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

BACKGROUND: 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.
OBJECTIVES: 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.
METHODS: 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.
RESULTS: Increases in the multimodel mean 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.
CONCLUSION: 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. https://doi.org/10.1289/EHP57.

Introduction

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 () 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 values. describes the propensity for a parasite or microparasite to survive and be propagated. Conditions that allow are those under which the tick population can persist, while conditions under which 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 of I. scapularis. Ogden et al. (2014b) produced maps of estimated 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 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 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 of I. scapularis, b) to quantify and communicate the implications of climate model variability on the projected values, and c) to provide 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.

Methods

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.

Model	Resolution (lat°×lon°)	Reference
CCSM4.0	0.9×1.25 (100km×89km)	Gent et al. (2011)
CESM1-CAM5	0.9×1.25 (100km×89km)	Meehl et al. (2013)
MIROC-ESM	2.81×2.81 (312km×201km)	Watanabe et al. (2011)
MIROC-ESM-CHEM	2.81×2.81 (312km×201km)	Watanabe et al. (2011)
MIROC5	1.4×1.4 (156km×100km)	Watanabe et al. (2010)
GFDL-CM3	2×2.5 (222km×179km)	Donner et al. (2011)
GFDL-ESM2G	2×2.5 (222km×179km)	Dunne et al. (2012)
GFDL-ESM2M	1.13×1.13 (126km×81km)	Dunne et al. (2012)
MRI-CGCM3	1.13×1.13 (126km×81km)	Yukimoto et al. (2012)
IPSL-CM5A-LR	1.88×3.75 (209km×268km)	Dufresne et al. (2013)
IPSL-CM5A-MR	1.25×2.5 (139km×179km)	Dufresne et al. (2013)
BCC-CSM1.1	2.81×2.81 (312km×201km)	Wu et al. (2014)
BCC-CSM1.1-M	1.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 and to assess the effect of climate change on 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 . 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 resolution (Watanabe et al. 2011) to interpolate the CMIP5 GCMs and the NARR.

Estimating Mean 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):

Equation 1 describes the relationship between and annual cumulative degree-days (). is defined as the annual sum of the daily temperatures in degrees Celsius above the freezing point of water. The were computed as the accumulation of daily temperatures for each year using the bias-corrected daily surface air temperatures at each grid point and for each model simulation. From the gridded data, we calculated the spatial and temporal evolution variability of for each model across North America. 30-yearmultimodel mean 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, intersects the degree-day domain axis at two points. This results in two threshold temperature conditions for tick persistence, one at and the other at . We show the spatial distribution only for the values that were reflective of temperature conditions equal to or above . 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 and 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 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 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 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 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.

Results

Values of I. scapularis in North America

For the historical period of 1971 to 2000, the multimodel mean yielded values ≥ 3 limited to areas in Ontario below 45°N, which is slightly higher than, but close to, the values obtained using observed climate data (Ogden et al. 2014b). These lower-latitude areas of Ontario were identified as having the highest values in eastern and central Canada (Figure 2). For the historical time period, the estimated values from all GCMs were in good agreement, displaying a 95% confidence interval of (Figure 2). Estimates of the geographic range where 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 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 toward or over 2 for several consecutive years (2000 to 2005 in Ontario and 1998 to 2003 in Nova Scotia) (Figures S1 and S2).

Figure 2.

Multimodel mean 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. values were mapped for areas east of the Rocky Mountains and for elevations below 500 m. The bar provides a color gradient showing the colors, which correspond to specific 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 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.

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 estimates for the period 1971–2000. The bar provides a color gradient showing the colors, which correspond to specific 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 and consequently on of I. scapularis; increases in the multimodel mean values estimated for the period 2011 to 2040 relative to 1971 to 2000 were statistically significant under all four RCP scenarios (Figures 2 and Figures S3b and S4b). Changes in were also statistically significant under all four RCP emission scenarios between the historical period 1971 to 2000 and the period 2041 to 2070 (Figures 2 and Figures S4c and S5c).

The northern expansion of the ticks’ geographic range as well as the increases in the multimodel mean 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 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 estimates (Figures 2 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 values and the greatest range expansion were obtained using the RCP8.5, followed closely by the RCP4.5. The lowest 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 estimates for the latter 30 years of the century when and therefore increases become greater under RCP6.0 compared with those under RCP4.5 (Figure S6). This being said, the differences in the estimated 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 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 below 1. However, by 2099, the mean under RCP2.6 and RCP8.5 was estimated at approximately and , respectively (Figure S3c and Figure 3). 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 from below 1 to above 1 in at least some locations in southern Canada.

Figure 3.

Multimodel mean 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. values were mapped for areas east of the Rocky Mountains and for elevations below 500 m. The bar provides a color gradient showing the colors, which correspond to specific 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 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.

values in Nova Scotia and Southern Ontario

The predicted evolution of 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.

Figure 4.

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

Figure 5.

Annual estimates of the 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 multimodel mean values experienced some annual fluctuation, but they remained relatively stable up to the early 1970s. multimodel mean values then began to increase with surface air temperature warming in the 1970s (IPCC 2013). By 2005, the annual multimodel mean values for Nova Scotia and southern Ontario were estimated at approximately and , respectively.

Between 2005 and 2040, the evolution of the values was of similar magnitude for all four RCP simulations. All RCP simulations displayed mean values for Nova Scotia and southern Ontario by 2040. Only after 2040 did the annual multimodel means and the climate model variability become increasingly dependent on the specified RCP simulation. By 2099, estimated mean values were the greatest under RCP8.5 simulations and the lowest under RCP2.6 simulations. In Nova Scotia, the annual 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 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, 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 values were below 1 and remained so until the 2040s for all RCPs. From the 2040s onward, mean values increased up to 4 using RCP8.5. However, for the other RCPs, remained close to 1 by 2099 (Figure 5).

Discussion

Changes in 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 to . 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 estimates for I. scapularis. To accomplish this, 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 . 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 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 increases will occur under future climate conditions. Our results also demonstrate that climate change may affect the of I. scapularis despite the existing uncertainties in climate change projections. For example, Figure 4 indicates multimodel mean values were estimated to increase under all four RCP scenarios (Figure 4). values 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, 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 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 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 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 to . 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 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. 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.

Conclusions

This study has shown that all CMIP5 RCP simulations project statistically significant changes in the of I. scapularis over time, as a result of temperature increases associated with climate change. The results suggest that increases in 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.

Click here for additional data file.