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

1Considerable attention has been devoted to the possibility that climate change will exacerbate the burden of mosquito-borne diseases such as malaria and dengue, with important impacts on public health. Early assessments of the effects of climate change on malaria and dengue used simplistic models to assess possible effects of climate change on their basic reproductive numbers (R0, the universally recognized metric of the capacity of a parasite or pathogen to reproduce given particular environmental conditions). However, these assessments were criticized for giving weight to future increases in R0 whether or not such increases resulted in R0 rising above the critical threshold of > 1 for disease persistence and for being oversimplistic by only accounting for climate effects rather than the full range of nonclimatic factors that impact the occurrence of these diseases. Any impact of climate on R0 of malaria and dengue is limited by the effects of variations in human host density, mosquito control, infection prevention and treatment in humans, and human management of the environment (e.g., agriculture, forest management, logging) that affect the ecology and epidemiology of the vectors, pathogens, and diseases. Consequently, the strength of evidence for recent climate warming effects on malaria risk has been questioned and much debated.

2Many vector-borne diseases of public health significance (e.g., Lyme disease, West Nile virus) are, however, maintained in transmission cycles that involve wild animal hosts. These cycles are independent of human cases, and the spatio-temporal risk of human disease is less dependent on the direct effects of human activities than is the risk from malaria and dengue. Nevertheless, despite some assessments, the effects of climate change on vector-borne zoonoses have also been downplayed mostly on the basis of limited evidence for recent effects of climate change. Lyme disease emerged (or likely reemerged) in the northeastern United States in the late 1970s due to the expansion of tick populations, which was generally thought to have been associated with changes in land use over some decades that resulted in reforestation and expansion of the population of the deer that are key hosts for the ticks. Lyme disease is now emerging in Canada and some northern U.S. states due to the northward expansion of the geographic range of the tick vector Ixodes scapularis (I. scapularis), which is dispersed from source populations by migratory birds and terrestrial hosts.

3A mechanistic simulation model of the I. scapularis life cycle has identified temperature effects on I. scapularis population survival in order to assist in assessment of current and future on-the-ground Lyme disease risk in Canada. Prospective field studies and retrospective analyses of surveillance data on tick and pathogen emergence in southeastern Canada validated the model findings and identified temperature as a statistically significant determinant and possible driver of emergence of the tick in Canada. The I. scapularis model was modified to permit the direct calculation of R0 for I. scapularis via the next generation operator approach, which, given the universal use of R0 and its estimation for a wide range of parasites and pathogens under many different conditions, allowed comparison of R0 variations in the present study with observed variations for other parasites and pathogens.

4Here, we have estimated projected effects of climate change on R0 of an arthropod vector using a model that has been extensively ground-truthed, and we have assessed the ecological and epidemiological significance of the projected changes in R0 by comparing them with ranges of R0 values observed for other parasites and pathogens.

5We estimated R0 under current and future projected climatic conditions at 30 sites in Canada that formed two roughly south–north transects in Ontario and Quebec, two Canadian provinces where I. scapularis ticks are becoming established. These transects were chosen to capture the climate variability that exists in the region. For simplicity in data presentation, the sites were grouped into clusters [Southern Ontario, Huron Ontario, Upper Southern Ontario, South-Western Quebec, and the Boreal region (see Supplemental Material, Table S1 and Figure S1)] according to geographic proximity and similarity in temperature conditions (see Supplemental Material, “Variation in temperature and R0 amongst sites,” pp. 2–6, Figures S2–S4), and mean values for clusters are presented. We also estimated R0 for two sites in the United States where Lyme disease is endemic in the Northeast and upper Midwest, respectively: Old Lyme (Connecticut), where the human Lyme disease cases were first recognized, and Fort McCoy (Wisconsin).

