Ticks and Tick-borne Diseases 6 (2016) 615–622 Contents lists available at ScienceDirect Ticks and Tick-borne Diseases journal homepage: www.elsevier.com/locate/ttbdis Original article Climate change influences on the annual onset of Lyme disease in the United States Andrew J. Monaghan a,∗ , Sean M. Moore b , Kevin M. Sampson a , Charles B. Beard c , Rebecca J. Eisen c a Research Applications Laboratory, National Center for Atmospheric Research, 3090 Center Green Dr., Boulder, CO 80301, USA Johns Hopkins School of Public Health, 615 N. Wolfe St., Baltimore, MD 21205, USA c Division of Vector-Borne Diseases, Centers for Disease Control and Prevention, 3150 Rampart Rd., Fort Collins, CO 80522, USA b a r t i c l e i n f o Article history: Received 2 October 2014 Received in revised form 23 April 2015 Accepted 7 May 2015 Available online 15 May 2015 Keywords: Lyme disease Climate change Ixodes scapularis Borrelia burgdorferi a b s t r a c t Lyme disease is the most commonly reported vector-borne illness in the United States. Lyme disease occurrence is highly seasonal and the annual springtime onset of cases is modulated by meteorological conditions in preceding months. A meteorological-based empirical model for Lyme disease onset week in the United States is driven with downscaled simulations from five global climate models and four greenhouse gas emissions scenarios to project the impacts of 21st century climate change on the annual onset week of Lyme disease. Projections are made individually and collectively for the 12 eastern States where >90% of cases occur. The national average annual onset week of Lyme disease is projected to become 0.4–0.5 weeks earlier for 2025–2040 (p < 0.05), and 0.7–1.9 weeks earlier for 2065–2080 (p < 0.01), with the largest shifts for scenarios with the highest greenhouse gas emissions. The more southerly midAtlantic States exhibit larger shifts (1.0–3.5 weeks) compared to the Northeastern and upper Midwestern States (0.2–2.3 weeks) by 2065–2080. Winter and spring temperature increases primarily cause the earlier onset. Greater spring precipitation and changes in humidity partially counteract the temperature effects. The model does not account for the possibility that abrupt shifts in the life cycle of Ixodes scapularis, the primary vector of the Lyme disease spirochete Borrelia burgdorferi in the eastern United States, may alter the disease transmission cycle in unforeseen ways. The results suggest 21st century climate change will make environmental conditions suitable for earlier annual onset of Lyme disease cases in the United States with possible implications for the timing of public health interventions. © 2016 The Authors. Published by Elsevier GmbH. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Introduction Lyme disease is a multisystem tick-borne bacterial zoonosis that is endemic in parts of North America, Europe and Asia. In the United States, Lyme disease is the most commonly reported vector-borne illness (CDC, 2008), with more than 25,000 Lyme disease cases reported annually since 2007 (CDC, 2014). The majority of Lyme disease cases are reported from Northeastern and north-central States where nymphal Ixodes scapularis ticks serve as the primary bridging vectors of the pathogenic bacterium Borrelia burgdorferi sensu stricto from zoonotic hosts to humans (CDC, 2008; Piesman, ∗ Corresponding author at: National Center for Atmospheric Research, P.O. Box 3000, Boulder, CO 80307-3000, USA. Tel.: +1 303 497 8424. E-mail addresses: monaghan@ucar.edu (A.J. Monaghan), semoore@jhsph.edu (S.M. Moore), ksampson@ucar.edu (K.M. Sampson), cbb0@cdc.gov (C.B. Beard), dyn2@cdc.gov (R.J. Eisen). 1989). Lyme disease transmission occurs seasonally, and the majority of human cases report onset of clinical signs of infection during the months of June, July and August, a period that corresponds with exposure to the nymphal life stage of I. scapularis (CDC, 2008; Piesman, 1989). The geographic distribution of Lyme disease is focal, and inter-annual variation in case counts and seasonal onset is considerable (Diuk-Wasser et al., 2012; Moore et al., 2014). Because Lyme disease cases can only occur in areas where humans encounter B. burgdorferi-infected ticks, much of the variability in where and when Lyme disease cases occur is attributable to the geographic distribution and seasonal host-seeking patterns of the ticks that serve as vectors of B. burgdorferi. Although at local scales host community structure plays a large role in determining the density of infected nymphs (Mather et al., 1989; Ostfeld et al., 2006), at regional scales, temperature, humidity and precipitation are robust predictors of spatial and temporal distributions of I. scapularis (Brownstein et al., 2003; Diuk-Wasser et al., 2006, 2010; Estrada-Pena, 2002). These variables have also been associated with http://dx.doi.org/10.1016/j.ttbdis.2015.05.005 1877-959X/© 2016 The Authors. Published by Elsevier GmbH. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 616 A.J. Monaghan et al. / Ticks and Tick-borne Diseases 6 (2016) 615–622 the geographical and temporal distributions of human cases of Lyme disease in the United States (Ashley and Meentemeyer, 2004; McCabe and Bunnell, 2004; Moore et al., 2014; Ogden et al., 2014; Subak, 2003; Tran and Waller, 2013). Understanding how meteorology impacts the seasonality of Lyme disease case occurrence can aid in targeting limited prevention resources and may shed light on how climate change could affect the seasonal occurrence of the disease (Gray, 2008). Based on cases reported through the National Notifiable Diseases Surveillance System (NNDSS) from 1992 to 2007, a companion study modeled the timing of the start, peak, duration and end of the Lyme disease season for 12 endemic States in the Northeast, mid-Atlantic and upper Midwest as a function of meteorological variables (Moore et al., 2014). Moore et al. (2014) found significant associations between meteorological variables and the timing of the onset, peak and duration of the Lyme disease season; however, meteorological variables did not predict the end of the season. The strongest associations were found for the onset of the Lyme disease season. Across all States and years, the beginning of the Lyme disease season ranged from week 16–26 of the calendar year, and 60% of the variation was attributable to the geographic and temporal variability of climatic and other environmental factors. The Lyme disease season began earlier in more southerly and coastal States compared with more northerly and inland States. Earlier onset of the Lyme disease season was positively associated with warmer and more humid conditions and lower rainfall amounts during the preceding winter and spring months. Other studies also indicate that warmer and/or more humid conditions are associated with I. scapularis characteristics including geographic distribution (Brownstein et al., 2003; Estrada-Pena, 2002; Ogden et al., 2008) and increased density of host-seeking nymphal ticks (Diuk-Wasser et al., 2006, 2010). Given that climate models project a temperature increase over the United States of 1.5–5.5 ◦ C by the end of the 21st century following a 0.8 ◦ C increase during the 20th century, and that rainfall amounts will likely continue rising over the Northeastern U.S. (USGCRP, 2014) where most Lyme disease cases occur, it is plausible that climate change may affect the annual onset of Lyme disease in forthcoming decades. While previous studies have primarily examined climate change impacts on the geographic distribution, host-seeking phenology and reproductive rate of Ixodes scapularis (e.g., Brownstein et al., 2005; Ogden et al., 2006, 2008, 2014; Simon et al., 2014; Levi et al., 2015), none have investigated climate change impacts on the seasonality of human Lyme disease cases in the United States. Here, the national model of Moore et al. (2014) is employed to investigate how projected 21st century climate changes may affect the timing of annual Lyme disease onset in the eastern United States. Development and implementation of such models can aid in determining the magnitude by which climate change may drive shifts in the annual onset of Lyme disease cases, allowing public health officials to gauge whether it will be necessary to adjust future interventions to account for altered seasonality of Lyme disease. Materials and methods National Lyme disease model The best-fit (adjusted R2 = 0.785) national-level model for Lyme disease onset presented in Moore et al. (2014) is: LOW = 17.56252 − 0.014 × GDDW20 + 0.945 × SDM5 + 0.009 × PRCPAW8 + 0.093 × DIST (1) where LOW is Lyme Onset Week (week 1 is defined at the beginning of the calendar year), GDDW20 is the cumulative growing degree days from week 1 to week 20, SDM5 is the mean saturation deficit in mmHg in the 5 weeks before the onset week, PRCPAW8 is the cumulative rainfall in mm from week 8 (approximately the beginning of spring) through the onset week, and DIST is distance in decimal degrees to the Atlantic Ocean coastline from the weighted mean center of each State’s total Lyme disease cases. LOW is defined as the week with the maximum percent increase in the number of Lyme disease cases over the previous week. The model indicates that Lyme disease season is expected to begin 1.4 weeks earlier for each additional 100 cumulative GDDs through week 20, about 1 week later for each 1 mmHg increase in saturation deficit (i.e., if humidity decreases with respect to the air temperature), and about 0.9 weeks later for each 100 mm increase in cumulative precipitation between week 8 and the beginning of the Lyme disease season. The time-invariant variable DIST provides a measure of the maritime or continental climate influences in a State. Compared to inland areas, near-coastal areas often have smaller climatic fluctuations due to the moderating influence of the ocean (Bailey, 1980). The model is applied to each State and year separately and then results are aggregated to the regional or national level as needed, or temporally averaged to obtain long-term averages of LOW. The LOW model was developed using human cases of Lyme disease reported to the Centers for Disease Control and Prevention (CDC) by State and territorial health departments as part of the National Notifiable Disease Surveillance System (NNDSS) from 1992 to 2007 (CDC, 2009). Over 95% of Lyme diseases cases in the United States occurred in 13 States in the east and north-central regions during the study period: Connecticut, Delaware, Maine, Maryland, Massachusetts, Minnesota, New Hampshire, New Jersey, New York, Pennsylvania, Rhode Island, Virginia, and Wisconsin. Reports from Delaware during the study period did not include an illness onset date and Delaware was subsequently excluded from the analysis. Therefore, the national model was developed using Lyme disease case data from the 12 States accounting for >90% of all United States cases reported for 1992–2007. Additional details on case data, model development, and the methodology for defining observed LOW are in Moore et al. (2014). Climate data Moore et al. (2014) describe in detail the historical climate data used in the development of the national model. The data are briefly summarized here for clarity. Historical meteorological variables were obtained or derived from the 1/8th degree primary forcing data for Phase 2 of the North American Land Data Assimilation System (NLDAS-2) (Xia et al., 2012). The observation-constrained meteorological variables of NLDAS-2 span 1979-present and are considered to be of suitable quality for use in climate-sensitive human health applications over North America (Luber, 2014). The NLDAS-2 variables were aggregated to the county-level using the Zonal Statistics spatial analysis tool in ArcGIS (Esri, Redlands, CA). State averages of the NLDAS-2 variables were then calculated annually for 1992–2007 using the county-level data, weighted by the percentage of cases in each county during 1992–2007. Future climate projections were selected from a multi-model ensemble of atmosphere-ocean global climate models (AOGCMs) that participated in phase five of the Coupled Model Intercomparison Experiment (CMIP5) (Taylor et al., 2012). The CMIP5 simulations support the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC, 2013) and the Third National Climate Assessment for the United States (USGCRP, 2014). Specifically used were AOGCM simulations from a database of CMIP5 climate and hydrology projections that have been empirically downscaled with the bias-corrected spatial disaggregation method (Archive Collaborators, 2014; Brekke et al., 2013; Maurer et al., 2007). The empirically downscaled projections were chosen A.J. Monaghan et al. / Ticks and Tick-borne Diseases 6 (2016) 615–622 617 Table 1 List of climate modeling centers and AOGCMs used in this study. Modeling center (or group) Institute ID Model name Community Earth System Model Contributors NOAA Geophysical Fluid Dynamics Laboratory NASA Goddard Institute for Space Studies Met Office Hadley Centre (additional HadGEM2-ES realizations contributed by Instituto Nacional de Pesquisas Espaciais) Atmosphere and Ocean Research Institute (The University of Tokyo), National Institute for Environmental Studies, and Japan Agency for Marine-Earth Science and Technology NSF-DOE-NCAR NOAA GFDL NASA GISS MOHC/INPE CESM1 (CAM5) GFDL-CM3 GISS-E2-R HadGEM2-ES MIROC MIROC5 because they are mapped on the same 1/8th degree domain as the historical NLDAS-2 data used in the original LOW model development, and because the database has a variety of AOGCMs and scenarios available, which facilitate an uncertainty analysis. Five AOGCMs were selected from the database (Table 1) according to the following three criteria: 1. They have at least one simulation available for all four CMIP5 climate change scenarios (described below); 2. They rank among the top of the CMIP5 AOGCMs in their ability to simulate observed temperature and rainfall globally according to Knutti et al. (2013); 3. They each come from a different model genealogy according to Knutti et al. (2013), ensuring each model is sufficiently unique from the others. t-test for the difference between the means of the histograms. The statistical significance tests address two types of uncertainty: (1) that due to the choice of global climate model (the uncertainty among the 5 different AOGCMs), and (2) that due to the interannual variability of the climatic variables (the uncertainty among each year of the baseline and future 16-year periods). A third type of uncertainty is not addressed in this paper: that due to the uncertainty of the statistical model that describes LOW (Eq. (1)). The LOW model uncertainty could hypothetically be addressed by applying several different plausible national LOW models, however the four leading national LOW models published in Moore et al. (2014) have nearly identical explanatory variables, with only slight deviations from one another. Therefore, only the ‘best’ national LOW model is employed here in order to avoid adding unnecessary complexity. Results AOGCM simulations from all four future emissions scenarios from CMIP5 are used. These are known as representative concentration pathway (RCP) scenarios (Moss et al., 2010) and include RCP2.6, RCP4.5, RCP6.0 and RCP8.5, with the numbers indicating the greenhouse-gas (GHG) radiative forcing near 2100 (e.g., 2.6 W m−2 , 4.5 W m−2 , 6.0 W m−2 and 8.5 W m−2 ). RCP2.6 is a low emissions scenario with aggressive reductions in GHG emissions representing a technically feasible trajectory for limiting the global mean temperature increase to 2 ◦ C or less (van Vuuren et al., 2011). RCP4.5 is a low-to-moderate emissions scenario representing a trajectory that may be plausible if, for instance, GHG emissions pricing were introduced in order to limit and stabilize radiative forcing (Thomson et al., 2011). RCP6.0 is a moderate GHG emissions scenario that is similar to RCP4.5 in that a variety of strategies for reducing GHGs would be applied to eventually stabilize radiative forcing near the end of the 21st century (Masui et al., 2011). RCP8.5 is a high GHG emissions scenario representing a plausible trajectory if little is done to curb greenhouse gas emissions (Riahi et al., 2011). For the first ensemble member of each AOGCM and RCP scenario combination, average monthly maximum temperature (TMAX ), minimum temperature (TMIN ), and precipitation (PRCP) data for the 12 State region were aggregated to the county level from their 1/8th degree grid, then averaged to the State level with the countybased weighting technique described above. This was done for the 16-year baseline period over which the LOW model was developed (1992–2007), as well as for two future 16-year periods: 2025–2040 and 2065–2080. The TMAX , TMIN and PRCP fields were then used to calculate GDDW20 , SDM5 , and PRCPAW8 for the two future periods via a delta-based method (e.g., Hay et al., 2000), as described in the Supplemental Material, Part 1. Statistical analysis Present-day and future histograms of annual LOW and climatic data are constructed for numerous combinations of AOGCMs, RCP scenarios, and regions (State and national levels). The frequency distributions are approximately normal and therefore, future changes of LOW or climatic variables for each case are tested for significance (p < 0.01 and p < 0.05) with a two-tailed Student’s National-level results for historical and future distributions of LOW are shown in Fig. 1. The multi-model future projection of Fig. 1. Box plots comparing the distributions of the national-level historical data for annual Lyme Onset Week (LOW) with the AOGCM multi-model mean distributions of LOW for each of the four RCP scenarios and two future periods. Each box plot shows the values of LOW for the maximum (top of dashed line), 75th percentile (top of box), mean (line through middle of box), 25th percentile (bottom of box) and minimum (bottom of dashed line) of the distribution. All distributions are comprised of values for 12 States and 16 years (N = 192). Circles along the top of each panel indicate whether the AOGCM multi-model mean is significantly different from the historical mean (see top legend). Box plot colors indicate different time periods (see middle legend). Black symbols on each box plot indicate the mean value of LOW from each individual AOGCM that contributes to the multi-model ensemble (see bottom legend). 618 A.J. Monaghan et al. / Ticks and Tick-borne Diseases 6 (2016) 615–622 national-average LOW for both the early (2025–2040) and late (2065–2080) 21st century periods is significantly earlier for all four RCP scenarios compared to the historical national average LOW for 1992–2007 of 21.2 weeks. The early 21st century changes are similar among the four RCP scenarios because their respective GHG emissions trajectories do not diverge substantially until after midcentury. On average, LOW is projected to become 0.4–0.5 weeks earlier for 2025–2040, and 0.7–1.9 weeks earlier for 2065–2080, depending on the scenario. The strongest changes for the late century scenario are for the highest GHG emissions scenario, RCP8.5. National-level results for the historical and future distributions of the climatic variables associated with LOW are shown in Fig. 2. The average temperatures for Jan-May (TJAN-MAY ), though not directly used in Eq. (1), give a sense for the temperature increases that occur during the winter and spring months leading up to LOW. National-level mean TJAN-MAY increases by 1.2–1.7 ◦ C for 2025–2040, and by 1.8–4.5 ◦ C for 2065–2080, depending on the scenario. Consistent with the warmer temperatures, GDDW20 increases by 54–76 GDDs for 2025–2040, and by 99–232 GDDs for 2065–2080. SDM5 increases by 0.25–0.34 mmHg for 2025–2040, and by 0.43–0.92 mmHg for 2065–2080 as a result of warmer temperatures under constant relative humidity conditions (see Supplemental Materials, Part 1). PRCPAW8 increases by 18–32 mm for 2025–2040, and by 30–53 mm for 2065–2080. The changes for all four climate variables in Fig. 2 are statistically significant for both future periods and all RCP scenarios. For all climatic variables, the smallest changes by the late 21st century are for the RCP2.6 scenario and the largest are for RCP8.5. The historical values of LOW and associated climate variables for 1992–2007 for each State and nationally are shown in Table 2. The differences from these historical values for 2065–2080 for the “best-case” RCP2.6 and “worst-case” RCP8.5 scenarios are shown in Tables 3 and 4 respectively (see Supplemental Materials, Part 2 for results from other RCP scenarios and for 2025–2040). The States are categorized into the four regions following Moore et al. (2014) (note that the “mid-Atlantic” region is identical to the “south” region in Moore et al., 2014). Some 1992–2007 State-level values of LOW are up to a few tenths of a week different than those presented in Moore et al. (2014) because missing values were present in the Moore et al. (2014) data, whereas Eq. (1) is used in this study to calculate LOW from the climatic data, so no missing values exist. This approach provides a more representative average of LOW for the 1992–2007 historical period. By 2065–2080, LOW is projected to become significantly earlier (p < 0.05) for 4-of-12 States for the RCP2.6 scenario, and for 11-of-12 States for RCP8.5. Maine undergoes the smallest changes, −0.2 weeks (RCP2.6) and −0.9 weeks (RCP8.5), neither being statistically significant at the 95% confidence interval. Virginia undergoes the largest changes, −1.5 weeks (RCP2.6) and −3.5 weeks (RCP8.5), both statistically significant. In general, the largest future changes in LOW occur in the comparatively warmer midAtlantic States where LOW is historically earliest and the increase in GDDW20 is largest on an absolute basis, followed by the Midwest which, despite having the coolest Jan-May temperatures, exhibits the strongest warming ( TJAN-MAY ) and the largest percentage increase in GDDW20 . LOW increases are smaller in the North and East where there are comparatively smaller increases in TJAN-MAY and GDDW20 , and larger increases in SDM5 and PRCPAW8 (both of which are favorable for later LOW), though it is noteworthy that PRCPAW8 changes are not statistically significant for any State for RCP2.6, and only for 6–12 States for RCP8.5. Changes in GDDW20 contribute most to LOW, about 3–5 times more than SDM5 and about 3–23 times more than PRCPAW8 (comparison of final three columns in Tables 3 and 4). Additional commentary on the Statelevel results is included in Section “Discussion”. Discussion Future projections based on five AOGCMs and four emissions scenarios suggest an earlier beginning to the Lyme disease season nationally by 0.4–0.5 weeks (2025–2040) and 0.7–1.9 weeks (2065–2080). The greatest changes were observed under the highest GHG scenario (RCP8.5). Notably, regional differences in LOW are expected. Larger changes in LOW are projected for the more southerly States of the mid-Atlantic region compared to the more northerly States of the North and Midwest. For example, for the RCP8.5 2065–2080 case LOW becomes 3.5 weeks earlier in Virginia compared to 0.9 weeks in Maine and 1.8 weeks in Minnesota, despite smaller increases in average winter-spring temperatures (TJAN-MAY ) in Virginia (4.0 ◦ C) compared to Maine (4.9 ◦ C) and Minnesota (5.4 ◦ C). This raises the question of whether the LOW model is overly sensitive to the choice of a threshold-based variable, GDDW20 , as a predictor. To explain, even though the GDDs increase more on a percentage basis in the northern States, the absolute increases are generally smaller because the base GDD threshold temperature of 10 ◦ C is exceeded for a shorter period of time in cooler areas. Are the differential changes in LOW between northern and midAtlantic States a model artifact resulting from absolute changes in GDDs in the cooler northern States being smaller than for the mid-Atlantic States? To address the differential changes in LOW between northerly and southerly States, observed historical interannual variability of LOW for 1992–2007 in the mid-Atlantic States was assessed to determine if it is disproportionately larger than in the northern States in comparison to the inter-annual variability in winter and spring temperatures. The two States with the coldest and warmest average Jan-May temperatures, Minnesota and Virginia, were compared. The standard deviation of TJAN-MAY for Minnesota is 1.95 ◦ C compared to 0.94 ◦ C for Virginia, therefore Minnesota exhibits more than double the inter-annual variability for winter and spring temperatures compared to Virginia. However, despite having much larger temperature variability the standard deviation of annual LOW is similar for Minnesota and Virginia, 1.4 weeks and 1.3 weeks respectively (note that these values differ slightly from the standard deviations for LOW presented in Table 2 because they are observed, whereas model results are presented in Table 2). Therefore, similar changes in temperature are associated with larger changes in LOW in Virginia compared to Minnesota, and in general for mid-Atlantic versus northern States (results are similar for other States, not shown). The reduced sensitivity of LOW to a given change in temperature in the northern States has two origins. First, GDDs increase less for a unit increase in temperature in the northern regions because they are based on a threshold temperature that is generally exceeded more often in the midAtlantic States. During model validation, the authors found the same north/south differences among GDDs even when lowering the GDD threshold to from 10 ◦ C to 6 ◦ C, given that previous studies suggest that a minimum temperature threshold required for I. scapularis or closely related I. ricinus to commence host-seeking activity following winter diapause may range from 6 to 10 ◦ C (Gray, 1984, 1985; MacLeod, 1935, 1936; Mount et al., 1997; Perret et al., 2000; Tagliapietra et al., 2011). Second, saturation deficit and precipitation contribute more toward offsetting the effects of GDDs on LOW in the northern States (Tables 3 and 4). In summary, the observed inter-annual variability of LOW and associated climatic drivers during the historical period are in agreement with the model results that show a differential sensitivity of LOW to temperature variations between northern and mid-Atlantic States. The results indicate that a threshold-based temperature proxy such as GDDs, which can account for differential sensitivity among States, is an appropriate explanatory variable. Additionally, the sensitivity analysis indicates that future changes for LOW may be larger in A.J. Monaghan et al. / Ticks and Tick-borne Diseases 6 (2016) 615–622 619 Fig. 2. As in Fig. 1, but for the national-level climate variables including Jan-May average temperature (A), cumulative GDDs through Week 20 (B), average saturation deficit in the 5 weeks before LOW (C), and cumulative precipitation from week 8 until LOW (D). Table 2 State- and national-level historical (1992–2007) mean ± standard deviation for LOW and associated climate variables. DIST is included for completeness. Region State LOW (weeks) TJAN-MAY (◦ C) GDDW20 (GDDs) SDM5 (mmHg) PRCPAW8 (mm) DIST (deg) Midwest MN WI 22.8 ± 0.8 21.9 ± 0.7 0.1 ± 1.8 0.5 ± 1.7 120 ± 51 106 ± 46 4.2 ± 0.5 3.1 ± 0.5 207 ± 61 205 ± 41 11.45 11.72 North ME MA NH 22.0 ± 0.9 22.0 ± 0.6 22.2 ± 0.8 1.8 ± 1.2 4.3 ± 1.1 2.6 ± 1.1 44 ± 27 75 ± 31 70 ± 36 2.3 ± 0.5 2.7 ± 0.4 2.8 ± 0.5 323 ± 96 322 ± 61 325 ± 83 0.18 0.31 0.60 East CT RI NJ NY PA 21.7 21.9 20.7 21.3 20.2 mid-Atl MD VA 19.2 ± 1.2 18.2 ± 1.0 8.3 ± 1.0 8.8 ± 1.0 224 ± 57 281 ± 59 National – 21.2 ± 0.8 4.4 ± 1.2 129 ± 42 ± ± ± ± ± 0.7 0.6 0.9 0.8 0.9 4.6 5.0 6.0 4.8 6.1 ± ± ± ± ± 1.1 1.0 1.1 1.1 1.2 105 89 154 105 176 ± ± ± ± ± 35 30 45 37 50 3.0 2.8 3.0 2.8 3.0 ± ± ± ± ± 0.5 0.4 0.5 0.4 0.4 309 330 276 274 252 ± ± ± ± ± 58 67 67 60 68 0.34 0.12 0.35 0.54 0.64 2.9 ± 0.5 3.0 ± 0.6 228 ± 77 189 ± 60 0.06 0.30 3.0 ± 0.5 270 ± 67 2.22 620 A.J. Monaghan et al. / Ticks and Tick-borne Diseases 6 (2016) 615–622 Table 3 The 2065–2080 minus 1992–2007 differences of RCP2.6 AOGCM multi-model mean values of LOW, and the differences and contributions of associated climate variables.a Region State Contribution of change to LOW Change: 2065–2080 minus 1992–2007 LOW (weeks) TJAN-MAY (◦ C) GDDW20 (GDDs) SDM5 (mmHg) PRCPAW8 (mm) of GDDW20 (weeks) of SDM5 (weeks) of PRCPAW8 (weeks) Midwest MN WI −0.8 −1.0 2.5 2.4 128 125 0.74 0.55 29 28 −1.8 −1.7 0.7 0.5 0.3 0.3 North ME MA NH −0.2 −0.3 −0.5 1.8 1.6 1.9 57 65 84 0.32 0.35 0.41 32 28 31 −0.8 −0.9 −1.2 0.3 0.3 0.4 0.3 0.2 0.3 East CT RI NJ NY PA −0.4 −0.3 −0.8 −0.5 −1.1 1.6 1.3 1.6 1.7 1.7 81 61 111 84 132 0.38 0.28 0.42 0.40 0.46 41 28 37 36 33 −1.1 −0.9 −1.5 −1.2 −1.8 0.4 0.3 0.4 0.4 0.4 0.4 0.3 0.3 0.3 0.3 mid-Atl MD VA −1.0 −1.5 1.5 1.6 109 149 0.39 0.46 22 18 −1.5 −2.1 0.4 0.4 0.2 0.2 National – −0.7 1.8 99 0.43 30 −1.4 0.4 0.3 a Future values with statistically significant change (p < 0.05) compared to 1992–2007 are underlined. the mid-Atlantic States compared to the northern States for similar increments of warming, especially given that combined increases in the offsetting saturation deficit and precipitation variables are projected to be larger in the northern States. As with all empirical models, there is always the possibility of the independent variables being serendipitously correlated with the dependent variables. However, numerous studies that have linked warmer and/or more humid conditions with I. scapularis characteristics in a manner conducive with earlier LOW support the associations between LOW and the explanatory climatic variables GDDW20 , and SDM5 and PRCPAW8 in equation (1) (Brownstein et al., 2003; Diuk-Wasser et al., 2006, 2010; Estrada-Pena, 2002). While the positive association between PRCPAW8 and LOW may seem counter-intuitive given that SDM5 is also positively associated with LOW – i.e., higher precipitation is associated with later LOW, whereas one might expect lower precipitation to be associated with later LOW since humidity is lower – it is important to note PRCPAW8 may be impacting different aspects of the transmission cycle than SDM5 , and at different timescales. PRCPAW8 covers a period of 2–4 months prior to LOW, whereas SDM5 covers a 7 day period 5 weeks prior to LOW. Greater PRCPAW8 is associated with cooler winter and spring temperatures as manifested by a negative correlation with the explanatory climatic variable GDDW20 (R = −0.35; p < 0.01); cooler temperatures associated with wetter conditions may delay the onset of Lyme disease by delaying interstadial tick development (Ogden et al., 2004), by delaying host-seeking activities (Eisen et al., 2002), or by delaying the springtime growth of vegetation that ticks exploit (Moore et al., 2014). Additionally, reduced human outdoor activity (and exposure to ticks) throughout the spring and early summer may also be associated with cooler, wetter conditions during that period (e.g., Fisman, 2007). In conclusion, the model projects that increasing temperatures during the 21st century are expected to result in an earlier onset of Lyme disease cases in the eastern United States. Notably, the focus is on the temporal occurrence of human infections, and therefore the model captures the complexity of how meteorological variables influence human behaviors resulting in contact with B. burgdorferi infected nymphal ticks, and host-seeking phenology of nymphal ticks. Although both human behavior and nymphal host seeking phenology are influenced by meteorological conditions (Randolph, 2004; Ogden et al., 2004), it is likely that they respond to different thresholds and at different time scales (Moore et al., 2014). As a result, the impacts of climate variability and change on each of these mechanisms individually cannot be elucidated within our modeling framework, and therefore the model can only estimate their Table 4 As in Table 3 but for the RCP8.5 scenario.a Region State Contribution of change to LOW Change: 2065–2080 minus 1992–2007 ◦ LOW (weeks) TJAN-MAY ( C) GDDW20 (GDDs) SDM5 (mmHg) PRCPAW8 (mm) of GDDW20 (weeks) of SDM5 (weeks) of PRCPAW8 (weeks) Midwest MN WI −1.8 −2.3 5.4 5.3 274 278 1.45 1.09 69 61 −3.8 −3.9 1.4 1.0 0.6 0.6 North ME MA NH −0.9 −1.1 −1.5 4.9 4.3 4.8 158 168 213 0.74 0.80 0.91 70 55 64 −2.2 −2.4 −3.0 0.7 0.8 0.9 0.6 0.5 0.6 East CT RI NJ NY PA −1.5 −1.3 −2.3 −1.5 −2.8 4.2 3.9 4.1 4.4 4.2 205 175 251 199 286 0.87 0.73 0.88 0.86 0.93 65 55 47 57 39 −2.9 −2.5 −3.5 −2.8 −4.0 0.8 0.7 0.8 0.8 0.9 0.6 0.5 0.4 0.5 0.4 mid-Atl MD VA −2.5 −3.5 3.9 4.0 250 326 0.83 0.93 25 22 −3.5 −4.6 0.8 0.9 0.2 0.2 National – −1.9 4.4 232 0.92 53 −3.2 0.9 0.5 a Future values with statistically significant change (p < 0.05) compared to 1992–2007 are underlined. A.J. Monaghan et al. / Ticks and Tick-borne Diseases 6 (2016) 615–622 combined impacts. In practice, how human behavior may change as a result of climate change may be relatively unpredictable. Although the LOW model provides a quantitative assessment of how the annual onset of Lyme disease cases may shift as climate changes during the 21st century, it does not explore how incidence changes in relation to meteorological variables and therefore does not address how an earlier onset could translate to changes in overall annual disease incidence. An earlier onset could result in more cases if increased temperatures result in greater abundance of ticks and increased contact rates with humans, or could result in the same number of cases or fewer cases if tick numbers and human-tick encounter rates remain the same, but tick activity shifts seasonally. Because the end of the Lyme disease transmission season could not be modeled (Moore et al., 2014), it is impossible to ascertain from this model whether an earlier LOW would translate to a longer season overall. It remains unclear if increased temperatures could condense the duration of the season as is observed in warmer Mediterranean climates (Eisen et al., 2003; Gray et al., 2008, 2009), or if the season might be extended into autumn and winter months due to increased nymphal and adult tick questing during these periods (Dautel et al., 2008; Gray et al., 2009). Alternatively, nymphal host seeking could shift abruptly from nymphs being active in the spring to the autumn. Such a situation was proposed to arise if increasingly warmer spring and summer conditions resulted in acceleration of the pre-oviposition and pre-eclosion periods, which would allow larvae to emerge and become active earlier in the summer. Additionally, warmer summer and autumn temperatures were hypothesized to facilitate faster larval development, resulting in increased numbers of larvae that were able to feed and engorge in summer and molt into nymphs that feed in late summer or autumn of the same calendar year (Ogden et al., 2008). The model does not account for the possibility that abrupt shifts in the life cycle of I. scapularis may alter the disease transmission cycle in unforeseen ways. The model results presented within should therefore be taken within the context of these and other potentially confounding factors, and should not be interpreted as more than what they are: a projection of enhanced climatic suitability for earlier LOW in the 21st century. Such knowledge may be useful for informing stakeholders of how the timing of Lyme disease prevention efforts may need to be altered, such as the application of acaracides in springtime to coincide with nymphal emergence (Rand et al., 2010), or public awareness campaigns to reduce human exposure to ticks (Connally et al., 2009; Hayes and Piesman, 2003). Conclusions A climate-based empirical model was driven with an ensemble of downscaled CMIP5 climate model simulations to project the impacts of climate change on the annual onset week of Lyme disease nationally and for 12 States. The average LOW is projected to become 0.4–0.5 weeks earlier nationally for 2025–2040 (p < 0.05), and 0.7–1.9 weeks earlier for 2065–2080 (p < 0.01). The smallest changes are projected for the lowest greenhouse gas emissions scenario, RCP2.6, while the largest changes would occur for the highest greenhouse gas emissions scenario (RCP8.5). Warming temperatures during the winter and spring months increase GDDW20 and cause earlier LOW in the future. Projected increases in SDM5 and PRCPAW8 both partially offset the effects of warming on LOW. Regionally, the mid-Atlantic States are projected to have larger changes in LOW compared to the northern States, a result the historical record supports. The results of the present study suggest that 21st century climate change, particularly increasing temperatures, will likely make environmental conditions suitable for earlier onset of Lyme disease cases in the United States. 621 Competing interest The authors declare they have no actual or potential competing financial interests. Acknowledgements We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for The Coupled Model Intercomparison Experiment, and we thank the climate modeling groups (listed in Table 1 of this paper) for producing and making their model output available. The U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led the development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. The Centers for Disease Control and Prevention funded this work. The National Science Foundation supports the National Center for Atmospheric Research. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.ttbdis.2015.05.005. References Archive Collaborators, 2014. Downscaled CMIP3 and CMIP5 Climate and Hydrology Projections, Available at: http://gdo-dc.ucllnl.org/downscaled cmip projections (accessed 17.06.14). Ashley, S.T., Meentemeyer, V., 2004. Climatic analysis of Lyme disease in the United States. Clim. Res. 27, 177–187. Bailey, R.G., 1980. Description of the Ecoregions of the United States (Misc Publ 1391). U.S. Department of Agriculture, Washington. Brekke, L., Thrasher, B.L., Maurer, E.P., Pruitt, T., 2013. Downscaled CMIP3 and CMIP5 Climate and Hydrology Projections: Release of Downscaled CMIP5 Climate Projections, Comparison with Preceding Information, and Summary of User Needs. U.S. Department of the Interior, Denver. Brownstein, J.S., Holford, T.R., Fish, D., 2003. A climate-based model predicts the spatial distribution of the Lyme disease vector Ixodes scapularis in the United States. Environ. Health Perspect. 111, 1152–1157. Brownstein, J.S., Holford, T.R., Fish, D., 2005. Effect of climate change on Lyme disease risk in North America. EcoHealth 2, 38–46. CDC (Centers for Disease Control and Prevention), 2008. Surveillance for Lyme disease – United States, 1992–2006. MMWR Surveill. Summ. 57, 1–9. CDC (Centers for Disease Control and Prevention), 2009. Summary of notifiable diseases: United States, 2007. MMWR Morb. Mortal. Wkly. Rep. 56, 1–94. CDC (Centers for Disease Control and Prevention), 2014. Reported Cases of Lyme Disease by Year, United States, 2003–2012, Available at: http://www.cdc.gov/ lyme/stats/chartstables/casesbyyear.html (accessed 05.07.14). Connally, N.P., Durante, A.J., Yousey-Hindes, K.M., Meek, J.I., Nelson, R.S., Heimer, R., 2009. Peridomestic Lyme disease prevention: results of a population-based case-control study. Am. J. Prev. Med. 37, 201–206. Dautel, H., Dippel, C., Kämmer, D., Werkhausen, A., Kahl, O., 2008. Winter activity of Ixodes ricinus in a Berlin forest. Int. J. Med. Microbiol. 298 (Supplement 1), 50–54. Diuk-Wasser, M.A., Gatewood, A.G., Cortinas, M.R., Yaremych-Hamer, S., Tsao, J., Kitron, U., et al., 2006. Spatiotemporal patterns of host-seeking Ixodes scapularis nymphs (Acari: Iodidae) in the United States. J. Med. Entomol. 43, 166–176. Diuk-Wasser, M.A., Hoen, A.G., Cislo, P., Brinkerhoff, R., Hamer, S.A., 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. Diuk-Wasser, M.A., Vourc’h, G., Cislo, P., Hoen, A.G., Melton, F., Hamer, S.A., et al., 2010. Field and climate-based model for predicting the density of host-seeking nymphal Ixodes scapularis, an important vector of tick-borne disease agents in the eastern United States. Glob. Ecol. Biogeogr. 19, 504–514. Eisen, L., Eisen, R.J., Lane, R.S., 2002. Seasonal activity patterns of Ixodes pacificus nymphs in relation to climatic conditions. Med. Vet. Entomol. 16, 235–244. Eisen, R.J., Eisen, L., Castro, M.B., Lane, R.S., 2003. Environmentally related variability in risk of exposure to Lyme disease spirochetes in northern California: effect of climatic conditions on habitat type. Environ. Entomol. 32, 1010–1018. Estrada-Pena, A., 2002. Increasing habitat suitability in the United States for the tick that transmits Lyme disease: a remote sensing approach. Environ. Health Persp. 110, 635–640. Fisman, D.N., 2007. Seasonality of infectious diseases. Ann. Rev. Public Health 28, 127–143. Gray, J.S., 1984. Studies on the dynamics of active populations of sheep tick, Ixodes ricinus L., in Co. Wicklow, Ireland. Acarologia 25, 167–178. Gray, J.S., 1985. Studies on larval activity of the tick Ixodes ricinus L. in Co. Wicklow, Ireland. Exp. Appl. Acarol. 1, 307–316. 622 A.J. Monaghan et al. / Ticks and Tick-borne Diseases 6 (2016) 615–622 Gray, J.S., 2008. Ixodes ricinus seasonal activity: implications of global warming indicated by revisiting tick and weather data. Int. J. Med. Microbiol. 298, 19–24. Gray, J.S., Dautel, H., Estrada-Pena, A., Kahl, O., Lindgren, E., 2009. Effects of climate change on ticks and tick-born dieses in Europe. Interdiscip. Perspect. Infect. Dis. 7, http://dx.doi.org/10.1155/2009/593232. Hay, L.E., Wilby, R.L., Leavesley, G.H., 2000. A comparison of delta change and downscaled GCM scenarios for three mountainous basins in the United States. J. Am. Water Resour. Assoc. 36, 387–397. Hayes, E.B., Piesman, J., 2003. Current concepts – how can we prevent Lyme disease? N. Engl. J. Med. 348, 2424–2430. IPCC (Intergovernmental Panel on Climate Change), 2013. Climate change 2013: the physical science basis. In: Stocker, T.F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S.K., Boschufng, J., Nauels, A., Xia, Y., Bex, V., Midgley, P.M. (Eds.), Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge/New York. Knutti, R., Masson, D., Gettelman, A., 2013. Climate model genealogy: generation CMIP5 and how we got there. Geophys. Res. Lett. 40, 1194–1199. Levi, T., Keesing, F., Oggenfuss, K., Ostfeld, R.S., 2015. Accelerated phenology of blacklegged ticks under climate warming. Philos. Trans. R. Soc. B 370, 20130556. Luber, G., 2014. Climate change and human health: NASA and Centers for Diesease Control and Prevention (CDC) collaboration. Geocarto Int. 29, 17–18. MacLeod, J., 1935. Ixodes ricinus in relation to its physical environment. II. The factors governing survival and activity. Parasitology 27, 123–149. MacLeod, J., 1936. Ixodes ricinus in relation to its physical environment. IV. An analysis of the ecological complexes controlling distribution and activities. Parasitology 28, 295–319. Masui, T., Matsumoto, K., Hijioka, Y., Kinoshita, T., Nozawa, T., Ishiwatari, S., et al., 2011. An emission pathway for stabilization at 6 W m−2 radiative forcing. Clim. Change 109, 59–76. Mather, T.N., Wilson, M.L., Moore, S.I., Ribeiro, J.M.C., Spielman, A., 1989. Comparing the relative potential of rodents as reservoirs of the Lyme disease spirochete (Borrelia burgdorferi). Am. J. Epidemiol. 130, 143–150. Maurer, E.P., Brekke, L., Pruitt, T., Duffy, P.B., 2007. Fine-resolution climate projections enhance regional climate change impact studies. Eos Trans. AGU 88, 504. McCabe, G.J., Bunnell, J.E., 2004. Precipitation and the occurrence of Lyme disease in the northeastern United States. Vector Borne Zoonotic Dis. 4, 143–148. Moore, S.M., Eisen, R.J., Monaghan, A.J., Mead, P., 2014. Meteorological influences on the seasonality of Lyme disease in the United States. Am. J. Trop. Med. Hyg. 90, 486–496. Moss, R.H., Edmonds, J.A., Hibbard, K.A., 2010. The next generation of scenarios for climate change research and assessment. Nature 463, 747–756. Mount, G.A., Haile, D.G., Daniels, E., 1997. Simulation of blacklegged tick (Acari: Ixodidae) population dynamics and transmission of Borrelia burgdorferi. J. Med. Entomol. 34, 461–484. Ogden, N.H., Lindsay, L.R., Beauchamp, G., Charron, D., Maarouf, A., O’Callaghan, C.J., Waltner-Toews, D., Barker, I.K., 2004. Investigation of relationships between temperature and developmental rates of tick Ixodes scapularis (Acari: Ixodidae) in the laboratory and field. J. Med. Entomol. 41, 622–633. Ogden, N.H., Maarouf, A., Barker, I.K., Bigras-Poulin, M., Lindsay, L.R., Morshed, M.G., O’callaghan, C.J., Ramay, F., Waltner-Toews, D., Charron, D.F., 2006. Climate change and the potential for range expansion of the Lyme disease vector Ixodes scapularis in Canada. Int. J. Parasitol. 36, 63–70. Ogden, N.H., Bigras-Poulin, M., Hanincova, K., Maarouf, A., O’Callaghan, C.J., Kurtenbach, K., 2008. Projected effects of climate change on tick phenology and fitness of pathogens transmitted by the North American tick Ixodes scapularis. J. Theor. Biol. 254, 621–632. Ogden, N.H., Radojevic, M., Wu, X., Duvvuri, V.R., Leighton, P.A., 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. Ostfeld, R.S., Canham, C.D., Oggenfuss, K., Winchcombe, R.J., Keesing, F., 2006. Climate, deer, rodents and acorns as determinants of variation in Lyme disease risk. PLoS Biol. 4, http://dx.doi.org/10.1371/journal.pbio.0040145. Perret, J.L., Guigoz, E., Rais, O., Gern, L., 2000. Influence of saturation deficit and temperature on Ixodes ricinus tick questing activity in a Lyme borreliosis-endemic area (Switzerland). Parasitol. Res. 86, 554–557. Piesman, J., 1989. Transmission of Lyme-disease spirochetes (Borrelia-Burgdorferi). Exp. Appl. Acarol. 7, 71–80. Rand, P.W., Lacombe, E.H., Elias, S.P., Lubelczyk, C.B., St Amand, T., Smith, R.P., 2010. Trial of a minimal-risk botanical compound to control the vector tick of Lyme disease. J. Med. Entomol. 47, 695–698. Riahi, K., Rao, S., Krey, V., Cho, C., Chirkov, V., Fischer, G., et al., 2011. RCP8.5 – a scenario of comparatively high greenhouse gas emissions. Clim. Change 109, 33–57. Randolph, S.E., 2004. Tick ecology: processes and patterns behind the epidemiological risk posed by ixodid ticks as vectors. Parasitology 129, S37–S65. Simon, J.A., Marrotte, R.R., Desrosiers, N., Fiset, J., Gaitan, J., Gonzalez, A., et al., 2014. Climate change and habitat fragmentation drive the occurrence of Borrelia burgdorferi, the agent of Lyme disease, at the northeastern limit of its distribution. Evol. Appl. 7, 750–764. Subak, S., 2003. Effects of climate on variability in Lyme disease incidence in the northeastern United States. Am. J. Epidemiol. 157, 531–538. Tagliapietra, V., Rosa, R., Arnoldi, D., Cagnacci, F., Capelli, G., Montarsi, F., et al., 2011. Saturation deficit and deer density affect questing activity and local abundance of Ixodes ricinus (Acari, Ixodidae) in Italy. Vet. Parasitol. 183, 114–124. Taylor, K.E., Stouffer, R.J., Meehl, G.A., 2012. An overview of CMIP5 and the experiment design. Bull. Am. Meteorol. Soc. 93, 485–498. Thomson, A.M., Calvin, K.V., Smith, S.J., Kyle, G.P., Volke, A., Patel, P., et al., 2011. RCP4.5: a pathway for stabilization of radiative forcing by 2100. Clim. Change 109, 77–94. Tran, P.M., Waller, L., 2013. Effects of landscape fragmentation and climate on Lyme disease incidence in the northeastern United States. Ecohealth 10, 394–404. USGCRP (United States Global Change Research Program), 2014. In: Melillo, J.M., Richmond, T.C., Yohe, G.W. (Eds.), Climate Change Impacts in the United States: The Third National Climate Assessment. U.S. Government Printing Office, Washington. van Vuuren, D.P., Stehfest, E., den Elzen, M., Kram, T., van Vliet, J., Deetman, S., et al., 2011. RCP2.6: exploring the possibility to keep global mean temperature increase below 2 ◦ C. Clim. Change 109, 95–116. Xia, Y.L., Mitchell, K., Ek, M., Sheffield, J., Cosgrove, B., Wood, E., et al., 2012. Continental-scale water and energy flux analysis and validation for the North American Land Data Assimilation System project phase 2 (NLDAS-2): 1. Intercomparison and application of model products. J. Geophys. Res. 117, D03109, http://dx.doi.org/10.1029/2011JD016048.