PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of tropmedLink to Publisher's site
 
Am J Trop Med Hyg. Jun 1, 2012; 86(6): 1062–1071.
PMCID: PMC3366524
Geographic Variation in the Relationship between Human Lyme Disease Incidence and Density of Infected Host-Seeking Ixodes scapularis Nymphs in the Eastern United States
Kim M. Pepin,* Rebecca J. Eisen, Paul S. Mead, Joseph Piesman, Durland Fish, Anne G. Hoen, Alan G. Barbour, Sarah Hamer, and Maria A. Diuk-Wasser
Fogarty International Center, National Institutes of Health, Bethesda, Maryland; Bacterial Diseases Branch, Division of Vector-Borne Infectious Diseases, Centers for Disease Control and Prevention, Fort Collins, Colorado; Yale School of Public Health, New Haven, Connecticut; Department of Community and Family Medicine, Dartmouth Medical School, Lebanon, New Hampshire; Departments of Microbiology and Molecular Genetics and Medicine, University of California, Irvine, California; Department of Fisheries and Wildlife, Michigan State University, East Lansing, Michigan; Department of Veterinary Integrative Biosciences, Texas A&M University, College Station, Texas
*Address correspondence to Kim M. Pepin, Fogarty International Center, National Institutes of Health, 31 Center Drive, MSC 2220, Bethesda, MD 20892-2220. E-mail: uzu2/at/cdc.gov
Received October 10, 2011; Accepted February 26, 2012.
Prevention and control of Lyme disease is difficult because of the complex biology of the pathogen's (Borrelia burgdorferi) vector (Ixodes scapularis) and multiple reservoir hosts with varying degrees of competence. Cost-effective implementation of tick- and host-targeted control methods requires an understanding of the relationship between pathogen prevalence in nymphs, nymph abundance, and incidence of human cases of Lyme disease. We quantified the relationship between estimated acarological risk and human incidence using county-level human case data and nymphal prevalence data from field-derived estimates in 36 eastern states. The estimated density of infected nymphs (mDIN) was significantly correlated with human incidence (r = 0.69). The relationship was strongest in high-prevalence areas, but it varied by region and state, partly because of the distribution of B. burgdorferi genotypes. More information is needed in several high-prevalence states before DIN can be used for cost-effectiveness analyses.
Lyme disease is the most common vector-borne disease in the United States, with more than 25,000 confirmed cases in 2009.1 Incidence of reported cases continues to rise. Between 2005 and 2009, mean incidence in the 13 states with highest incidence had increased from 29.6 ± 10.6 per 100,000 in 2005 to 49.6 ± 15.5 per 100,000 in 2009, whereas in 11 states with lower incidence, mean incidence has increased from 1.3 ± 0.7 to 2.3 ± 1.7 per 100,000.2 Efforts to stop this emergence have been hampered by the complex ecology of the disease and the lack of effective and affordable interventions.37
In the Northeastern and Midwestern United States, Lyme disease is caused by the bacterium Borrelia burgdorferi sensu stricto and transmitted to humans by Ixodes scapularis, primarily nymphs.8 Humans acquire the disease mainly during the months of May and July after contact with the vector during outdoor activities in tick habitat, which is characterized by closed-canopy deciduous and mixed forest. Three previous studies have shown a positive relationship between human cases and the density of B. burgdorferi-infected ticks, but these studies were limited in geographic scope to Connecticut and Rhode Island, which are both states with very high prevalence.911 Another two studies showed a positive correlation between tick abundance and human cases, again focusing on specific regions (Wisconsin12 and Westchester, New York13). Two more location-specific studies have examined the link between numbers of ticks sampled by passive surveillance and cases in Rhode Island and Maine.14,15 To our knowledge, these studies are the only studies that have attempted to quantify and/or characterize the relationship between acarological risk and human cases in the eastern United States. However, an understanding of this relationship and how it varies geographically is important for models used in identifying optimal strategies for implementing vector- and reservoir-targeted interventions.
In 2009, 95% of Lyme disease cases in the United States were reported from the Northeast and Midwest.2 However, there is large variation in incidence among states where transmission of B. burgdorferi occurs, possibly because of differences in reservoir host ecology,1618 climate,1923 rates of human contact,13 genetic variation among strains,2326 or differences between states in reporting practices. Thus, relationships between acarological risk and incidence may differ among geographic regions. Investigating the differences is an important stepping stone to a mechanistic understanding of the underlying drivers linking density of infected nymphs (DIN) and human incidence. Improved knowledge of the relative importance of individual drivers in different locations might aid in tailoring prevention and control recommendations to different geographical regions.
We have previously reported on a large-scale study to measure and develop models of the density of I. scapularis nymphs and B. burgdorferi prevalence in 36 states (2,411 counties) at an 8 × 8-km spatial scale.2729 Here, we used the raw data (rDIN and raw density of nymphs [rDON]) to examine the strength of the relationship between DON, DIN, and human incidence at the county level to evaluate whether the model-based predictions of acarological indices (mDON and mDIN) significantly correlated with reported human Lyme disease incidence and identify factors that modify the relationship between DIN and human cases. We then used the estimated indices to investigate geographical differences in the relationship between acarological risk and county-level incidence. To our knowledge, this study is the first study to examine the relationship over such a broad spatial extent and evaluate geographic differences in how DIN translates to human incidence.
Human case data.
Human case data were obtained from case reports to the Centers for Disease Control and Prevention (CDC) by state and territorial health departments as part of the National Notifiable Diseases Surveillance System (NNDSS); detailed information on methods and case definitions has been published previously.30 For each county in the United States, the number of cases between 2004 and 2006 was averaged to get the mean number of cases per county during the time period that tick data were collected (Figure 1C shows these data as incidence per 100,000). Note that cases are reported based on the county of residence, and thus, it is possible that transmission occurred outside the county of residence if a person had traveled recently.
Figure 1.
Figure 1.
Map of study sites, estimated acarological risk, and human incidence. (A) Counties are color-coded by mDON/km2; unique study sites are shown as dark red dots. (B) Counties are color-coded by mDIN/km2. (C) Counties are color-coded by incidence/100,000 (more ...)
Density of nymphs.
Host-seeking I. scapularis nymphs were collected by drag sampling at 301 sites in the United States east of the 100th meridian between 2004 and 200627,28 (Figure 1A). All sites were closed-canopy deciduous forests in state parks, state forests, or other natural areas with public access. The spatial distribution of sites included 36 states (286 counties), which we classified based on the incidence of human cases as follows: > 10/100,000, DE, CT, ME, MD, MA, NH, NJ, NY, PA, RI, VT, WI, and MN (high); 0.5–5/100,000, VA, WV, OH, SC, NC, IL, IN, IA, MI, and ND (low); < 0.5/100,000, AL, AR, FL, GA, KS, KY, MS, MO, NE, OK, SD, TN, and TX (rare). Each site was dragged three to six times in summer between late May and September, each time covering 1 km2 per visit. Visits were spaced as uniformly as possible throughout this time frame in each location to prevent detection bias because of sampling time (Supplemental Figure 1) (note that there were several sites where this spacing was not possible, but they were not localized in a single area and thus, should not consistently bias the results). Daily nymphal density per 1,000 m2 (i.e., DON) was calculated as the area under the nymphal frequency curve between the first and last sampling dates divided by the total number of days during this time period. Although nymphs tend to be most active in May and June in the northeast, we continued to sample through late summer (August and September), because (1) peak nymphal activity dates could differ across a large geographic area sampled, (2) human cases are most prevalent in late summer, and (3) humans tend to be most active in forested areas during summer holidays (May and September). In most counties, the same site was dragged each year except in three locations (Cumberland, ME; Suffolk, NY; and Chippewa, WI), where different sites were visited in different years. Because each site was considered to be representative of the county in which it was situated, daily nymphal densities from multiple years and/or all sites within the same county were averaged to get the nymphal density per 1,000 m2 during the entire study period. Density estimates were scaled to the county level by multiplying by the proportion of deciduous and mixed forest (i.e., I. scapularis habitat31) for the county where each site was located. This area was estimated using the ArcMap Zonal Statistics tool (ESRI, Redlands, CA) based on the US Geological Survey National Land Cover Database (http://www.mrlc.gov/nlcd.php). Hereafter, we refer to this estimate of DON as rDON. Note that accounting for land cover in the calculation of rDON was necessary to compare the tick data with human case data (which were available at the county level), because there is high variability between counties in forest coverage (amount of I. scapularis habitat).
A previous study28 developed a regression model for predicting DON at an 8 × 8-km resolution throughout the eastern United States. This model assumed a zero-inflated negative binomial error structure and included climate, landscape, and spatial structure as covariates. Specifically, elevation, monthly vapor pressure deficit, annual amplitude of the maximum daily temperature, and annual phase of the minimum daily temperature were covariates in the zero-inflated component of the model, whereas the annual amplitude of the normalized vegetation index and spatial autocorrelation were covariates in the negative binomial component. For the current study, we calculated the average DON by county using the previously developed model. We refer to model-derived estimates of DON as mDON. Nine counties were excluded because of missing climate or landscape data that were used as covariates. This lack of data was because the climate/landscape data, obtained from National Aeronautics and Space Administration (NASA), consisted of 8-km grids that were clipped to minimize the number of pixels that contained both land and water. Some counties that were coastal and small did not either contain data or have a sufficient number of pixels to calculate an average value. These 10 counties included Dukes, Nantucket, Suffolk, and Barnstable in Massachusetts; Bronx, Kings, New York, and Richmond in New York; and Bristol and Newport in Rhode Island.
Density of infected nymphs.
Detection and quantification of B. burgdorferi in nymphal tick specimens was performed by species-specific quantitative polymerase chain reaction (PCR) as previously described.23 Individual ticks were screened for the presence of B. burgdorferi. Nymphal infection prevalence (NIP) was calculated as the proportion of nymphs that tested positive for B. burgdorferi. Using a binomial probability distribution, it was expected with 95% confidence that at least 1 of 14 ticks would be infected if the site harbored B. burgdorferi-infected nymphs with an infection prevalence of 0.20. Thus, the lower limit for pathogen detection was 14 ticks, and sites with fewer ticks were excluded from analyses with rDIN. In the high and low prevalence areas (for which we conducted the analyses using rDIN and mDIN), the average number of ticks collected per site was 77.6 ± 21.6 (2 standard errors of the mean [SEMs]), meaning that most sites had at least 55 ticks. Yearly prevalence per county was calculated as the mean prevalence for all sites within the same county during the study period. The density of infected nymphs (rDIN) was calculated as the product of rDON and prevalence.
Similar to mDON, model DIN (mDIN) was estimated at the county level using a model that was developed previously for predicting DIN at a finer spatial scale29 (Figure 1B). The model used for predicting DIN was similar to the DON model, except that maximum daily temperature was excluded from the zero-inflation component and the largest forest patch index (which is a surrogate for fragmentation describing the size of the largest forested patch in a county) was substituted for the normalized difference vegetation index in the negative binomial component. Predictions were made for all counties within the study region, except for the nine mentioned above for which there was no climate/landscape data.
Additional covariates in the model of rDIN and human incidence.
We included B. burgdorferi genotype frequency, climate, and forest fragmentation variables in a model selection routine that included rDIN as a fixed component of the model (described in Statistical analyses) to examine whether these variables accounted for some of the variation in incidence left unexplained by rDIN. Our logic was that these variables could impact human behavior or B. burgdorferi transmission from vectors to humans. For example, the 16S-23S rRNA intergenic spacer restriction fragment length polymorphism sequence types (RSTs) have been previously linked to differential capacity for hematogenous dissemination in humans and severity of Lyme disease symptoms26,32 as well as disease severity in animal models and transmission rates in reservoir mice.25,33 Because the frequency of these types varies geographically,23 we hypothesized that differences in pathogenicity and transmissibility between RST types may explain part of the geographic variation in Lyme disease incidence. As for the fragmentation and weather variables, we are unaware of work that has directly tested whether they could impact human contact with vectors or transmission rates, but weather has been hypothesized to impact transmission.34 Furthermore, the direction of the relationship between fragmentation and DIN compared with fragmentation and human incidence is conflicting,17 suggesting that fragmentation could act directly on incidence through some unexplained mechanism. Thus, we incorporated weather and fragmentation variables from a hypothesis-generating standpoint. These variables included the amplitude of county-level maximum temperature (TMAXmag) derived through temporal Fourier transformation, the mean area of forest patches (mean patch area [MPA]), and the mean area occupied by the largest forest patch (largest patch index [LPI]).
The proportion of each of the three RST types was estimated for 29 sites (thus, there was some missing data in this covariate, because the others were N = 40) as described previously.35 In brief, real-time PCR with probes specific for RST 1 or RST 2 strains was carried out for all individual B. burgdorferi-positive tick extracts. A nymph was characterized as containing an RST 3 strain if it was positive by species-specific PCR for B. burgdorferi but negative by PCR specific for RST 1 or RST 2. These data were previously reported in the work by Gatewood and others,23 which found that the frequency of RST types was associated with variability in the synchronization of I. scapularis life stages.
Statistical analyses.
All statistical analyses were conducted in R and Matlab. We evaluated the relationship between rDON and cases (N = 286 counties in the full raw data) by negative binomial regression on different subsets of the data (based on human incidence). County population size was included as an offset, and thus, the model essentially predicts incidence (frequency of infection per county with weighted effects from population size). Inspection of residuals showed that a negative binomial model was better for rDON, but a Poisson model with square-root transformed case data was equally good for rDIN (N = 40; only sites with > 14 ticks were included) analyses. Thus, we used a Poisson model for all DIN analyses; the negative binomial model tended to largely overestimate case numbers in high prevalence areas, and the Poisson model is simpler with one less parameter.
We evaluated performance of mDON and mDIN by fitting them to incidence in counties for which there were observed data (N = 40) using a Poisson model. We compared these fits to those fits with rDON and rDIN by Akaike Information Criterion (AIC). Lower scores indicate better model fits (a two-point difference is significant). Because both mDON and mDIN fitted the incidence data better than rDON and rDIN, we used the mDIN data to investigate spatial differences in the relationship between DIN and human cases. Thus, we assumed that mDIN in the 40 counties where rDIN was measured was similarly representative of mDIN in the other 1,141 counties in high- and low-prevalence regions where rDIN was not measured. We fitted data from all counties in high- and low-incidence states to cases (N = 1,181) using mDIN as a continuous covariate and state as a factor. We also fitted mDIN to cases for each state individually and compared predictions from the full model to those predictions from the state-specific models by using the likelihood function to calculate AIC for each set of model predictions in each state. Performance of the model predictions against observed cases was evaluated by Spearman's rank correlation r.
We investigated geographical effects on the relationship of mDIN and incidence by including region (east versus west) as a factor. East included states east of and including Ohio (classification of the 23 states included is indicated in Results). To identify factors that were responsible for modifying the relationship between DIN and incidence regionally, we conducted model selection with rDIN in the model to start and six covariates: RST types 1–3, LPI, MPA, and TMAXmag (described above). Only main effects were examined. Variables were selected by forward selection and 10-fold cross-validation (where coefficients are estimated from 90% of the data and used to predict the other 10%) using the χ2-distributed difference in deviance for determining whether a variable is kept. This method of model selection decreases overfitting by requiring that candidate variables contribute to significantly accurate out of sample predictions for their incorporation into a model. Competing models generated by the selection procedure were compared by AIC. Last, to understand the hierarchical structure of the effects of region and other factors on the relationship between rDIN and incidence, we conducted a recursive partitioning analysis by conditional inference36 implemented using the party package in R. This method is non-parametric, and it identifies covariates that explain the most significant amount of response variable (i.e., incidence) variation by repeatedly splitting the response variable into two groups. The resulting tree indicates the hierarchical nature of the covariate effects on rDIN.
Relationship of rDON and cases in different regions.
Overall, rDON was significantly correlated to reported incidence (Spearman's r = 0.66, N = 286). Next, we compared the relationship of rDON with cases in regions where the incidence is rare, low, or high to investigate whether vector density could explain case variability in any or all of these areas (Figure 1A versus C). We found that, within rare areas, rDON did not explain variation in human Lyme disease incidence across counties (P = 0.99; AIC change over the null model with only an intercept was -0.7, which is less than 2) (Table 1). However, rDON explained a significant amount of variation in human incidence in each of the low- and high-incidence areas (P < 0.03 and P < 0.0001, respectively; ΔAIC = 80 and 13, respectively) (Table 1).
Table 1
Table 1
Fit of DON in different regions
Addition of B. burgdorferi prevalence data.
Next, we examined whether the incorporation of B. burgdorferi prevalence data explained more of the variation in incidence relative to vector abundance alone (Figure 1B versus C and Table 2). For the raw data, all covariates explained a significant amount of the variation in incidence (ΔAIC was greater than two above the null model) (Table 2), but much was left unexplained (Spearman's r = 0.57 for rDON and 0.60 for rDIN). rDIN performed only slightly better than NIP, despite the apparently larger variation in rDIN by location (Supplemental Figure 2).
Table 2
Table 2
Comparison of raw data with model data covariates
Dissecting the roles of additional factors.
Because much of the variation in human incidence remained unexplained by the model with rDIN (r = 0.6) (Table 2), we identified other factors that may be needed for interpreting human incidence in terms of rDIN. First, we examined whether accounting for geographical region (Northeast versus Midwest) improved the fit of rDIN and incidence. For this examination, we included region as a two-level factor and examined if either the intercept or slope of the relationship of rDIN and incidence differed by region. We found that there was no difference in the slope of the relationship between regions but that the relationship between rDIN and incidence had a significantly higher intercept in the Northeast relative to the Midwest (P < 0.0001). Figure 2A and Table 3 show statistics for the reduced model that excludes the non-significant interaction term between region and rDIN, and more of the variation was explained (r = 0.77). Figure 2A and B highlights that similar levels of rDIN occurred in both regions, but mean levels of incidence are lower in the Midwest relative to the Northeast region (Figure 2A).
Figure 2.
Figure 2.
Factors influencing the relationship of DIN and cases. (A) Relationship of rDIN and incidence with geographic region included in the model (model details in the text). (A–C) Eastern counties are in black, and western counties are in grey. (B) (more ...)
Table 3
Table 3
Comparison of parameter estimates and model fits between a model that takes regional differences into account and the best model without region effects
When region was excluded from the model, the best model of rDIN and incidence included the frequency of RST 3 and TMAXmag (r = 0.74) (Figure 1C and Table 3). Both variables were negatively related to incidence (Table 3), indicating that a higher frequency of RST3 strains and a larger difference between summer and winter temperatures are associated with lower incidence in humans. Although RST 1, RST 2, and the fragmentations indices (MPA and LPI) were not selected as the most significant modifiers of the rDIN–human incidence relationship, LPI and RST1 were significant in models that contained only rDIN and one other variable (LPI, P < 0.031; RST 1, P < 0.0002). The regression trees in Figure 2D and E illustrate the hierarchical structure of the most significant combination of factors explaining the variation in incidence. The major factors that separate high from low incidence irrespective of rDIN are region (Figure 2D, top node) and RST 3 frequency (Figure 2E, top node). In the Midwest, where RST 3 frequencies are higher, incidence is lower, although rDIN occurs at similar levels as in the Northeast (Figure 2B). After variation because of region or RST 3 is segregated, in the Midwest, there is more variation in incidence that can be significantly classified into two groups (high versus low) according to levels of rDIN (second nodes). The second node does not occur for the Northeast, because the number of points with low rDIN was not high enough (minimum number for division was 10). Although TMAXmag was also significant, with a lower value in the Northeast relative to the Midwest (Figure 2B), it did not show significant partitioning in the regression tree analysis.
Evaluation of model-estimated acarological indices: mDON and mDIN.
We compared the performance of previously estimated acarological indices, mDON and mDIN, with rDON and rDIN by fitting all four covariates individually to the counties for which there were empirical observations. Both estimated indices performed significantly better than the raw data (Table 2) (pseudo-R2 showed 11% and 26% improvements in fit over rDON and rDIN, respectively), and mDIN (r = 0.69) fit the incidence data better than mDON (r = 0.58). This finding shows that the estimated indices performed better at explaining variation in incidence, likely because they account for effects of other factors (i.e., landscape fragmentation, spatial autocorrelation, and weather) that are known to affect I. scapularis and B. burgdorferi distributions1618,28 and their interactions.23
Spatial variation in the relationship between mDIN and human incidence at the state level.
Because mDIN fit the incidence data best and provided acarological risk data for all counties within the 23 states in high- and low-prevalence areas (listed in Tables 4 and and5),5), we used this index to examine whether differences in the relationship of mDIN and human incidence occurred at a finer spatial scale: the state level (rather than the east versus west). We investigated state-specific differences in the relationship between mDIN and incidence through two model structures: a national-level model that included state as a factor versus separate models of mDIN for individual states. We found that there was a significant effect of state in the full model (P < 0.0001) and that the state-level models performed better than the full model in all states (Tables 4 and and5,5, compare AIC for full with state). Lack of fit of the full model was not specific to any geographic location (Figure 3B –I) or levels of mDIN or incidence (Table 5, columns 6–9).
Table 4
Table 4
Model including all 23 states
Table 5
Table 5
Models for individual states
Figure 3.
Figure 3.
Relationship of DIN and cases by state. The full model (A) included DIN as a continuous predictor and state as a 23-level factor. Because state was highly significant, we also modeled the relationship of DIN and human cases individually for each state (more ...)
Of 23 states, 12 states showed a significant relationship between mDIN and human incidence, with an average correlation between model predictions and incidence of 0.55 (range = 0.2–0.86) (Tables 4 and and5).5). The other 11 states showed no significant relationship. Three of these states, Delaware, Rhode Island, and Connecticut (Table 4), had a very low number of counties, and thus, their lack of fit is likely because of small sample size. For states where the relationship was significant, the fit was good (r > 0.6) in 4 of 12 states (New York, Pennsylvania, Massachusetts, and Minnesota), weak (r = 0.4–0.6) in 5 more states (Wisconsin, Virginia, Iowa, Illinois, and Michigan), and poor (r < 0.40) in the other 3 states (Maine, West Virginia, and Indiana), despite the high significance of mDIN in the models (Tables 5, bold P values and Supplemental Figure 3). Furthermore, the slopes of the relationship were significantly different among the states where the relationship was significant, including an unexpected negative relationship for Massachusetts (Figure 3B–I and Table 5, bold b-values and SE), highlighting that there are quantitative differences between states, even where the relationship is highly significant.
We quantified the relationship of acarological risk and human Lyme disease incidence for all of the known US range of I. scapularis-borne B. burgdorferi. Our study is the first to examine this association over such a large geographic range. We found that both rDON and rDIN explained a statistically significant amount of variation in human incidence in low- and high-incidence areas but that neither explained variation in incidence in rare areas. Restricting the analysis to low- and high-incidence areas, we found that the relationship of mDIN and incidence varied quantitatively by geographic location, that mDIN alone explained a statistically significant amount of variation in human incidence at the county level in only 12 of 20 states (including only states with more than eight counties), that the fit of the relationship was good in only four states (New York, Pennsylvania, Massachusetts, and Minnesota), and that fit was reversed in one of these states (Massachusetts). A subanalysis of 40 high incidence counties with NIP data identified some of the factors that contribute to geographic differences in the relationship of rDIN and incidence, most notably, genetic composition of B. burgdorferi.
One factor that has been proposed to modify the relationship between DIN and reported cases is variation in B. burgdorferi genotype frequency.26,32 Differences in virulence37 and transmission efficiency25,38 may lead to differences in infection risk and/or reporting rates by humans infected with different strains. Thus, areas with low incidence but high rDIN could be explained by an overrepresentation of strains with a low virulence in or transmission to humans. Although more studies are needed to confirm this hypothesis, our regression analyses support it for RST 3. We found that (1) incidence data could be partitioned similarly by region or RST 3, (2) RST 3 was more frequent in the region with lower incidence but similar rDIN, and (3) there was a negative relationship between RST 3 and incidence estimated by the Poisson regression model. Although the other two RST types each showed different regional distributions, neither of these two types alone caused deviations from the mean relationship of rDIN and reported incidence. Nevertheless, because higher levels of RST 3 (a genetically diverse group defined as non-RST 1 or 2) led to lower than expected incidence based on rDIN, by extrapolation, the sum of RST 1 and RST 2 frequencies was correlated with higher than expected incidence based on rDIN. Thus, our results suggest that the monophyletic RST 1 and 2 groups39 cause higher than expected reported cases based on rDIN, whereas a high prevalence of other strains results in fewer reported cases. More generally, our results highlight that genotypic differences are important for predicting human incidence spatially through acarological risk indices. However, more analyses of the genetic determinants of virulence (and transmission) in humans and the spatial dynamics of B. burgdorferi genetic diversity are needed to understand the effects of B. burgdorferi genotype on reported incidence.
In addition to B. burgdorferi genotype, there are other factors not measured in our study that could contribute to geographical differences and weak associations in the relationship between mDIN and human incidence. One possibility is human movement. Although peridomestic exposure accounts for most cases of Lyme disease,12 some cases are a consequence of travel to other counties for outdoor work or recreation. This movement could explain the negative relationship between mDIN and cases that was found for Massachusetts. For example, if Massachusetts residents tend to get infected in counties that are densely forested but live in counties (i.e., where cases would be reported) that are mainly urban (i.e., low estimates of mDIN), then such a negative relationship could be observed. More accurate spatial risk predictions might be obtained if case reports included probable place of exposure, although accurate estimates are very challenging to obtain.40 Also, we found that spatial scale was crucial for interpreting the relationship between mDIN and incidence, because we showed that there are quantitative differences regionally (Midwest versus Northeast) as well as between states within these regions. It is possible that the county scale is appropriate in states with larger forested patches, such as New York and Pennsylvania (which showed relatively strong positive relationships between mDIN and incidence), whereas a finer scale is necessary in highly residential states with very dispersed forest cover, such as New Jersey and Maryland (which showed no significant relationship). For example, a study that measured DIN in forested areas in six Rhode Island towns found that DIN explained 97% of the variation in incidence among those towns.9 However, most of our study sites were not located in residential areas, and the incidence data corresponded to a much broader range relative to the acarological data, which increases the variance in the relationship between mDIN and incidence. Last, geographical differences in Lyme disease control efforts at the residential or personal level could also contribute to variation in the relationship between mDIN and human incidence.
There are also some experimental design caveats that should be noted when interpreting our results. First, our estimates of mDIN were based on data collected exclusively in public forested areas. Thus, our estimates of mDIN in counties that are mainly residential could be inaccurate. Second, sampling at different sites was done at different times between May and early September, and although most sites were sampled at least five times, some sites were only sampled three times. Although these sampling design caveats likely account for some of the unexplained variation in the relationship between DIN and human incidence, they were not spatially systematic (i.e., they occurred randomly with respect to state and region) and thus, should not introduce systematic bias in our results. Third, this analysis is based on original field efforts that were optimized for collections of I. scapularis nymphs based on phenology in the northeastern United States (emergence in late May); therefore, areas in which nymphs become active earlier would have inaccurate estimates of nymphal densities. For example, a previous study found that nymphs were observed as early as April in the south.41 However, even with intensive sampling, only a few specimens were recovered at this time; thus, this difference should not introduce much bias in our results. There are also two other factors unaccounted for in our sampling design that could cause geographic variation in human incidence between northern and southern areas. First, because we only sampled I. scapularis, the existence of an alternative vector might uncouple the relationship between mDIN and incidence. However, this possibility would only occur if the alternative vector caused a significant proportion of human cases, and there is no evidence to support this case. Second, it has been reported that nymphs in southern areas very rarely feed on humans,41 although it is still unclear whether this difference from northern areas is because of the low densities of I. scapularis nymphs or an actual feeding preference. Future environmental surveillance should aim to quantify DIN in residential areas as well as other habitats and potential vectors to better understand how DIN translates to incident human cases across varying ecological conditions.
Reservoir- or vector-targeted control of human Lyme disease relies on understanding how DIN corresponds to human incidence geographically. Additionally, a mechanistic understanding of the link between DIN and incidence should be sought, because it is important for model frameworks that anticipate changes in the geographic range of human incidence and evaluate the cost efficacy of all types of prevention methods over time. Our multistate analysis showed large variation in the relationship of DIN and incidence between states, with only three states showing a reliable quantification of the relationship. To more accurately quantify how acarological risk translates to human Lyme disease risk spatially, longitudinal environmental surveillance of several factors at spatial scales finer than the county level would be needed. Unfortunately, such studies are not likely to be cost-effective if conducted using a random sampling design. Thus, future research to determine how acarological risk translates to human incidence spatially should first aim to identify probable sites of transmission to humans over a broad spatial scale through epidemiological studies that monitor human movement (preferably to and from areas where acarological risk is quantified). Intensive environmental sampling in some of these sites would help in quantifying the relationship between DIN and human incidence more cheaply.
Supplementary Material
Supplemental Figures.
ACKNOWLEDGMENTS
The authors thank Graham Hickling, Jean Tsao, Uriel Kitron, Edward Walker, Roberto Cortinas, Jonas Bunikis, and Michelle Rowland, who were key collaborators in the data collection and development of the acarological risk map. The authors thank the 80 field assistants who participated in the data collection; Tim Andreadis, Katherine Watkins, Laura Krueger, Jessica Payne, Elizabeth Racz, Kelly Liebman, Liza Lutzker, and David Boozer for tick identification and logistic support; Paul Cislo, Gwenaël Vourc'h, and Robert Brinkerhoff for statistical and data analyses support; Carlos Diuk for database support; Forrest Melton, Brad Lobitz, and Andrew Michaelis for technical assistance; and Dennis Grove, Lindsay Rollend, and Russell Barbour for arranging collection permits. This project was funded by Centers for Disease Control and Prevention–Division of Vector-Borne Infectious Diseases Cooperative Agreement no. CI00171-01. K.M.P. was supported by the Research and Policy in Disease Dynamics (RAPIDD) program of the Science and Technology Directorate, the US Department of Homeland Security, and the Fogarty International Center, National Institutes of Health. A.G.B. was supported by National Institute of Allergy and Infectious Disease Grant AI-065359.
Footnotes
Authors' addresses: Kim M. Pepin, Fogarty International Center, National Institutes of Health, Bethesda, MD; and Bacterial Diseases Branch, Division of Vector-Borne Infectious Diseases, Centers for Disease Control and Prevention, Fort Collins, CO, E-mail:uzu2/at/cdc.gov. Rebecca J. Eisen, Paul S. Mead, and Joseph Piesman, Bacterial Diseases Branch, Division of Vector-Borne Infectious Diseases, Centers for Disease Control and Prevention, Fort Collins, CO, E-mails:dyn2/at/cdc.gov,pfm0/at/cdc.gov, and jfp2/at/cdc.gov. Durland Fish and Maria A. Diuk-Wasser, Yale School of Public Health, New Haven, CT, E-mails: durland.fish/at/yale.edu and maria.diuk/at/yale.edu. Anne G. Hoen, Department of Community and Family Medicine, Dartmouth Medical School, Lebanon, NH, E-mail: anne.g.hoen/at/dartmouth.edu. Alan G. Barbour, Departments of Microbiology and Molecular Genetics and Medicine, University of California, Irvine, CA, E-mail: abarbour/at/uci.edu. Sarah Hamer, Department of Fisheries and Wildlife, Michigan State University, East Lansing, MI; and Department of Veterinary Integrative Biosciences, Texas A&M University, College Station, TX, E-mail: shamer/at/cvm.tamu.edu.
1. CDC Reported Cases of Lyme Disease by Year, United States, 1996–2010. 2011. http://www.cdc.gov/lyme/stats/chartstables/casesbyyear.html Available at. Accessed January 2011.
2. CDC Lyme Disease Incidence Rates by State, 2005–2010. 2011. http://www.cdc.gov/lyme/stats/chartstables/incidencebystate.html Available at. Accessed January 2011.
3. Connally NP, Durante AJ, Yousey-Hindes KM, Meek JI, Nelson RS, Heimer R. Peridomestic lyme disease prevention results of a population-based case-control study. Am J Prev Med. 2009;37:201–206. [PubMed]
4. Gould LH, Nelson RS, Griffith KS, Hayes EB, Piesman J, Mead PS, Cartter ML. Knowledge, attitudes, and behaviors regarding lyme disease prevention among Connecticut residents, 1999–2004. Vector Borne Zoonotic Dis. 2008;8:769–776. [PubMed]
5. Hayes EB, Piesman J. Current concepts—how can we prevent Lyme disease? N Engl J Med. 2003;348:2424–2430. [PubMed]
6. Piesman J, Eisen L. Prevention of tick-borne diseases. Annu Rev Entomol. 2008;53:323–343. [PubMed]
7. Shen AK, Mead PS, Beard CB. The lyme disease vaccine—a public health perspective. Clin Infect Dis. 2011;52:S247–S252. [PubMed]
8. Fish D. In: Ecology and Environmental Management of Lyme Disease. Ginsberg H, editor. New Brunswick, NJ: Rutgers University Press; 1993. pp. 25–42. (Population ecology of Ixodes dammini).
9. Mather TN, Nicholson MC, Donnelly EF, Matyas BT. Entomologic index for human risk of Lyme disease. Am J Epidemiol. 1996;144:1066–1069. [PubMed]
10. Stafford KC, Cartter ML, Magnarelli LA, Ertel SH, Mshar PA. Temporal correlations between tick abundance and prevalence of ticks infected with Borrelia burgdorferi and increasing incidence of Lyme disease. J Clin Microbiol. 1998;36:1240–1244. [PMC free article] [PubMed]
11. Connally NP, Ginsberg HS, Mather TN. Assessing peridomestic entomological factors as predictors for Lyme disease. J Vector Ecol. 2006;31:364–370. [PubMed]
12. Kitron U, Kazmierczak JJ. Spatial analysis of the distribution of Lyme disease in Wisconsin. Am J Epidemiol. 1997;145:558–566. [PubMed]
13. Falco RC, McKenna DF, Daniels TJ, Nadelman RB, Nowakowski J, Fish D, Wormser GP. Temporal relation between Ixodes scapularis abundance and risk for Lyme disease associated with erythema migrans. Am J Epidemiol. 1999;149:771–776. [PubMed]
14. Johnson JL, Ginsberg HS, Zhioua E, Whitworth UG, Markowski D, Hyland KE, Hu RJ. Passive tick surveillance, dog seropositivity, and incidence of human Lyme disease. Vector Borne Zoonotic Dis. 2004;4:137–142. [PubMed]
15. Rand PW, Lacombe EH, Dearborn R, Cahill B, Elias S, Lubelczyk CB, Beckett GA, Smith RP. Passive surveillance in Maine, an area emergent for tick-borne diseases. J Med Entomol. 2007;44:1118–1129. [PubMed]
16. Allan BF, Keesing F, Ostfeld RS. Effect of forest fragmentation on Lyme disease risk. Conserv Biol. 2003;17:267–272.
17. Brownstein JS, Skelly DK, Holford TR, Fish D. Forest fragmentation predicts local scale heterogeneity of Lyme disease risk. Oecologia. 2005;146:469–475. [PubMed]
18. Estrada-Pena A. Diluting the dilution effect: a spatial Lyme model provides evidence for the importance of habitat fragmentation with regard to the risk of infection. Geospat Health. 2009;3:143–155. [PubMed]
19. Ashley ST, Meentemeyer V. Climatic analysis of Lyme disease in the United States. Clim Res. 2004;27:177–187.
20. Brownstein JS, Holford TR, Fish D. A climate-based model predicts the spatial distribution of the Lyme disease vector Ixodes scapularis in the United States. Environ Health Perspect. 2003;111:1152–1157. [PMC free article] [PubMed]
21. Ogden NH, Maarouf A, Barker IK, Bigras-Poulin M, Lindsay LR, Morshed MG, O'Callaghan CJ, Ramay F, Waltner-Toews D, Charron DF. Climate change and the potential for range expansion of the Lyme disease vector Ixodes scapularis in Canada. Int J Parasitol. 2006;36:63–70. [PubMed]
22. Ostfeld RS, Canham CD, Oggenfuss K, Winchcombe RJ, Keesing F. Climate, deer, rodents, and acorns as determinants of variation in Lyme-disease risk. PLoS Biol. 2006;4:1058–1068. [PMC free article] [PubMed]
23. Gatewood AG, Liebman KA, Vourc'h G, Bunikis J, Hamer SA, Cortinas R, Melton F, Cislo P, Kitron U, Tsao J, Barbour AG, Fish D, Diuk-Wasser MA. Climate and tick seasonality are predictors of Borrelia burgdorferi genotype distribution. Appl Environ Microbiol. 2009;75:2476–2483. [PMC free article] [PubMed]
24. Dykhuizen DE, Brisson D, Sandigursky S, Wormser GP, Nowakowski J, Nadelman RB, Schwartz I. Short report: the propensity of different Borrelia burgdorferi sensu stricto genotypes to cause disseminated infections in humans. Am J Trop Med Hyg. 2008;78:806–810. [PMC free article] [PubMed]
25. Hanincova K, Ogden NH, Diuk-Wasser M, Pappas CJ, Iyer R, Fish D, Schwartz I, Kurtenbach K. Fitness variation of Borrelia burgdorferi sensu stricto strains in mice. Appl Environ Microbiol. 2008;74:153–157. [PMC free article] [PubMed]
26. Wormser GP, Brisson D, Liveris D, Hanincova K, Sandigursky S, Nowakowski J, Nadelman RB, Ludin S, Schwartz I. Borrelia burgdorferi genotype predicts the capacity for hematogenous dissemination during early Lyme disease. J Infect Dis. 2008;198:1358–1364. [PMC free article] [PubMed]
27. Diuk-Wasser MA, Gatewood AG, Cortinas MR, Yaremych-Hamer S, Tsao J, Kitron U, Hickling G, Brownstein JS, Walker E, Piesman J, Fish D. Spatiotemporal patterns of host-seeking Ixodes scapularis nymphs (Acari: Iodidae) in the United States. J Med Entomol. 2006;43:166–176. [PubMed]
28. Diuk-Wasser MA, Vourc'h G, Cislo P, Hoen AG, Melton F, Hamer SA, Rowland M, Cortinas R, Hickling GJ, Tsao JI, Barbour AG, Kitron U, Piesman J, Fish D. 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. 2010;19:504–514.
29. Diuk-Wasser M, Gatewood A, Cislo P, Brinkerhoff R, Hamer S, Rowland M, Cortinas R, Melton F, Hickling G, Tsao J, Bunikis J, Barbour A, Kitron U, Piesman J, Fish D. Human risk of infection with Borrelia burgdorferi, the Lyme disease agent, in eastern United States. Am J Trop Med Hyg. 2012;86:320–327. [PMC free article] [PubMed]
30. Bacon RM, Kugeler KJ, Mead PS. Surveillance for Lyme disease—United States, 1992–2006. MMWR Surveill Summ. 2008;57:1–9. [PubMed]
31. Guerra M, Walker E, Jones C, Paskewitz S, Cortinas MR, Stancil A, Beck L, Bobo M, Kitron U. Predicting the risk of Lyme disease: habitat suitability for Ixodes scapularis in the north central United States. Emerg Infect Dis. 2002;8:289–297. [PMC free article] [PubMed]
32. Wormser GP, Liveris D, Nowakowski J, Nadelman RB, Cavaliere LF, McKenna D, Holmgren D, Schwartz I. Association of specific subtypes of Borrelia burgdorferi with hematogenous dissemination in early Lyme disease. J Infect Dis. 1999;180:720–725. [PubMed]
33. Wang G, Ojaimi C, Iyer R, Saksenberg V, McClain SA, Wormser GP, Schwartz I. Impact of genotypic variation of Borrelia burgdorferi sensu stricto on kinetics of dissemination and severity of disease in C3H/HeJ mice. Infect Immun. 2001;69:4303–4312. [PMC free article] [PubMed]
34. Estrada-Pena A. Tick-borne pathogens, transmission rates and climate change. Front Biosci. 2009;14:U2674–U2780. [PubMed]
35. Tsao JI, Wootton JT, Bunikis J, Luna MG, Fish D, Barbour AG. An ecological approach to preventing human infection: vaccinating wild mouse reservoirs intervenes in the Lyme disease cycle. Proc Natl Acad Sci USA. 2004;101:18159–18164. [PubMed]
36. Hothorn T, Hornik K, Zeileis A. Unbiased recursive partitioning: a conditioxnal inference framework. J Comput Graph Statist. 2006;15:651–674.
37. Jones KL, Glickstein LJ, Damle N, Sikand VK, McHugh G, Steere AC. Borrelia burgdorferi genetic markers and disseminated disease in patients with early Lyme disease. J Clin Microbiol. 2006;44:4407–4413. [PMC free article] [PubMed]
38. Derdakova M, Dudioak V, Brei B, Brownstein JS, Schwartz I, Fish D. Interaction and transmission of two Borrelia burgdorferi sensu stricto strains in a tick-rodent maintenance system. Appl Environ Microbiol. 2004;70:6783–6788. [PMC free article] [PubMed]
39. Bunikis J, Garpmo U, Tsao J, Berglund J, Fish D, Barbour AG. Sequence typing reveals extensive strain diversity of the Lyme borreliosis agents Borrelia burgdorferi in North America and Borrelia afzelii in Europe. Microbiology. 2004;150:1741–1755. [PubMed]
40. Eisen L, Eisen RJ. Need for improved methods to collect and present spatial epidemiologic data for vectorborne diseases. Emerg Infect Dis. 2007;13:1816–1820. [PMC free article] [PubMed]
41. Goddard J, Piesman J. New records of immature Ixodes scapularis from Mississippi. J Vector Ecol. 2006;31:421–422. [PubMed]
Articles from The American Journal of Tropical Medicine and Hygiene are provided here courtesy of
The American Society of Tropical Medicine and Hygiene