Geographical and environmental factors driving the increase in the Lyme disease vector Ixodes scapularis

. The population densities of many organisms have changed dramatically in recent history. Increases in the population density of medically relevant organisms are of particular importance to public health as they are often correlated with the emergence of infectious diseases in human populations. Our aim is to delineate increases in density of a common disease vector in North America, the blacklegged tick, and to identify the environmental factors correlated with these population dynamics. Empirical data that capture the growth of a population are often necessary to identify environmental factors associated with these dynamics. We analyzed temporally- and spatially-structured field collected data in a geographical information systems framework to describe the population growth of blacklegged ticks ( Ixodes scapularis ) and to identify environmental and climatic factors correlated with these dynamics. The density of the ticks increased throughout the study ’ s temporal and spatial ranges. Tick density increases were positively correlated with mild temperatures, low precipitation, low forest cover, and high urbanization. Importantly, models that accounted for these environmental factors accurately forecast future tick densities across the region. Tick density increased annually along the south-to-north gradient. These trends parallel the increases in human incidences of diseases commonly vectored by I. scapularis . For example, I. scapularis densities are correlated with human Lyme disease incidence, albeit in a non-linear manner that disappears at low tick densities, potentially indicating that a threshold tick density is needed to support epidemiologically-relevant levels of the Lyme disease bacterium. Our results demonstrate a connection between the biogeography of this species and public health.


