Estimated Effects of Projected Climate Change on the Basic Reproductive Number of the Lyme Disease Vector Ixodes scapularis

Background: The extent to which climate change may affect human health by increasing risk from vector-borne diseases has been under considerable debate. Objectives: We quantified potential effects of future climate change on the basic reproduction number (R0) of the tick vector of Lyme disease, Ixodes scapularis, and explored their importance for Lyme disease risk, and for vector-borne diseases in general. Methods: We applied observed temperature data for North America and projected temperatures using regional climate models to drive an I. scapularis population model to hindcast recent, and project future, effects of climate warming on R0. Modeled R0 increases were compared with R0 ranges for pathogens and parasites associated with variations in key ecological and epidemiological factors (obtained by literature review) to assess their epidemiological importance. Results: R0 for I. scapularis in North America increased during the years 1971–2010 in spatio-temporal patterns consistent with observations. Increased temperatures due to projected climate change increased R0 by factors (2–5 times in Canada and 1.5–2 times in the United States), comparable to observed ranges of R0 for pathogens and parasites due to variations in strains, geographic locations, epidemics, host and vector densities, and control efforts. Conclusions: Climate warming may have co-driven the emergence of Lyme disease in northeastern North America, and in the future may drive substantial disease spread into new geographic regions and increase tick-borne disease risk where climate is currently suitable. Our findings highlight the potential for climate change to have profound effects on vectors and vector-borne diseases, and the need to refocus efforts to understand these effects. Citation: Ogden NH, Radojević M, Wu X, Duvvuri VR, Leighton PA, Wu J. 2014. Estimated effects of projected climate change on the basic reproductive number of the Lyme disease vector Ixodes scapularis. Environ Health Perspect 122:631–638; http://dx.doi.org/10.1289/ehp.1307799


Table of Contents Page
Variation in temperature and R 0 amongst Canadian sites 2 Table S1. Basic characteristics of sites used in this study 4 Figure S1. The locations of sites used in the study 5 Figure S2. Model-estimated R 0 and DD > 0°C using observed daily temperature 5 Figure S3. Climate classification tree for the sites 6 Figure S4. Daily normal of ANUSPLIN daily maximum temperatures for each cluster 6 Model sensitivity analysis 7 Table S2. The parameters used in sensitivity analysis 9 Validation of climate model output 11 Table S3. Basic information about RCM simulations 13 Figure S7. Mean DD > 0°C for each cluster using observed and model temperature 14 Sites at which we estimated R 0 were selected in previous studies (Ogden et al. 2005;Wu et al. 2013) ( Table S1) and form rough south-north transects in eastern Canada ( Figure S1). They aimed to capture the range of climatic zones in this region (Natural Resources Canada 2013a), where I. scapularis is currently expanding its range, which is due to complex orography, water proximity and forest structure (Natural Resources Canada 2013b) that influence local temperature regimes. We re-investigated the degree of climate variability amongst these sites in order to come to a compromise between inter-site variability and the need for some simplification for data presentation. In earlier studies (Ogden et al. 2005) the best index of the suitability of seasonally-variable temperature, at a particular location, for I.
scapularis survival was annual cumulative degree-days above 0°C (DD > 0°C) because temperaturedependent intersstadial tick development (a key variable through which temperature affects R 0 ) does not occur below 0°C (Ogden et al. 2004). We grouped the sites into five clusters on the basis of geographic proximity and climatological similarity. Both DD > 0°C from observed data and R 0 estimated ANUSPLIN data from 1971 to 2010 at the sites in each cluster are shown in Figure S2. The climatological similarity of sites within the clusters was supported by estimation of the Euclidian Distance of DD > 0°C calculated from ANUSPLIN data from 1971 to 2010 as a measure of similarity and Ward's criteria (Derry 2008;Munoz-Diaz and Rodrigo 2004;Romesburg 1990) for linking grid points with the highest similarity ( Figure S3). The clusters were: 1) Southern Ontario -the north shore of Lake Erie, 2) Huron Ontario -Ontario adjacent to the south-east extent of Lake Huron, 3) Upper Southern Ontario -Ontario adjacent to the eastern extent of Lake Huron close to Georgian Bay, 4) South-Western Quebec, and 5) Boreal region -the most northern sites in both Ontario and Quebec. The significance of variations amongst clusters in the relationship between modelled R 0 and DD > 0°C was tested in a linear regression model (STATA SE11.0: Statacorp, College Station, Tx). R 0 estimated by the tick model was the outcome while DD > 0°C and clusters (as a categorical variable) were explanatory variables. Adjusted for DD > 0°C (coefficient = 2.1E-03, SE = 4.5E-05, P < 0.001), R 0 varied 2 significantly amongst clusters (F = 163.4, df = 4, P < 0.001) suggesting that presentation of simulations from all clusters separately would be more rigorous than combining all data. Differences amongst sites and clusters in the relationship between R 0 and DD > 0°C can be explained by differences in seasonal patterns of daily temperatures: R 0 is slightly higher in South-Western Quebec than Huron Ontario even though DD > 0°C values are similar because, despite longer winters (when temperatures are below 0°C), mid-summer temperatures in South-Western Quebec are higher resulting in faster tick development ( Figure S4).     and 5. Boreal region. The Y-axis is the Euclidian distance as a measure of climate similarity.
Daily maximum temperature (°C) Figure S4. Daily normal of ANUSPLIN daily maximum temperatures for each cluster for the period 1971 to 2010.