6Modeling R0. The I. scapularis model is a deterministic model consisting of 12 ordinary differential equations as described by, based on the mechanistic simulation model described by. This model captures the effects of temperature on host-seeking activity and the rates of development from one life stage to the next (effects common to, but variable among, all arthropod vectors), parameterized from field and laboratory studies on I. scapularis. Mortality rates of nonfeeding I. scapularis in Canada and the northeastern United States are similar in summer and winter, presumably due to the insulating effects of the litter layer in woodland habitats. Our analyses operated on the hypothesis that the effects of ambient temperature on I. scapularis population survival are indirect via effects on temperature-dependent rates of development of ticks from one life stage to the next. The lower the temperature, the longer is the tick life cycle and, due to constant daily per capita mortality, the fewer larval ticks survive to become mated adult female ticks. At a threshold temperature, mortality outstrips reproduction and the tick populations die out (or fail to become established)—that is, at this temperature, threshold R0 falls below unity. At the latitudes under study here, effects of climate change on I. scapularis are expected to be the effect of climate warming on shortening the life cycle, resulting in increasing R0. Quadratic effects of temperature on arthropod vector life history traits are common, and quadratic effects of temperature on tick activity are included in the model. High temperatures may impact tick survival, causing a northward contraction of the southern range of I. scapularis, resulting in a northward shift, rather than overall expansion, of the geographic range of climatic suitability for I. scapularis. Here we confined our study to Canada and the main regions of Lyme disease risk in the United States north of 40°N. Impacts of rainfall on off-host tick survival and on host-seeking activity are considered accounted for in the model in assuming a) tick populations only become established in woodlands where the microclimate is suitable for tick survival, and b) most temperate woodlands types occur where rainfall is sufficient for I. scapularis survival, which is supported by studies in Canada. Future projections for increased precipitation across much of Canada with climate warming are already being seen to occur, so rainfall changes are not expected to limit the northward spread of I. scapularis. For the present study, we modified the simulation model of in order to calculate R0 by the next generation operator as described by. Apart from the temperature values used to calculate tick development and host-finding rates, the values for host numbers (20 deer and 200 rodents) and all other parameter values were those used as starting values by. R0 was estimated using mean monthly temperature data for each year and location as described in the following sections. Variations in host abundance affect the final size of the tick population but not the temperature threshold. The temperature threshold would be affected by variations in the mortality rates of ticks in the environment, and slight variations in this have been observed in the field. For the sensitivity analysis of R0 to variations in model parameter values, see Supplemental Material, “Model sensitivity analysis,” pp. 7–11, Table S2 and Figures S5 and S6.

7Modeling R0 under current climate. For observed temperatures, we used Australian National University Splines (ANUSPLIN) of 10-km gridded daily time-series data, which were obtained by thin-plate smoothing spline interpolation of daily climate station observations while accounting for latitude, longitude, and elevation. ANUSPLIN data cover the 40 years (1971–2010) that encompass the period of Lyme disease emergence in North America, have coverage across northern North America, and account for missing data by temporal and spatial interpolation. A mean daily near-surface temperature was assigned to the 32 study sites, which were weather stations located over a wide range of the orographic and forest ecosystems of Ontario and Quebec or interpolations of ANUSPLIN data for locations in the United States. Monthly mean near-surface air temperatures were used to parameterize the I. scapularis population model for estimating annual values of R0 for each site, for each year from 1971 to 2010 (see Supplemental Material, “Variation in temperature and R0 amongst sites,” pp. 2–6). Values for annual cumulative degree days > 0°C (DD > 0°C), contemporaneous for each estimated R0 value, were computed as the accumulation of daily temperature > 0°C for each year for each site. The tick model calculated R0 for each year using temperature data for that year, but in reality R0 will depend on the temperature conditions over the 2- to 3-year life cycle of the tick. Therefore, we used the moving average of R0 over 3 years (that year, the previous year, and the subsequent year) to describe R0 for each year for each site; so for each site, we obtained a time series for R0 and DD > 0°C for 38 years (1972–2009) and DD > 0°C for 40 years (1971–2010).

8Modeling R0 using projected climate data. An ensemble of modeled temperature data available from three regional climate model (RCM) and two global climate model (GCM) runs were used to estimate future changes in R0 (see Supplemental Material, “Validation of climate model output,” pp. 11–14; Table S3). We extracted time series of daily temperature for the 30 Canadian sites from each model. Simulated temperature at a given site was defined as the mean of the closest grid values to that site, increasing confidence in the physical representativeness.