INTRODUCTION
Vector-borne diseases are one of the most common classes of emerging and re-emerging infectious diseases and constitute a major public health threat (Taylor et al. 2001, Jones et al. 2008).
The frequent emergence of these diseases in human populations is often presumed to result v www.esajournals.org from anthropogenic changes to the environment. These environmental changes can promote increases in population sizes or geographic distributions of the pathogens or their vectors (e.g., Guerra et al. 2002, Lounibos 2002, Ogden et al. 2006, Ogden et al. 2009). Current and past population densities of these disease systems are often constrained by the environmental limits of the vectors due to their more direct interactions with abiotic factors (e.g., Brown et al. 1996, Krasnov et al. 2005, Härkö nen et al. 2010, Khatchikian et al. 2010a. Identifying the environmental and climatic features that allow or promote increases in density is essential to understand the basic biological and ecological processes of vectors. Further, the application of empirically validated biogeographic models allows accurate predictions of future population expansions and can constitute a critical component for the development of public health policies (Guerra et al. 2002, Gage et al. 2008, Diuk-Wasser et al. 2010, Kaplan et al. 2010, Khatchikian et al. 2010b.
Anecdotal and correlational evidence suggests that the population density of the blacklegged tick (Ixodes scapularis) has increased in the Northeastern and Midwestern United States over the previous few decades (e.g., Daniels et al. 1998, Estrada-Peña 2002, Ogden et al. 2008b). The presumed increase in I. scapularis density parallels the increase in disease incidence of several human diseases commonly vectored by this tick including anaplasmosis, babesiosis, and Lyme disease. Further, reports of increases in Ixodid tick densities are catholic across the northern hemisphere (Materna et al. 2005, Ogden et al. 2009, Parkinson and Berner 2009, Jaenson and Lindgren 2011, Jaenson et al. 2012 suggesting that common factors may be involved. Multiple hypotheses have been proposed to account for changes in the density of I. scapularis including changes in climate, land-use, and other factors that may have recently provided the ecological conditions necessary for colonization or population growth (e.g., Estrada-Peña 2002, Brownstein et al. 2005, Connally et al. 2006, Diuk-Wasser et al. 2010. However, assessing the relevance of fine-scale factors to regional-scale demographic patterns requires a wealth of fine temporal-and spatial-scale data. Further, the data must encompass the timeframe in which the demographic processes occurred to accurately assess the dynamics and to correlate environmental features with those dynamics. Temporally as well as spatially structured data are often necessary to identify the environmental and climatic factors that affect the abundance of populations undergoing dynamic shifts in density (Guisan andThuiller 2005, Sutherst andBourne 2009). Data collected during the dynamic phases of population growth allow for the development of models that explicitly incorporate time and space in addition to environmental and climatic variables. Such models cope with the inherent spatial and temporal heterogeneity of dynamic populations allowing for more accurate assessments of factors affecting the abundance of populations.
Here we present empirical data detailing the population expansion of the blacklegged tick over the previous decade. We use these finescaled, temporally and spatially structured data to build and empirically validate biogeographic models that identify the environmental factors correlated with the apparent population density increase. We analyze tick density data in a GIS framework using samples collected throughout the dynamic phase of population expansion. Tick density and environmental data were collected from 44 locations over a seven year period that corresponds to the timeframe in which ticks (and the most common disease that it vectors) noticeably emerged in the northern Hudson River Valley of the State of New York, USA ( Fig. 1, Appendix: Tables A1 and A2). We describe the dynamic pattern of tick population growth and identify environmental variables that account for these patterns. In this way, we provide a biogeographical framework that can serve as the foundation to further explore the basic ecology and demographic history of this medically important vector species.

Tick sampling
Nymphal and adult I. scapularis ticks were collected from a total of 44 sites across 22 counties in the Hudson River Valley, New York State (USA), from 2004 through 2010 (Fig. 1). Collection sites were chosen based on forested areas for which access was obtained. Not all collection sites were visited each year. Hostseeking ticks were collected by standardized dragging using a 1-m 2 piece of white flannel or canvas (Ginsberg and Ewing 1989). Tick density estimates were calculated as the number of ticks collected by drag sampling divided by the total area surveyed (m 2 ). Sampling locations were visited during June-August to estimate nymphal densities and during March-May and October-November to estimate adult tick densities. Density estimates of adults co-occurring during nymphal collection were included in the analyses. Adult ticks collected during the spring were considered the same cohort as those collected the previous fall in all analyses (see Tick phenology). Each sampling visit obtained a single density estimate for the appropriate tick stage. Tick collection efforts resulted in the collection of 7,140 nymphs (from 171 sampling visits) and 12,462 adults (from 238 sampling visits). Collection effort was evenly distributed across years except 2005 where only 20 density estimates were recorded.

Tick phenology
In northeastern North America, larval ticks hatch from egg masses and actively seek their first blood meal host, often a small mammal or bird, in late summer. After successfully feeding, larvae drop from hosts, molt to the nymphal stage, and overwinter in a quiescent state until they actively seek a second blood meal host, often a small mammal or bird, in early summer. Nymphal ticks that successfully feed drop from their host and molt to the adult stage which seeks a third host, often a white-tailed deer, in the fall. Adult females feed and mate on hosts and detach prior to producing eggs that hatch the following year. Adult females that fail to find a host in the fall overwinter and resume activities the following spring. The average life cycle of the blacklegged tick occurs over two years, with the overwhelming majority of time spent detached from hosts and directly subject to environmental elements and with little ability to move distances greater than 1 meter. The days spent on animal hosts-approximately three, five, and twelve for the larvae, nymphs, and adults, respectivelyprovide an opportunity for dispersal. More extensive descriptions of tick phenology are available in prior publications (e.g., Daniels et al. 1996, Ostfeld et al. 1996, Randolph 1998, Ogden et al. 2006, Ogden et al. 2008a).

Covariate selection
Tick density estimates may be affected by environmental conditions that modify questing behavior during sampling (e.g., Hubálek et al. 2003, Perret et al. 2003, Perret et al. 2004, Crooks and Randolph 2006, Devevey and Brisson 2012. To account for this source of error, environmental factors including temperature, humidity, wind conditions, and cloud cover were recorded at each sampling session. We included collection week to consider tick seasonal patterns (see Tick phenology). To assess factors that affect tick densities, we selected variables previously hypothesized to affect tick survival, rates of development, and/or habitat suitability. Mean annual temperature was utilized to consider tick development and survival rate (Lindsay et al. 1995, Brownstein et al. 2003, Ogden et al. 2004. Temperature and precipitations during the warmest quarter were included to consider desiccation risk during summer months (Ber- v www.esajournals.org trand and Wilson 1996, Brownstein et al. 2003). Precipitation during the coldest quarter was used to assess the potential effect of snow cover or excessive moist conditions at the soil level (Guerra et al. 2002). Temperature during the coldest quarter was incorporated to consider overwinter survival (Lindsay et al. 1995, Brownstein et al. 2003. Landcover indices were included to consider habitat quality (Lindsay et al. 1998, Guerra et al. 2002, Allan et al. 2003, Brownstein et al. 2005. Importantly, the geographic position (latitude and longitude) and collection year were included in all models to explain the variance over time and space that can be attributed to population growth (see also Spatial and temporal autocorrelation). The spatial and temporal variables are essential to capture the dynamic nature of the changing population.

Covariate data sources
Climatic data were extracted from Worldclim with 30-arcseconds resolution (Hijmans et al. 2005). Landscape data with 90 m resolution were obtained from the Land Cover Trends Project (Trends study year 2000; US Geological Survey, http://landcovertrends.usgs.gov/) to estimate the proportion of total landcover corresponding to forest, urban, and shrubs using a 7 3 7 pixels kernel. The total area of forest patches (minimum area 2 3 2 pixels) and ratio of forest perimeter to area were calculated using 50 m, 100 m, and 500 m search radii in Fragstats 3.3 (McGarigal et al. 2002). Landscape datasets were resampled to 30arcseconds resolution to match the climate data resolution for spatial predictions using cubic interpolation. All datasets were windowed to match our study area (75834 0 30 00 W, 738W; 418N, 458N). Point values for sample locations were obtained by point to raster interception. GIS processing was performed using IDRISI Taiga (Clark Labs, Worcester, Massachusetts, USA) and ArcGis 9.3 (ESRI, Redlands, California, USA) with the Geospatial Modeling Environment 0.4.0 plug-in (Spatial Ecology LLC, http://www. spatialecology.com/gme).

Model development and validation
Multiple regression models were developed independently for each tick developmental-stage (nymph and adult). All models were developed using a subset of the data that included all estimates from 2004 through 2008, while the remaining years (2009)(2010) were utilized only for model validation. Models were validated by linear regression of predicted to observed values (Rykiel 1996). We selected covariates, interactions between covariates (up to second order), and non-linear effects (up to second order) a priori to develop and evaluate competing models. Higher degree interactions and polynomial responses were not considered because of uncertainty in the biological interpretation. Akaike's Information Criteria (Akaike 1974) corrected for small sample size (Burnham and Anderson 2002) (AIC c ) was employed to select the best performing models. Significant effects of covariates were evaluated using Bonferonni-corrected multiple regression analyses (Bland andAltman 1995, Abdi 2007).
Nymphal and adult tick density estimates from collections (individuals per m 2 ) were extrapolated to individuals per 30-arcseconds square (;63 hectare) and ln transformed to achieve normally distributed data (for an extensive review of advantages and requirements of this approach see O'Hara and Kotze 2010). Four data points that had density estimates of zero could not be transformed and were removed from the analysis. Models that use count data (Poisson regressions and zero inflated Poisson regressions [Zeileis et al. 2007]) were explored to estimate biases introduced by the data transformation using the R package PSCL (R Development Core Team 2008). JMP 7.0 (SAS 2007) and SPSS 17 (SPSS 2008) were used for all remaining statistical analysis using Type IV sum of square.

Spatial and temporal autocorrelation
The geographic coordinates and year were included as covariates in all models to account for the spatial and temporal autocorrelation that is inherent in directional processes. Autocorrelation that is not accounted for by these covariates can reduce the effective degrees of freedom and can make the models prone to type I error. The residuals of the models were tested for autocorrelation using Moran's I in Passage 2.0 (Rosenberg and Anderson 2011) and checked for departures from normality using the D'Agostino Omnibus test Pearson 1973, D'Agostino et al. 1990).

Human disease incidence as a function of tick densities
We used the models with higher explanatory power to estimate predicted mean nymph and adult density estimates per county and year by averaging the predicted densities within pixels included in each county. The county-wide mean densities were correlated with the human Lyme disease incidence (per 100,000 habitant), the most commonly-reported disease vectored by I. scapularis (Appendix: Tables A1 and A2). Two additional diseases that I. scapularis vectorsanaplasmosis and babesiosis-have a very low incidence and thus could not be correlated with tick densities. After initial data exploration, we used piecewise correlations to identify hinge values (or break points) a posteriori in the correlations between the predicted nymph and adult density estimates (McGee andCarleton 1970, Neter et al. 1996). The method identifies the parameters of two correlations, one above and one below an iteratively determined hinge value, by finding the threshold value that minimizes the sum of squares of the errors around each correlation. We verified that the piecewise correlations minimize the total sum of squares of errors compared with a single first or second degree polynomial correlation.

RESULTS
Direct estimates of both nymphal and adult tick densities demonstrate that tick densities were greatest in lower latitudes of the study area and gradual decreased toward northern latitudes of the study area ( Fig. 2A). The observed patterns in tick density estimates across space are similar to the observed patterns in reported human Lyme disease cases (Fig. 2B). Additionally, tick density estimates increased annually across the region with the exception of 2010. Examination of the dataset revealed that observed densities in 2010 were unusually low. For example, the average density of nymphs was 0.022 per m 2 in 2010 in contrast with an average density of 0.036 nymphs per m 2 detected in all the other years combined. Despite the pattern irregularity in 2010, the trend of increasing densities was apparent in the study region.
The dynamic patterns in the tick density estimates observed in the raw data were cap-tured in our models (Tables 1 and 2). The statistical models combining geographical, temporal, seasonal, environmental, climatic, and landscape covariates accounted for the majority of the variance in nymphal and adult tick density data (R 2 ¼ 0.64 and 0.62). Models omitting covariates from any of these categories accounted for substantially less of the variance in the observed density data (Table 1; see Appendix:  Table A3). Variation in nymphal density estimates were best explained by a model that included geographic location, year, week of sampling (season), summer precipitation, minimum winter temperature, and proportion of forest cover (Table 1). Adult density estimates were best explained by a similar model that substituted winter precipitation for summer precipitation and included urbanization of the landscape (Table 1). Despite the simplicity of v www.esajournals.org these models, each explained more than 60% of the variation in the observed tick densities. Several competing models for both nymphal and adult density had a similar fit to the observed data (DAIC c , 2; see Appendix: Tables A3 and A4). These models differ from the bestfitting model only by alternative covariates that are highly-correlated with covariates in the bestfitting model. The coefficients estimated for each covariate, as well as the model predictions, are nearly identical among these best models indicating robustness to covariate selection. For example, replacing minimum winter temperature with annual mean temperature, covariates that are highly correlated (R 2 ¼ 0.71), in the bestfitting models for adult tick density estimates results in comparable coefficient estimates and similar fit to the data (R 2 ¼ 0.622 vs. R 2 ¼ 0.613, DAIC c ¼ 1.72). The models were insensitive to uneven stratified sampling effort by latitude (see Appendix: Table A5) and to effects of the log transformation of the response variable; models obtained with count data in Poisson regressions identified the same covariates as log-transformed models with equal probabilities, directions, and intensities (see Appendix: Table A6).  v www.esajournals.org All covariates included in the best-fitting models explaining nymphal density estimates are statistically significant in multiple regression analysis, and all but urbanization are significant in models explaining adult density estimates (Table 2). Interestingly, while both nymphal and adult tick densities increase annually, the year covariate is the predominant factor only for nymphal density estimates (explaining 30.79% of variance) and a significant although minor explanatory factor for adult density estimates (4.33%). Seasonal variation (week) accounts for the majority of variation in adult tick density estimates (59.17%) and a substantial fraction of the nymphal density estimates (19.76%). Geographic location explains 20% of observed variance in nymphal and 13% in adult density estimates. Precipitation and winter temperatures are correlated with tick density estimates (negatively and positively, respectively) and account for a small proportion of the variation in nymphal and adult tick density estimates (10% and 3%, respectively). Landscape covariates are also correlated with tick densities, explaining 10% and 15% of the variation in nymphal and adult density estimates, respectively.
The models that best explain nymphal density data from 2004 through 2008 successfully predicted the nymphal density estimates collected in 2009 and 2010 with 80% and 74% accuracy, respectively (Fig. 3A, B). The prediction for 2010 was consistently greater than the observed values although the general trend is maintained. The models developed with the 2004-2008 adult density estimate dataset had lower predictive power for 2009 adult density estimates (48%), primarily due to poor predictions for two locations (Fig. 3C). The model predictions for 2010 adult density estimates presented similar performance to those for nymphs (67% accuracy) and were similarly biased toward over-prediction (Fig. 3D).
The nymphal (Fig. 4A, B) and adult (Fig.  4C, D) density estimate models demonstrated good predictive power of trends over wide geographical areas with only the single bias in the 2010 dataset described above (Fig. 4D). The predicted surfaces of nymphal and adult tick density estimates differ in subtle details. For example, the nymphal density estimate model predicts a smooth surface with progressive transition between low and high predicted estimates (variance to mean ratio ¼ 0.61 for 2009 and 0.59 for 2010 [ Fig. 4B], suggesting a regular pattern). In contrast, the adult density estimate model predicts abrupt changes across the landscape (variance to mean ratio ¼ 1.32 for 2009 and 1.28 for 2010 [ Fig. 4D], suggesting a clustered pattern). The residuals from the linear models showed no departure from normality (D'Agostino omnibus test, a ¼ 0.05) and no evidence of autocorrelation, indicating that the spatial and temporal autocorrelation inherent in the dataset was accounted for by the covariates in the models. v www.esajournals.org Human Lyme disease incidence, the tick-borne disease with the greatest reported incidence in the region, was correlated with both nymphal and adult tick density estimates in a non-linear, threshold-type function (Fig. 5A, B). That is, human Lyme disease incidence remains at moderate levels at tick-collection densities below the optimal hinge values estimated using piece-wise correlation analyses (hinge values equal 1.13 and 0.71 individuals per hectare for nymphs and adults respectively). When densities are greater than the hinge values, human Lyme disease incidence is linearly correlated with local tick density estimates and explains 45% and 42% (for nymphs and adults respectively) of the variation in human Lyme disease rates.

DISCUSSION
In the northern hemisphere, the northern distribution and abundance of many organisms is limited by a complex combination of harsh environmental conditions and specific ecological characteristics. However, changes in one or several environmental factors can lead to dramatic increases in population densities. In this study, we used temporally-and spatially-structured data to characterize the population dynamics of the blacklegged tick, I. scapularis. Much of the heterogeneity in the density estimates of nymphal and adult ticks results from the annual population increases and to the marked latitudinal gradient (higher densities in southern areas and gradual decrease toward the north; Fig. 2A). The temporally-and spatially-structured data collected during dynamic changes in population size allowed confident identification of environmental and climatic factors that correlate with these changes. Including precise estimates of these factors into models accurately predicted tick density estimates in the near future (i.e., observed density estimates in the excluded data-sets of 2009 and 2010). These models can allow the assessment of future risks of human contact with diseases commonly vectored by this tick species.
The density estimates of both tick life-stages were greatest in the southern regions and decreased with increasing latitudes (Fig. 2A). Further, tick density estimates appeared to increase annually at most locations across the region. These trends, which are apparent in the sampling data, were formalized and quantified in statistical models that incorporated time, space, and environmental factors as covariates to account for the dynamics in tick density estimates. Approximately half of the variation in tick density estimates of both life stages was accounted for by the covariates year, geographic location, and seasonality indicating large population dynamics on these temporal and spatial scales. The remaining explainable variation in tick density estimates could be attributed to environmental factors that were not correlated with the temporal or spatial factors.
The climatic and landscape covariates in our model account for a substantial proportion of the variation in tick density estimates in our dataset Fig. 5. The relationship between model-predicted nymphal and adult density estimates (original units prior to log transformation are individuals/hectares) and the observed incidence of Lyme disease per county from 2004 to 2009. (A) Lyme disease incidence is significantly correlated with nymph density estimates greater than the hinge value (dotted line; hinge value ¼ 1.13 individuals hectare À1 ; R 2 ¼ 0.45, n ¼ 270, P , 0.0001) but not when estimates are lower than the hinge value (R 2 ¼ 0.02, n ¼ 40, P . 0.05). (B) Similarly, the correlation between Lyme disease incidence and local adult densities is only significant for values greater that the hinge value (hinge value ¼ 0.71 individuals hectare À1 ; R 2 ¼ 0.42, n ¼ 281, P , 0.0001) but not when the estimates are lower than the hinge value (R 2 ¼ 0.02, n ¼ 29, P . 0.05).
v www.esajournals.org (Table 2). Our results support the hypothesis that habitat suitability and climate are important factors determining tick population growth rates and densities, as previously suggested (e.g., Guerra et al. 2002, Brownstein et al. 2005). This hypothesis is further supported by a previous study that classified the northeastern portion of our study area as unsuitable for I. scapularis populations prior to 1990 (Estrada-Peña 2002) suggesting that climate change and anthropogenic changes to the landscape altered habitat suitability, allowing I. scapularis to establish and increase. However, the underlying causes of the recent discovery of I. scapularis populations in the northern Hudson River Valley remain unclear. Although improvements in habitat suitability may have allowed previously-established populations to grow from very low densities to detectable levels, it is also possible that I. scapularis populations have only recently colonized the region. Future analyses of the migratory patterns of I. scapularis in the Hudson Valley will help to discriminate between these hypotheses and further clarify the effects of climate change and landscape variation on tick demography.
Environmental factors such as extreme winter temperatures, summer or winter precipitation, and landscape variables were correlated with densities estimates suggesting that these factors directly regulate tick population dynamics. This hypothesis is supported by experimental studies demonstrating that mean and extreme temperatures, precipitation, and landscape features can directly affect survival, development rate, or reproduction rate in I. scapularis ticks (e.g., Lindsay et al. 1995, Lindsay et al. 1998). However, these environmental factors are also likely to affect the populations of vertebrate species that are the primary food sources for ticks (Ostfeld et al. 2006, Brunner et al. 2011. The composition and densities of vertebrate species have an important role in regulating tick densities, suggesting that environmental factors may also act indirectly. Decoupling direct effects of environmental factors on tick population dynamics from indirect effects requires a comprehensive catalog of the distribution and abundance of each of the tens to hundreds of potential vertebrate hosts species (Anderson 1988) across the study area and throughout the collection duration.
Unfortunately, the densities of few wildlife species have been estimated at local scales (Ostfeld et al. 2006) and none have been systematically estimated across the region of interest even at coarse-scales. The results of the current models can direct future work aimed at discriminating between the direct and indirect effects of the identified environmental factors.
The seasonal activity patterns of the tick lifestages dramatically affected the local density estimates for both the adult and nymphal ticks. The seasonality covariate accounted for this variation such that neither the estimates of annual increases in density across the region nor the estimates of the effects of the geographic location on population densities were affected by the sampling scheme. The effect of the seasonal activity pattern on immediate density estimates was expected as the seasonal activity patterns of the I. scapularis life-stages have been extensively documented (e.g., Daniels et al. 1996, Ostfeld et al. 1996, Randolph 1998, Ogden et al. 2006, Ogden et al. 2008a, Devevey and Brisson 2012. Variance in adult density estimates was strongly affected by the week of sample collection due to the two annual peaks of adult tick activity (Daniels et al. 1996). Surprisingly, excluding spring density estimates from the adult tick dataset results in a model that is nearly indistinguishable from that of the nymphs. In these models, less than 0.5% of the variation is accounted for by seasonality and other variables such as year (12.55%), geographical location (38.35%), and climatic covariates (32.1%) explain substantially more variance. Thus, our analyses suggest that the factors that affect the tick population densities are nearly identical regardless of the life-stage investigated, lending confidence to the robustness of these results.
The models that best described the data collected in 2004-2008 resulted in highly accurate predictions of tick densities estimates in 2009 and 2010 (Figs. 3 and 4). Thus, models that include broad-scale climatic and landscape variables are sufficient to predict the near-term densities of I. scapularis populations. Although all predictions were generally accurate, the model consistently predicted greater density estimates than were observed in the 2010 field-collection dataset suggesting that an important regional-scale factor in 2010 was not incorporated into the best v www.esajournals.org fitting models produced from the 2004-2008 dataset. We expect that a climatic factor that showed little variation between 2004-2008, and thus was not included in the best models, fluctuated dramatically in 2010 resulting in the low tick density estimates observed. Despite the over-prediction in 2010, the spatial trend in the model predictions was consistent with the observed data suggesting that these models can have a direct public health translation. However, extrapolation of the model predictions to areas that are beyond the geographic range used to develop and calibrate the models should be regarded with caution. For example, the northwestern region of our prediction map is projected by both statistical models to have the lowest tick density estimates. Although these expectations are compatible with the general model of tick density estimates and anecdotal evidence of current tick densities, the use of specific numerical values obtained in this region should be discouraged. Further data with longer time series and wider geographical extent are needed to determine the long-term and coarse scale accuracy of these models.
The density estimates of blacklegged ticks explain a substantial proportion of the variation in the human risk of contracting diseases commonly vectored by I. scapularis (Fig. 5). Interestingly, these results suggest that the correlation between human incidence and tick densities is non-linear and disappears at low tick density estimates. Considering that the majority of the human cases of Lyme disease are attributable to transmission from the nymphal stage, the detected correlation with adult densities estimates should be considered as a proxy of general tick densities in the area. These data could suggest that there is a threshold population density of I. scapularis needed to support epidemiologically-relevant levels of B. burgdorferi in vector tick populations as previously suggested ). However, B. burgdorferi can be detected at low prevalence in animal hosts suggesting that low-level transmission may occur even at low I. scapularis densities (Oliver et al. 2006), potentially explaining the sporadic Lyme disease cases in areas with few I. scapularis ticks. This result reveals the importance of a non-trivial parameter that needs to be incorporated in the development of epidemiological models that aim to predict human risk based on vector densities. Incorporating microbial pathogen population dynamics along with vector densities in a spatial and temporal framework will likely increase the accuracy and predictive power of models that assess human disease risk.
Studies to accurately determine areas of high human disease risks are critical to developing effective vector control and disease mitigation strategies. Such assessments must consider the complex combination of environmental conditions and ecological requirements that affect the distribution and abundance of medically relevant organisms. In this study, we provide three explicit contributions: First, we show that statistical models including temporal, geographic, climatic, and land-use variables accurately determine tick density estimates. In this way, we promote the inclusion of these factors into models built to predict future tick densities which may be generalizable to other habitats. Second, we correlate the increases in tick density with human disease risk over time and space, demonstrating the direct public health relevance of such models. Third, we provide explicit maps that highlight high density estimates tick areas. These maps allow the general population to exercise additional preventive measures in highrisk areas and allow vector control agencies to target these areas for intervention. Effective vector control not only reduces disease risk within the targeted area, but can also slow the rate of progression of the tick into new regions.

ACKNOWLEDGMENTS
C. E. Khatchikian and M. Prusinski contributed equally to this work. The authors would like to express their sincere gratitude to various NY State, county, local municipality, and private landowners for granting us use of their properties to conduct this research. We would also like to extend thanks to the following individuals for their assistance in collection and/or identification of ticks: J. Kokas, R. Falco, S. Kogut, J. H. Lee, M. VanDeusen, J. Hallisey, and a multitude of Entomological Assistants, student interns and county health department staff. Additional thanks to G. Lukacik for compiling human case numbers. We thank two reviewers for helpful comments on the manuscript. This work was supported by the Centers for Diseases Control and Prevention grant CK000170 and the National Institute of Health grants AI076342 and AI097137. Week 2 , air temp during collection int(1) , humidity during collection int(1) , winter prec, min winter temp, forest int(2) , urban int (2) 15 1.8734 0.659 Notes: Abbreviations are: prec, precipitation; min, minimum; temp, temperature. Superscript int marks interactions between variables. Superscript 2 indicates additional quadratic term.
v www.esajournals.org Table A5. In order to assess the robustness of our models to uneven sampling, we randomly selected subsets from each dataset to create datasets with locations evenly stratified by latitude. These subsets represent 60.9% and 41.42% of the dataset originally used for model development of nymphs and adults, respectively. The parameters estimates obtained with such unbiased subsets are remarkably similar to those obtained with the full dataset. The unbiased nymph model (Regression model for nymph density estimates using a randomly selected subset evenly stratified by latitude; R 2 ¼ 0.604, F ¼ 12.874, df ¼ 58, P ,0.0001), with roughly 60% of the data points, is practically identical to the original (see Table 2). Two of the covariables (longitude and summer precipitations) are not significant after Bonferroni corrections when compared with the original model. The unbiased adult model (Regression model for adult density estimates using a randomly selected subset evenly stratified by latitude; R 2 ¼ 0.76, F ¼ 21.151, df ¼ 60, P , 0.0001), with roughly 40% of the data points, is similar to the original (see Table 2). Although exact parameters estimates and significance differ slightly, the relative effects and directions of each one are similar. v www.esajournals.org Table A6. Equivalent Poisson multiple regression models to the linear multiple regression selected for nymph (df ¼ 101) and adult (df ¼ 159) densities. Statistics presented include the standard error term (SE), Z-value, and the probability value (P).