Model sensitivity analysis
The I. scapularis model is sensitive to a greater or lesser degree to changes in model parameter values, particularly changes in temperature (which is explicitly modelled here), but also host abundance and offhost tick mortality (Wu et al. 2013). Host abundance (deer and mice in this model) will vary from location to location on both local and national scales depending on the specific suitability of habitat in each location for those species. Off-host tick mortality varies with the capacity of the environment of the surface layers of the soil, which provide the environment for ticks during the long period of development from one instar to the next, to protect the ticks from degrees of temperature, desiccation or flooding that result in greater mortality (Ogden et al. 2006). This results in many environments being unsuitable for tick survival, and partly explains the association between I. scapularis and woodland habitats. Woodland habitats do vary in their capacity to provide refugia for developing ticks, which is reflected in different mortality rates of ticks undergoing intersstadial development (Ogden et al. 2006).
To assess the sensitivity of our projected changes in R 0 to all model variables we performed local and global sensitivity analyses. For these analyses, monthly temperature data for 1971-2070 obtained from the regional climate model CRCM 4.2.3 were used to determine the temperature-dependent tick intersstadial development rates and tick activity rates (Wu et al. 2013). All model variables were incorporated in the sensitivity analyses (Table S2). All baseline values were those in (Wu et al. 2013), and all simulations were performed in Matlab R2010a (MathWorks, Natick, MS, USA).

Local sensitivity analysis
Local sensitivity analysis was used to quantify the impact of variation of individual parameters on estimated R 0 (Benton and Grant 1999). Local sensitivity analyses were performed for data from one site of each of the five Canadian clusters: Point Pelee, London, Wiarton, Hemmingford and Timmins. All parameters were varied one at a time by 5% (above and below baseline values, with all other parameter values set at baseline values) and 3-year mean R 0 values for each site for each year from 1972 to 2069 7 were obtained. A comparative index S i was used to assess the sensitivity of the change in R 0 , from one time slice to another, to a 5% change in the i th parameter: where r i 2 is the average value for R 0 for the later time period, and r i 1 is the average value for R period. This analysis indicated that the degree of change in R 0 from one time period to the next was most sensitive to development rates and mortality rates of eggs and larvae (Fig S5), although the effects of changing these parameter values by 5% had relatively modest effects on increases in R 0 from one time period to another (changes of 0.15 or less) (Fig S5). Simultaneously changing mortality rates of all offhost tick stages by 5% (as would be expected to occur in different woodland types: Ogden et al. 2006) had a larger effect on increases in R 0 (changes of 0.28 or less) from one time period to another ( Fig S5).

Global sensitivity analysis
We used Latin Hypercube sampling (LHS) to estimate the possible range of R 0 caused by simultaneous variation of multiple parameters. LHS is an efficient and widely used statistical sampling method of uncertainty analysis that treats each parameter as a random variable within a pre-determined range (Matser et al. 2009;Marino et al. 2008). In this global sensitivity analysis, we considered only one site  Note that U(--,--) indicates the uniform distribution in the process of allocation using LHS.  Table S2) numbers of deer and rodents, density-dependent effects on fecundity, duration of larval hardening, development rates of larvae, nymphs and adults, mortality rates of eggs, and hardening and questing larvae, density dependent mortality of feeding larvae, engorged larva and questing nymph mortality, density-dependent mortality of feeding nymphs, mortality of engorged nymphs and questing adults, density-dependent mortality of feeding adults and mortality rates of engorged adults. The final bars show the combined effects of changes in mortality rates of all off-host ticks simultaneously.

Validation of climate model output
To validate CRCM4.2.3 temperature output used in modelling R 0 , and to compare with an ensemble of climate models (Table S3), we compared DD > 0°C output from the models with observed (ANUSPLIN) data. In addition to CRCM4.2.3 (Laprise et al. 1997;Music et al. 2007), two RCM runs were available from the North American Regional Climate Change Assessment Program (NARCCAP; Mearns  A2 emission scenario (Nakicenovic and Swart 2000). The A2 scenario, which is described as mid-high Green-House-Gas (GHG) emission scenario with continuously increasing global population and regionally oriented economic growth, was chosen because of the availability of RCM model output using this emissions scenario. Daily temperatures from output of each RCM and GCM were aggregated onto common grid projection with ANUSPLIN at approximately 45 km horizontal grid resolution, true at 60°N (Heum et al. 2012 the accuracy of temperature prediction ( Fig S7) and bias-corrected temperature data were used in estimations of R 0 by the I. scapularis model. R 0 estimates for the Canadian sites using observed temperature data and climate models output are presented in Fig S8.   under recent climate for the 30 Canadian sites over 40 years of observed temperature data  to obtain a relationship between R 0 and DD > 0°C. We only have R 0 values for the study sites but we have current and future projected DD > 0°C data for most of North America. Using the 'Basic fitting' tool of Matlab® R2010a (MathWorks, Natick, MS, USA), we obtained two regression functions that