9We chose bias-corrected output from Canadian RCM CRCM4.2.3, version 4 to provide temperature data for the I. scapularis model because it relatively accurately and conservatively hindcasted observed ANUSPLIN data in comparison with other climate models and because it provided projected temperatures at a local scale. For full details justifying the climate model selection, see Supplemental Material, “Validation of climate model output,”pp. 11–14, Figures S7 and S8. Like other RCMs, CRCM4.2.3 dynamically downscales output from a coarser resolution GCM and produces data at a horizontal resolution of approximately 50 km. CRCM4.2.3 is driven by initial and boundary conditions of the Canadian GCM CGCM3.1 T47 (McFarlane et al. 2005; ). Hindcasting up to 2000 used greenhouse gas emissions for CGCM3.1 T47 as in the Coupled Model Intercomparison Project (CMIP) 20th century experiment. For future projections starting in 2001, we chose the A2 scenario (mid-high Green-House-Gas emission scenario) of the Intergovernmental Panel on Climate Change Special Report on Emission Scenarios because of the availability of regional climate model output using this scenario and because current actual trajectory of emissions corresponded best to this emissions scenario.

10Mapping R0 under current and future climate. We mapped 30-year mean values of R0 for I. scapularis using observed and projected values. Maps of DD > 0°C for North America north of 40°N and east of the Rocky Mountains were generated using observed data for 1971–2000, and using DD > 0°C projected by bias-corrected CRCM4.2.3 output for the period 2001–2070. We then computed R0 for each year from the gridded DD > 0°C data using the formula R0 = 1.072 × 10–6 DD > 0°C–4.658 × 10–3 DD > 0°C2 + 5.556, obtained using observed temperature data as described in Supplemental Material, “Mapping R0,”pp. 15–16, Figure S9 (the threshold for R0 ≥ 1 was DD > 0°C ≥ 2859.6°C). We then mapped mean 30-year values for R0 for the periods 1971–2000, 2011–2040, and 2041–2070 (Figure 2). Regions west of the Rocky Mountains were masked because it was assumed that I. scapularis will not cross the Rocky Mountains, west of which Lyme disease risk will continue to depend on transmission of Borrelia burgdorferi by the tick I. pacificus.

11Literature search on R0 ranges for parasite and pathogen systems. There are, to our knowledge, no equivalent estimates of how environmental changes may affect R0 of vectors. However, to better comprehend the ecological or epidemiological importance of projected changes of R0 of I. scapularis, we performed a literature review to obtain published estimates of how R0 for parasites and pathogens varies due to changes in factors already recognized as having ecological or epidemiological importance. These include variations in geographic location, host density, strain or genotype, disease control effort, and variations among different epidemics. Articles were searched in the National Centre for Biotechnology Information PubMed search site (http://www.ncbi.nlm.nih.gov/pubmed/) using search terms a) “basic reproduction number,” without specifying pathogens or parasites, and then b) “basic reproductive number,” repeated with one of the following terms: tick, mosquito, chagas, malaria, dengue, nematode, “seasonal influenza,” “pandemic influenza”, pH1N1, “avian influenza,” measles, HIV (human immunodeficiency virus), and fluke. Abstracts were reviewed, and relevant articles were reviewed in full. Relevant articles were those in which R0 for parasites and pathogens was calculated to explicitly estimate its value under field, rather than theoretical, conditions. This meant articles that employed simulation models using field data, fitting of epidemiological data (e.g., age–seroprevalence or age–infection prevalence), or other methods such as estimates from phylogenetic analysis. We did not review and use R0 ranges obtained in model-based sensitivity analyses, variations in R0 associated with seasonal variations in mosquito abundance in one location (which may vary from zero to very high values), estimates where control methods effectively eradicated disease resulting in almost infinite values for changes in R0, or model-predicted variations across whole potential geographic ranges that range from a theoretical high to zero values (e.g., ). We also did not use some very high modeled ranges of R0 for malaria [e.g., 1–11,000] when the modeling of empirical age–infection prevalence data produced strongly contrasting single-digit estimates of R0. The goal of the literature search was to provide an illustration of how important projected changes in R0 of I. scapularis could be, compared with R0 of other parasites and pathogens; it was not intended to be an exhaustive cataloging of all literature in this field. Further, we recognized that R0 estimates are not precise and vary according to the estimation method used.

12Using observed (ANUSPLIN) temperature data, R0 for I. scapularis in the late 1970s—when Lyme disease emerged in the northeastern United States—was estimated at approximately 3 and 1.9 in Old Lyme and Fort McCoy, respectively; at between 2 and 3 in Southern Ontario; approximately 1.5 in Huron Ontario and South-Western Quebec, but mostly < 1 in Upper Southern Ontario and the Boreal region (Figure 1). In Old Lyme, R0 increased almost linearly to approximately 3.5 by 1999 during the first period of expansion of I. scapularis in the northeastern United States. In Fort McCoy, R0 increased slightly, but this increase was small compared with interannual variations. In Southern Ontario, R0 increased to 4 by the early 2000s; during which time I. scapularis populations emerged at a number of locations in this region (Point Pelee National Park, Turkey Point, Rondeau Provincial Park; Figure 1A). In Huron Ontario and South-Western Quebec, R0 increased from 1.5 to 2.5 by the early 2000s; and subsequent to this (mostly from 2000 onward), I. scapularis populations began to emerge in South-Western Quebec (Figure 1D). In Upper Southern Ontario, R0 increased to > 1 in the late 1990s, but in the Boreal region R0 remained below unity for the whole 1971–2010 period (Figure 1).

13R0 values for I. scapularis obtained in model simulations using projected climate data were similar for an ensemble of climate models, and we used bias-corrected output from the regional climate model CRCM4.2.3 as a representative of the ensemble because of its spatial resolution and predictive accuracy. R0 for I. scapularis in Canada was projected to increase 1.5 to 2.3 times from the period 1971–2000 to 2001–2050, and 2.2 to 4.6 times from the period 1971–2000 to 2051–2069 depending on location (Figure 1, Table 1); and in the United States, R0 was predicted to approximately double to 7.1 and 5.2 in Old Lyme (Figure 1) and Fort McCoy, respectively, by 2051–2069 (Table 1). Increases in R0 to values > 1 predicted in regions where R0 was < 1 during the period 1970–2000 would be expected to facilitate range expansion of I. scapularis northward and possibly westward (Figure 2).

14The projected increases in R0 are equivalent, for the most part, to ranges of values of R0 estimated for other globally important parasites and pathogens associated with variations in major determinants of their ecology and epidemiology such as geographic location, pathogen genotype, different epidemics, reservoir host or vector density, and control efforts (Table 1).

15These findings suggest that increasing temperatures in northern North America that support an R0 for I. scapularis of > 1.5 have been coincident with, or in advance of, but not subsequent to, expanding numbers of locations where I. scapularis populations have become established. In Canada, where we have tracked the spread of I. scapularis, temperature has remained a statistically significant determinant of I. scapularis occurrence in field studies and analyses of surveillance data that accounted for alternative environmental determinants (e.g., host abundance, altitude, rainfall, habitat types, tick immigration rates). These observations supported a key role for temperature in I. scapularis populations becoming established at the northern edge of the tick’s range. Also, I. scapularis population expansion in Canada is occurring despite an overall deforestation rather than the reforestation thought to have driven the initial reemergence of Lyme disease in the United States. Forest fragmentation may enhance Lyme disease risk for a variety of reasons; however, I. scapularis ticks are invading Canada where forest fragmentation occurred over time scales long predating current I. scapularis invasion. Together, these findings suggest that even if recent warming in the region (5–10% increases in DD > 0°C; Figure 1) was not associated with global warming, a future warming climate will increase R0 of I. scapularis in northern North America. Increases in R0 may drive increased Lyme disease risk where it is already endemic (within limits determined by density-dependent regulation of the tick) and drive range expansion into more northern regions where it is currently absent. Furthermore, they support the hypothesis that climate warming in northeastern North America may have co-driven the emergence of Lyme disease risk, alongside other hypothesized factors such as reforestation and burgeoning deer populations, by facilitating the spread of I. scapularis from refuges. Expansion of I. scapularis in the northern United States has then provided the source of ticks to fuel its northward expansion into Canada.

16The immediate importance of future increased R0 of I. scapularis in the northeast and upper Midwest of North America is that a) regions currently climatically unsuitable become suitable for I. scapularis establishment (i.e., R0 changes from < 1 to > 1; Figure 2); b) in regions currently suitable for I. scapularis (where R0 > 1), tick invasion speed will accelerate as the likelihood of stochastic tick population fade-out reduces, and tick-borne pathogen invasion speed increases due to increasing tick abundance; and c) risk from I. scapularis–transmitted pathogens may increase where the tick and pathogen are already established due to increased tick abundance up to a point at which this is limited by density-dependent regulation. To date, I. scapularis invasion in the northern United States and Canada has been followed by invasion of the agent of Lyme disease, B. burgdorferi sensu stricto, hence we assume northward I. scapularis range expansion is synonymous with expansion in Lyme disease risk.

17The magnitude of projected increases in R0 of I. scapularis in the present study is of importance for the ecology and epidemiology of vector-borne diseases in general. This was illustrated by R0 ranges estimated for other globally important parasites and pathogens associated with variations in known key determinants of their ecology and epidemiology (Table 1). These ranges were of a similar magnitude to projected increases in R0 of I. scapularis with climate change. Tick species are likely to respond to the central tendency of increasing temperatures due to a) the long periods of inter-stadial development that take place in the surface layers of the soil where fluctuating air temperature are buffered, b) latency in responses of development rates to temperature changes minimizing effects of very short-term temperature fluctuations, c) their ability to return to soil-level refugia during extremes of heat, cold, drought, or rainfall while host seeking, d) their associations with woodland habitats within which a microclimate is buffered from extremes of temperature occurring in treeless areas, and e) ticks have no nonparasitic immature feeding stages whose survival is susceptible to short-term changes in weather as do dipteran vectors such as mosquitoes. Therefore, the increases in R0 projected here represent a possible magnitude of increase in mean R0 values arthropod vectors may experience with climate change. Around this mean, dipteran population R0 may fluctuate seasonally and annually over a much wider range due to the rapid effects of rainfall and temperature on reproduction and mortality rates. Abundance and geographic distributions of many mosquito-borne diseases are currently driven primarily by control efforts that superimpose on any climate effects. Large increases in vector R0 may, however, render current vector and vector-borne disease control methods ineffective as vector multiplication outstrips control efforts.

18R0 increases associated with climate change may be limited in some circumstances. Host population densities and habitat do not seem to be currently limiting on I. scapularis range expansion, but they may be in the future. Mosquitoes and other dipteran vectors can be dispersed by wind, but ticks need hosts for their dispersal to effect range expansion. I. scapularis are dispersed over long distances by migratory birds, but ticks that are not carried by migratory animals would be expected to have less capacity to invade climatically suitable environments. Ticks of public health importance such as I. ricinus and I. scapularis are mostly host and woodland habitat generalists, which facilitates range changes, whereas the more highly specialized the niche of a species, the less likely it will be to be dispersed and/or capable of becoming established in new locations.

19We estimated the degree to which projected climate change may impact the ecology of arthropod vectors and, by inference, vector-borne diseases. The emergence of Lyme disease in North America may itself have been partly driven by recent climate change. Confidence in our projections is increased by observed changes in temperature and estimated R0 for the vector associated with the actual emergence of the vector and the vector-borne diseases it transmits. Our findings suggest that effort should be refocused on assessing the health risks due to vector-borne disease, particularly vector-borne zoonoses, associated with our changing climate.