An Acarological Risk Model Predicting the Density and Distribution of Host-Seeking Ixodes scapularis Nymphs in Minnesota.

Ixodes scapularis is the vector of at least seven human pathogens in Minnesota, two of which are known to cause Lyme disease (Borrelia burgdorferi sensu stricto and Borrelia mayonii). In Minnesota, the statewide incidence of Lyme disease and other I. scapularis-borne diseases and the geographic extent over which cases have been reported have both increased substantially over the last two decades. These changes correspond with an expanding distribution of I. scapularis over a similar time frame. Because the risk of exposure to I. scapularis-borne pathogens is likely related to the number of ticks encountered, we developed an acarological risk model predicting the density of host-seeking I. scapularis nymphs (DON) in Minnesota. The model was informed by sampling 81 sites located in 42 counties in Minnesota. Two main foci were predicted by the model to support elevated densities of host-seeking I. scapularis nymphs, which included the seven-county Minneapolis-St. Paul metropolitan area and counties in northern Minnesota, including Lake of the Woods and Koochiching counties. There was substantial heterogeneity observed in predicted DON across the state at the county scale; however, counties classified as high risk for I. scapularis-borne diseases and counties with known established populations of I. scapularis had the highest proportion of the county predicted as suitable for host-seeking nymphs (≥ 0.13 nymphs/100 m2). The model provides insight into areas of potential I. scapularis population expansion and identifies focal areas of predicted suitable habitat within counties where the incidence of I. scapularis-borne diseases has been historically low.

INTRODUCTION

In the United States, Lyme disease, caused by infection with Borrelia burgdorferi sensu stricto or less commonly by the newly discovered organism Borrelia mayonii, is the most commonly reported vector-borne disease.[1,2] An average of 30,000 Lyme disease cases are reported each year, representing a 3-fold increase from 1992 to 2013, and most Lyme disease cases are reported from the Northeast, mid-Atlantic, and upper Midwestern United States.[1] Since the mid-1990s, the reported distribution of Ixodes scapularis, the primary vector of Lyme disease spirochetes, has expanded in the United States with the greatest increases observed in the upper Midwest and in the northeastern states.[3] In parallel to the observed expansion of the vector’s geographic range, the numbers of counties considered high incidence for Lyme disease by Kugeler et al.[4] also expanded in the same geographic regions over a similar time period. Some of the most dramatic increases in the occurrence of Lyme disease and other I. scapularis–borne illnesses have been observed in Minnesota, where, from 1992 to 2006, the percentage of counties reporting at least one human case of Lyme disease increased from 33% to 74%.[5] Moreover, from 1996 to 2011, there was a substantial and significant increase in the incidence of Lyme disease, anaplasmosis, and babesiosis.[6] Along with a statewide increase in the incidence of I. scapularis–borne diseases, counties reporting human cases of I. scapularis–borne diseases[6] and established populations of I. scapularis have expanded geographically.[3]

In Minnesota, I. scapularis serves as the vector of at least seven human pathogens including B. burgdorferi sensu stricto (henceforth referred to as B. burgdorferi), Anaplasma phagocytophilum, Babesia microti, Powassan virus, Borrelia miyamotoi, and two novel I. scapularis–borne pathogens (B. mayonii and Ehrlichia muris subs. eauclarensis) both of which appear to be restricted to the upper Midwest.[2,7-14] To more clearly define the current distribution of areas considered suitable for the establishment of I. scapularis populations in Minnesota, we developed a habitat suitability model for I. scapularis in Minnesota.[7] Nymphal ticks pose the greatest risk to humans because they are difficult to detect owing to their small size, and the seasonal timing of peak nymphal questing activity coincides with a period of increased human outdoor activity.[1,15-17] However, recognizing that 1) human risk of exposure to pathogens has been correlated with density of host-seeking nymphal I. scapularis, rather than simple measures of presence or absence,[18-20] and 2) habitat suitability for ticks is heterogeneous within counties,[21] we conducted an extensive field survey of host-seeking I. scapularis nymphs and developed a subcounty (30 × 30 m) resolution model of the DON to better inform estimates of acarological risk of human exposure to I. scapularis.

METHODS

Site selection for tick sampling.

We recently modeled the distribution of suitable habitat for I. scapularis in Minnesota at a subcounty scale (30 × 30 m) using the machine learning ecological niche model Maxent, which was fit using tick occurrence records previously collected by the Minnesota Department of Health.[7] Suitable areas in that study were best described by seven variables: land cover (especially cool temperate forests), maximum temperature during the warmest month, the amount of precipitation during the wettest quarter, the annual temperature range, the mean diurnal temperature range, elevation, and the amount of precipitation during the coldest quarter of the year. To estimate the density of host-seeking nymphs within a suitable habitat for the tick and to refine estimates of the suitable habitat based on a larger sample size, in this study, we directed tick sampling efforts to areas considered suitable for tick survival based on our previous model.[7]

To be considered for inclusion as a sampling site for the present study, a site had to be 1) classified as suitable habitat for I. scapularis by the Maxent model,[7] 2) publicly owned, 3) larger than 500 × 500 m, 4) accessible from maintained public access roads, and 5) have suitable vegetation, for example, contiguous forest with at least 50% closed canopy, as assessed by reviewing recent aerial imagery (GoogleEarth v 7.1.8.3036). These criteria were intended to ensure a reasonable expectation of finding ticks if they are present, and adequate area that was accessible to sample. To extrapolate model findings across Minnesota while minimizing instances where we were making predictions about areas with environmental characteristics that were not included in our build set, in addition to the five inclusion criteria already described, we selected sites with climatic and environmental variables characteristic of the rest of the state (Table 1) and also sought vast geographic coverage of the state. To accomplish the latter, we specified that all potential sampling sites had to be separated by at least 10 km. To accomplish the former, we sampled sites by using a scheme based on the concept of conditional Latin hypercube sampling (cLHS).

Table 1

Range of each variable used for predicting the distribution of suitable habitat

Variable	Statewide range	Maxent range	Selected range
Mean diurnal temperature range (°C × 10)	98–135	99–133	110–132
Maximum temperature during the warmest month (°C × 10)	218–301	230–296	247–290
Annual temperature range (°C × 10)	402–505	402–495	424–490
Mean temperature of the coldest quarter (°C × 10)	−158 to −70	−143 to −70	−100 to −83
Precipitation during the wettest quarter (mm)	214–336	246–336	278–327
Precipitation during the coldest quarter (mm)	39–116	43–116	51–81
Elevation (m)	175–668	179–600	227–525

Ranges are shown for the entire state of Minnesota, for the entire distribution of the Maxent model, and the range at the points selected for sampling in 2015.

GAP = U.S. Geological Survey Gap Analysis Program. May 2011, National Land Cover version 2.

The cLHS method works by maximally stratifying the distribution of each variable into as many strata as sampling sites desired. The method allows for efficient sampling of existing data and has been shown to be the most effective way to replicate the distribution of variables, proving superior to random sampling, principal components, and equal spatial strata sampling schemes.[22,23] Using cLHS, or similar approaches, results in a sample set that covers the full range of conditions encountered among the design variables. This minimizes the potential for extrapolation to a set of predictors not used to fit the model. Conditional Latin hypercube sampling has been used extensively to guide sampling site selection for soil and vegetation modeling.[24-28]

To identify the range of each variable across suitable habitats in Minnesota (Table 1), we resampled each of the seven variables used in our original suitability model[7] to a resolution of 90 × 90 m to create a manageable number of data points. These variables represent several key predictors of the density of host-seeking Ixodes spp. nymphs from other areas and include measures of forest type or coverage, temperature, precipitation, and elevation.[19,29-40] Moreover, these seven variables are not strongly correlated with each other, but are strongly correlated with several other potential predictive variables, as described previously.[7] The value of each of these seven variables was extracted to each 90 × 90 m pixel falling within the distribution of Maxent-predicted suitable habitat meeting the aforementioned selection criteria using the Extract Multiple Values to Points tool in ArcGIS (ArcMap 10.2; ESRI, Redlands, CA).

The site selection algorithm used in this study is not strict cLHS as we sample from the marginal distributions of the covariates and incorporate a proximity constraint. These modifications yield sample sites that span the domains of the covariates, individually, and provide more complete geographic coverage. We divided the distribution of each of the seven variables into five classes or bins based on quantiles, except land cover, which was binned according to the proportion of each land cover type predicted by the Maxent model; for example, 67.1% of binned values were cool temperate forests and 15.7% lowland and montane boreal forest,[7] using the cut2 package in R (R v 3.2.1; R Foundation for Statistical Computing, Vienna, Austria). We then randomly chose two points (potential sampling locations) without replacement from each bin. To minimize geographic clustering of sampling sites, we limited the proximity of selected points to 10 km from the nearest selected point by placing a 10-km buffer around a selected point and removing all other points falling within the buffer. The order of selection was based on the importance of each variable in the Maxent model as listed previously. Because of the potential for some selected sites to be inaccessible as a result of unforeseen circumstances such as flooding or closed roads, we chose a second set of sampling sites, here referred to as backup sites. Backup sites were chosen from the points that remained after the first dataset was chosen and points within the 10-km buffer had been removed from the dataset. Eighty primary and 80 secondary sampling sites were chosen, yielding a total of 160 candidate sampling sites. Each of the 80 primary sites was visited in April 2015, before sampling to ensure habitat suitability, that is, wooded and brushy mesic areas with at least 50% canopy coverage. If the primary site was found to be unsuitable or inaccessible, the secondary site was assessed for suitability.

Random site selection using the modified cLHS resulted in a distribution of 80 points located in 41 counties. We ultimately sought to predict the distribution and density across the state; however, potential sampling sites identified by the cLHS were extremely limited in the north, with only a single sampling site chosen in far northeastern Minnesota. To increase the spatial extent of sampled sites, we chose to include an additional sampling site at Voyageurs National Park. The site sampled at Voyageurs National Park met all of the aforementioned inclusion criteria. Ultimately, 81 sites located in 42 counties were chosen for sampling (Figure 1). All initial site selection was performed in R (R v 3.2.1).

Figure 1.

Locations of 81 study sites sampled for ticks in Minnesota shown overlaid with county-level Ixodes scapularis population status (this research yielded two counties with new established populations [Benton County and Sibley County] and two counties with new records of I. scapularis [Rice County and Nicollet County]; Eisen et al.[3]) and county-level I. scapularis–borne disease risk were defined as the number of cases per 100,000 population (high risk: ≥ 25 cases, moderate risk: 10–24.9 cases, low risk: < 10 cases) (Minnesota Department of Health, http://www.health.state.mn.us/divs/idepc/diseases/lyme/highrisk.html). The inset map shows the 2015 human population by county.

Tick drag sampling.

Eighty-one sites were sampled on two occasions, each separated by a range of 6–22 days, from May 31 to June 30, 2015. This sampling period was timed to coincide with the predicted peak of I. scapularis nymphal questing activity and just before the reported onset of illness for most Lyme disease cases in Minnesota.[41-43] The density of host-seeking I. scapularis nymphs (DON) was estimated by drag sampling 750 m2 of forest floor with a 1-m2 tick drag made of rubber-bonded cotton sheet (JoAnn Fabric no. 1491315) with a rope attached to a 48″ dowel inside the top edge. Weighted “fingers” were sewn to the bottom half of the drag to ensure sampling occurred near the forest floor and leaf litter layer. To minimize the likelihood of ticks falling off the drag before being collected, samplers stopped every 15 m, removed all ticks from the drag and themselves, and placed them in prelabeled vials containing 70% ethanol. All ticks were sent to the Centers for Disease Control and Prevention, Fort Collins, CO, identified to species[44-46] and life stage, and tested for pathogens.

Predictive variables for modeling the number of ticks.

The variables in the Maxent model informing site selection were based on average monthly climate data for 1960–1990.[7,47] We sought to produce a risk model which included current climate conditions that may have affected nymphal tick collections in 2015 and, therefore, chose to use 35-year climate averages (1980–2014) from version 2 of the 1-km spatial resolution Daymet dataset[48,49] to derive the BIOCLIM variables used to select sampling sites. This likely had minimal effect on the model as common variables among the two datasets were highly correlated (Pearson’s r ≥ 0.80; data not shown). In addition to the 19 BIOCLIM variables, we also considered eight other climate variables, including quarterly and annual measures of vapor pressure, growing degree days, soil water depletion, and annual precipitation accumulation (Table 2).

Table 2

Climate and landscape variables included for initial consideration in the zero-inflated negative binomial (ZINB) model

ZINB model variable	Correlated variable(s)
Mean diurnal range (mean of monthly (max temp − min temp))	Average vapor pressure (summer)
	Average vapor pressure (annual)
Temperature annual range	Temperature seasonality (standard deviation × 100)
	Mean temperature of the driest quarter
	Min temperature of the coldest month
	Precipitation of the driest month
	Precipitation seasonality (coefficient of variation)
	Precipitation of the driest quarter
	Precipitation of the coldest quarter
	Annual precipitation accumulation
	Average vapor pressure (annual)
	Average vapor pressure (winter)
	Average vapor pressure (fall)
Mean temperature of wettest quarter	Annual mean temperature
	Mean temperature of the driest quarter
	Average vapor pressure (summer)
	Min temperature of the coldest month
	Precipitation of the driest month
	Growing degree days (1 January–30 June)
	Soil water depletion (annual)
	Average vapor pressure (annual)
	Average vapor pressure (winter)
	Average vapor pressure (spring)
	Average vapor pressure (fall)
Precipitation of wettest quarter	Annual precipitation
	Precipitation of the wettest month
	Precipitation of the warmest quarter
	Annual precipitation accumulation
	Average vapor pressure (spring)
	Average vapor pressure (fall)
Elevation*	Isothermality
	Max temperature of the warmest month
	Mean temperature of the coldest quarter
Distance (m) to the nearest stream or river corridor (m)
Percent of agricultural land cover within 5.25 km2 buffer
Percent of forest land cover within 5.25 km2 buffer

Pearson’s correlation coefficient was calculated for all pairs of variables and highly correlated variables (correlation coefficient > 0.75) were represented by a single variable in the model.

National Elevation Dataset; https://lta.cr.usgs.gov/NED, last visited October 2014.

In addition to the eight climate candidate variables, we included measures of vegetation, elevation, and distance to streams which have been shown to be associated with tick abundance.[19,50,51] All sites chosen for sampling were located in closed canopy deciduous forest; thus, to simplify vegetation into categories, we chose two general measures of land cover surrounding each study site. The percentage of cool temperate forest and the percentage of agricultural land within a 5.25-km2 square buffer were calculated for each sampling location as well as the percentage of each land cover type within each 5.25 km2 pixel across the state. Because white-tailed deer serve as the primary host for adult ticks and are likely a considerable dispersal mechanism for ticks to inhabit new areas, the buffer area of 5.25 km2 was chosen to represent the home range of a male white-tailed deer (Odocoileus virginianus Zimmerman) in Minnesota (Minnesota Department of Natural Resources; www.dnr.state.mn.us/index.html). We calculated the distance to a stream or river corridor (http://www.arcgis.com, last visited July 2015) for each site and for each pixel in the state using the Near tool in ArcMap (v10.3; ESRI), which calculates the distance from each study site or centroid of each pixel to the nearest river or stream. For example, if a stream or river runs through a pixel, then the pixel value is zero. We extracted values for the variables for each of the 81 sampling sites using the aforementioned Extract Multiple Values to Points tool. We identified correlations among the 31 variables, that is, Pearson’s correlation coefficient ≥ 0.75 (Table 2), and ultimately included eight variables, four of which were also represented in the habitat suitability model.[7] All ecological data were preprocessed using ArcMap (v10.3; ESRI). Preprocessing included projecting all layers to Albers Equal Area North American Datum 1983, resampling all layers to a resolution of 30 × 30 m, and extracting the values of variables for each sampling site.

Model development.

We considered both the nature of our site selection and our data structure when choosing a modeling framework. All sites selected for tick sampling were located within the predicted distribution of suitability (see Johnson et al.[7]) and visually confirmed to have wooded and brushy mesic areas with at least 50% canopy coverage in April 2015. We observed zero ticks at 17 of 81 sites (21%), more than would be expected under a Poisson distribution. For sites where only a single nymph was collected per 750 m2 sampled, DON was equal to 0.13 nymphs/100 m2; therefore, we defined the lowest risk category as < 0.13 nymphs/100 m2 (Figure 2). To represent the highest potential risk of encountering nymphal I. scapularis in a given area, we identified the single sampling visit that yielded the highest number of nymphs and used this number of ticks to create a predictive model of the peak number of ticks likely encountered per 750 m2. Overdispersion, relative to a Poisson distribution, in the data (variance > mean) was tested for and identified using the AER package in R (R v 3.2.1; [dispersiontest: P = 0.001]).

Figure 2.

Predicted density of I. scapularis nymphs/100 m2 (DON) derived from the number of nymphs predicted by the zero-inflated negative binomial (ZINB) model compared with observed densities from field collections. The inset shows the standard deviation of the predicted density of I. scapularis nymphs/100 m2 and the black dots show the 11 study sites where no nymphs were collected but where the model predicted low densities of nymphs. This figure appears in color at www.ajtmh.org.

To account for excess zeros and overdispersion in the data, a zero-inflated negative binomial (ZINB) model was fit to the peak number of nymphs observed. As a mixture model, the ZINB model incorporates two sources of zeros: structural zeros, which represent sites where ticks are truly absent, and random zeros, which represents low-density sites where, on a given day, the observed count was zero when sampling occurred. The ZINB comprises a discrete, nonnegative distribution describing tick counts (Equation S1) and a binary distribution describing the source of true zeros or absences in the dataset (Equation S2). The two portions of the model are then combined to estimate the number of ticks expected at each site based on the aforementioned ecological variables (Equation S3). To ensure the ZINB model represented the best fit to the data, we used Vuong statistic[52] to determine if the ZINB distribution fit the data better than a negative binomial distribution.

To select variables to include in the final model, we first fit a full model containing all possible variables represented in both portions of the ZINB model (Table 2) using package pscl in R (R v 3.2.1). We then fit a reduced model including variables contributing significantly to model fit as identified by backward stepwise regression. We tested for spatial autocorrelation of observed DON/100 m2 and among model residuals given by the best fit reduced ZINB model using Moran’s I (Spatial Analyst, ArcMap v 10.3; ESRI). Model performance (fit) was assessed by examining the scaled Pearson’s residuals and using 5-fold cross-validation (CV); CV was performed 10 times and the mean absolute prediction error was computed. Examination of the distribution of residuals allowed us to further evaluate the model and identify outliers in the predicted number of ticks among the 81 sites. We used the parameter estimates given by the ZINB model and Raster Calculator (ArcMap v10.3; ESRI) to predict the number of nymphs (Equation S3) and evaluated statewide model fit using the delta method to estimate variance (Equation S4). Based on the predicted number of nymphs per site (750 m2), we calculated and produced a continuous risk map of the DON/100 m2 by dividing the number of nymphs predicted by 7.5. To present a more realistic distribution, we masked out areas deemed uninhabitable by ticks, for example, pixels defined as open water, swamp, agriculture fields, and high-intensity developed areas.

Voyageurs National Park was not selected using the cLHS framework but was included to extend the northern extent of sampling sites in central Minnesota. To verify that the inclusion of this point did not influence the predicted DON, we reran the ZINB model without that site and then compared parameter estimates and slope between models with or without the site included. Predictive maps were compared between those developed with versus without Voyageurs National Park using Raster Calculator (ArcMap v 10.3; ESRI).

Assessment of acarological risk model relative to human I. scapularis–borne disease risk and I. scapularis distribution records.

We assessed the agreement between DON, as estimated by the ZINB model with Lyme disease incidence from 2000 to 2015 (http://www.cdc.gov/lyme/stats/), and categorical county-level tick-borne disease risk based on the average incidence of reported I. scapularis–borne disease cases, including Lyme disease, anaplasmosis, and babesiosis from 2007 through 2015. Risk categories were defined as cases of I. scapularis–borne disease per 100,000 population: low risk: < 10 cases, moderate risk: 10–24.9 cases, and high risk: ≥ 25 cases, to effectively create a risk map of the state that generally incorporated the potential exposure risk to I. scapularis based on the natural biomes of Minnesota (http://www.health.state.mn.us/divs/idepc/diseases/lyme/highrisk.html). Within these categories, 32 counties are classified by the Minnesota Department of Health as high risk, 22 counties are classified as moderate risk, and the remaining 33 counties are classified as low risk (Figure 1). In addition to epidemiological outcomes, model-predicted DON was compared with previously described county-level status of I. scapularis populations, that is, no records, reported, and established.[3]

For each county, we calculated the percent area classified in each of three categories: 1) predicted risk (≥ 0.13 nymphs/100 m2), 2) predicted low risk (< 0.13 nymphs/100 m2), and 3) unsuitable habitat (masked out), using the Tabulate Area Spatial Analyst tool (ArcMarp v 10.3; ESRI). Differences in percent area predicted as risk and county-level tick-borne disease risk and county-level tick establishment were analyzed using the nonparametric Kruskal–Wallis rank sums test with the Dunn method for correction of multiple comparisons using the PMCMR package in R. The association of percent area predicted as risk with county-level Lyme disease incidence was assessed using linear regression and controlled for county population size. All statistical tests were carried out at a significance level of α = 0.05 performed in R (R v 3.2.1).

RESULTS

Tick collections.

Sampling sites were located in 42 of 87 Minnesota counties. Most sampling sites (70%; 57/81) were located in counties classified as high risk for I. scapularis–borne disease, 20% (16/81) of sites were located in counties classified as moderate risk, and the remaining 10% (8/81) were located in counties classified as low risk (Figure 1). Using the definitions of established (≥ 6 ticks, or two life stages) and reported (presence of at least one tick) populations described by Dennis et al.,[53] established populations of I. scapularis were identified at 80% of the sites sampled. Nine sampling sites were located in counties where previous records of I. scapularis were lacking.[3] From these nine sites, we documented new records of established populations in two counties, Benton and Sibley, and new records of reported populations in two counties, Rice and Nicollet (Figure 1).

Among the 81 sites sampled between May 31 and June 30, 2015, a total of 1,386 I. scapularis nymphs were collected from 64 of the 81 (79%) sites sampled; 35% of all I. scapularis collected were nymphs. The number of nymphs collected per site ranged from 0 to 77 on a single visit, with a median of five nymphs (0.66 nymphs per 100/m2). The largest number of nymphs collected during a single sampling session per site was used to calculate the DON for each site; DON ranged from zero to as high as 10.3 nymphs/100 m2. The highest DON was observed near the Minneapolis-St. Paul metropolitan area, and the geographic distribution of other high-density sites trended in a northwest direction away from the metropolitan area (Figure 2). Incidental captures included 831 adult I. scapularis from 66 sites, at least 1,721 I. scapularis larvae from 53 sites, 1,236 Dermacentor variabilis adults and two nymphs from 71 and two sites, respectively, and four Haemaphysalis leporispalustris nymphs from three sites. One male Ixodes marxi was collected from Otter Tail County and a single Ixodes texanus nymph was collected in Wright County, which is a new record of this species in Minnesota (L. Beati, Personal Communication; National Tick Collection, Georgia Southern University).

Zero-inflated negative binomial model of DON.

The zero-inflated (predicted absence) portion of our ZINB model showed a strong positive association with tick absence as the proportion of agricultural land increased (Table 3). More intuitively, the probability of tick presence increased as the proportion of agricultural land decreased (Figures 2 and 3). Because of the distribution of agricultural land in Minnesota, large portions of western and southern Minnesota were classified as unsuitable (Figures 2 and 3). The ZINB correctly predicted the presence of I. scapularis nymphs with high accuracy (79%). Most of the sites where the model predicted suitability but ticks were not collected were in areas where the model predicted low abundance of host-seeking nymphs (Figure 2).

Table 3

Parameter estimates for variables included in the zero-inflated negative binomial model

	Estimate	Standard error	z-value	PR > |z|
Zero-inflated
Intercept	−4.37	1.07	−4.08	< 0.001
% Agricultural land	7.24	2.05	3.52	< 0.001
Negative binomial
Intercept*	−27.92	7.98	−3.50	< 0.001
% Agricultural land	−3.32	0.67	−4.97	< 0.001
Mean diurnal temperature range*	−1.60	0.41	−3.95	< 0.001
Annual temperature range*	0.91	0.20	4.43	< 0.001
Precipitation of the wettest quarter	0.03	0.01	4.73	< 0.001
Elevation	< −0.01	< 0.001	−2.64	< 0.001

The zero-inflated portion of the model predicts the probability of no ticks (absence site); the negative binomial portion predicts the number of ticks.

When the model was run excluding Voyageurs National Park, parameter estimates and standard errors remained unchanged for all variables except the following, which are listed with revised parameter estimates (standard error) shown: Intercept −26.33 (7.97), mean diurnal temperature range −1.57 (0.41), and annual temperature range 0.86 (0.23).

Figure 3.

The distribution of each of the variables contributing significantly to the zero-inflated negative binomial (ZINB) model for predicting the number of I. scapularis nymphs/750 m2. The color ramps match the effect of the variables in the ZINB model (Table 3). That is, darker colors of each variable indicate association with higher tick counts as predicted by the ZINB model.

Five variables contributed significantly to the negative binomial (number of nymphs) portion of the ZINB model (Table 3). Higher amounts of agricultural land in the surrounding area had the strongest negative impact on the predicted number of nymphs. Mean diurnal temperature range and elevation also had a negative association with increasing numbers of nymphs (although their influence was less than the percentage of agricultural land), whereas annual temperature range and the amount of summer precipitation were positively associated with increasing numbers of nymphs (Table 3). Mean diurnal temperature range was highly correlated with annual and summer vapor pressure, and precipitation during the wettest quarter was highly positively correlated with annual mean temperature, annual and summer precipitation, and spring and fall vapor pressure (Table 2). Annual temperature range was significantly correlated (Pearson’s r > 0.75) with measures of extremes in temperature and precipitation, including minimum temperature of the coldest quarter, mean temperature and precipitation during summer and winter months, and annual precipitation (Table 2).

Examination of Pearson’s χ2 residuals produced by the ZINB model showed that all absolute values were less than three and only four sites had absolute values greater than two, fewer than would be expected by chance. Pearson’s scaled residuals > |2| represent sites where the predicted number of ticks is higher or lower than would be expected, which may lead to increased uncertainty of these predictions. Positive residuals greater than two indicate that the model is underestimating the number of nymphs at these four sites (Figure 2). Despite sites with large residuals being located near each other, no spatial autocorrelation was observed among the ZINB scaled residuals (Moran’s I: z-score = −1.01, P = 0.31) nor was spatial autocorrelation noted in observed DON (Moran’s I: z-score = 0.82, P = 0.41). That is, sites located in close proximity to one another were not significantly more likely to have similar DON. Five-fold CV of the model consistently predicted DON within 1.3 nymphs/100 m2 and the mean standard deviation across the state was 0.93 nymphs/100 m2 (Figure 2). As would be expected for a negative binomial distribution, the variance estimated using the delta method shows the highest variance in the metropolitan area and in Lake of the Woods and Koochiching counties in northern Minnesota, where the predicted number of nymphs was high.

Statewide predicted DON ranged from 0 to 19.6 nymphs/100 m2, with the highest densities predicted in Hennepin, Anoka, Ramsey, and Dakota counties (within the Minneapolis-St. Paul metropolitan area) and in northern Minnesota in St. Louis, Koochiching, Lake of the Woods, and Roseau counties (Figure 2). Overall, 31% of the state was predicted to support peak host-seeking nymphal densities of at least 0.13 nymphs/100 m2 (Figure 2). The inclusion of Voyageurs National Park had minimal effects on the ZINB model (Table 3). The predicted distributions were examined by subtracting the distribution of the model containing Voyageurs National Park from the distribution produced after removing Voyageurs National Park. All estimates of DON were slightly higher when Voyageurs National Park was removed. However, > 99% of pixels in the state differed by less than 1 nymph/100 m2. Among the remaining pixels comprising < 1% of the state, the maximum difference at any given location in the predicted nymphal density between the model with and without Voyageurs National Park was 1.6 nymphs/100 m2.

Comparison of DON with county-level epidemiologic and entomologic data.

The percent area predicted to have ≥ 0.13 nymphs/100 m2 in each county was significantly higher in “high” I. scapularis–borne disease incidence counties (median area predicted as risk per county = 43%) than “moderate” incidence counties (median area predicted as risk per county = 11%) and “low” incidence counties (median area predicted as risk per county < 0.5%) (Kruskal–Wallis rank sums test, P < 0.02; Figure 4). In addition, the percent area predicted as risk in each county was significantly higher for counties with established populations of I. scapularis (median area predicted as risk per county = 39%) than counties with reported populations (median area predicted as risk per county = 5%) or no records of I. scapularis (median area predicted as risk per county = 1%) (Kruskal–Wallis rank sums test, P < 0.001; Figure 4). There was no significant association between the percent area predicted as risk in each county and the cumulative number of Lyme disease cases reported per county from 2000 to 2015 when county population was accounted for (linear regression, P = 0.11).

Figure 4.

Comparison of the percent area predicted as risk in each county with (A) I. scapularis county-level population status (Eisen et al. 2016) and (B) county-level risk of I. scapularis–borne disease given as number of cases per 100,000 population (high risk: ≥ 25 cases, moderate risk: 10–24.9 cases, low risk: < 10 cases) (Minnesota Department of Health; http://www.health.state.mn.us/divs/idepc/diseases/lyme/highrisk.html). Groups with the same alphabetical designation are not significantly different (P ≥ 0.01).

DISCUSSION

In our model, DONs in Minnesota are predicted primarily by the amount of agricultural land in close proximity to the location of interest, a variable that is inversely related to the amount of forest. However, extreme temperatures and measures of precipitation are also identified as significant predictors. Based on our acarological risk model, we predict that ecological and climatic conditions are suitable to support peak densities of at least 0.13 host-seeking I. scapularis nymphs per 100 m2 in approximately one-third of the state. The highest densities are predicted to occur in the Minneapolis-St. Paul metropolitan area and in far north-central counties bordering Canada (Lake of the Woods, Koochiching, and St. Louis), with moderate density areas spanning between these foci. Large portions of western and southern Minnesota are considered unsuitable. Our predictions are consistent with epidemiologic trends of reported I. scapularis–borne diseases in Minnesota. Specifically, the percentage of the county classified by our model as suitable for supporting at least 0.13 nymphs/100 m2 is lowest in counties with low incidence of I. scapularis–borne diseases and incrementally greater in counties with moderate and high incidence of I. scapularis–borne diseases[6] (http://www.health.state.mn.us/divs/idepc/diseases/lyme/highrisk.html, last visited April 2017). Low-risk counties with high percent risk area predicted include Rice, Le Sueur, Sibley, Pope, and McLeod in descending order. In addition, counties with no records of I. scapularis but with high percent elevated risk area predicted include Roseau, Le Sueur, Meeker, and McLeod in descending order.

Despite nuanced differences that likely arise from the use of different tick records, predictive variables, and/or modeling approaches, our model predictions are largely consistent with previous assessments of acarological measures of the spatial distribution of I. scapularis in Minnesota. Specifically, our model predicts higher nymphal densities for counties with established I. scapularis populations compared with those where tick records were lacking[3,7] and follows spatial trends similar to previous predictive tick distribution models developed for Minnesota[7] or the eastern United States including Minnesota.[19,29,37]

Counties in central Minnesota and eastern counties bordering Wisconsin are consistently predicted by our model and previously published models to harbor I. scapularis[7,19,29,37] and correspond with areas of high incidence for I. scapularis–borne diseases.[6] Likewise, in agreement with our model, each of these models classified western Minnesota as unlikely to harbor I. scapularis and risk for I. scapularis–borne diseases is uniformly low in the western and southwestern parts of the state (Robinson et al.[6], http://www.health.state.mn.us/divs/idepc/diseases/lyme/highrisk.html, last visited April 2017) where most models predict lower habitat suitability and low DON[7,19,37] or no I. scapularis[29] (current model). Several counties in the northwestern portion of the state are classified as moderate or high risk. Notably, our model predicts some focal areas to be suitable within most of these counties and could represent areas of potential exposure. However, it is also possible that many cases reported from these counties were not exposed within their county of residence. Indeed, many of these counties are in close proximity to recreational areas classified as suitable for supporting very high densities of I. scapularis nymphs. Despite the success in collecting ticks at the Voyageurs National Park site, which was located in the northeast corner of Koochiching County, and the high density of I. scapularis predicted by the model, much of the rest of Koochiching county is swampy and rich in peat deposits, which is likely an unsuitable habitat for I. scapularis populations.

Our model identified a limited number of focal areas in southwestern Minnesota to be suitable for I. scapularis, a finding consistent with most previous models that considered this region to be low risk.[19,29] Southern Minnesota is characterized by vast expanses of agricultural lands and few forests (except riparian woodlands and planted shelter belts), resulting in the classification of this area by our model as unsuitable (Figure 3). However, a previous model for Minnesota[7] identified a potentially more expansive range of predicted habitat suitability in southeastern Minnesota where cool temperate forests are interspersed among agricultural lands. Notably, most of the sites from which we were unable to find ticks, despite the model classifying them as suitable, are located in the southern half of the state, generally in forested areas fragmented by agriculture. The incidence of I. scapularis–borne diseases is high in southeastern counties but typically lowest in the south-central and southwestern portions of the state (Figure 1).[6] As the tick continues to expand its range in Minnesota, many of these isolated forested patches may become colonized or tick densities might increase to detectable levels.

Previous models were inconsistent in predictions for the northeast and far northern counties. Correspondingly, uncertainty in our model is highest in the far northern counties and underscores a need for additional vector surveillance in that region. These differences may simply reflect differences in the input data used to develop the models. Specifically, the earliest species distribution models[29,37] were based on county-level tick records from the mid-1990s,[53] before I. scapularis became established in northern Minnesota. Lack of records from the north might explain why the models did not predict suitable habitat in far northern counties in Minnesota. Consistent with our model, Diuk-Wasser et al.[19] predicted extensive suitable tick habitat in north and northeastern Minnesota. This model was based on tick data collected in 2004 and used the same modeling framework (ZINB) as our study. Notably, our models predict that northern Minnesota is suitable to support much higher nymphal densities than the model presented by Diuk-Wasser et al.[19] This is likely attributable to differences in observed nymphal densities. Based on drag sampling, we observed nymphal densities up to 10 times greater than that observed by Diuk-Wasser et al.[19] when they sampled between 2004 and 2006. This difference may reflect more time since establishment and is consistent with increasing incidence of I. scapularis–borne disease cases.[6] Alternatively, these differences could simply reflect interannual variation in tick abundance and, thus, emphasize a limitation of our study: existing model predictions are based on a single year of drag sampling. Areas predicted to be suitable but where zero nymphs were observed may represent areas of future colonization by the ticks. Many of these areas are at the agricultural–forest intersection which may delay colonization because of increased fragmentation by agricultural land and less contiguous forest habitat, which may impede host movement to other suitable habitat patches.[54]

Our model is biologically plausible. The proportion of agricultural land, a variable inversely related to forest cover, was a strong significant predictor of the absence of ticks and was the most influential predictor of the number of nymphs. This finding is consistent with other studies that found measures of forest cover to be important predictors for tick distribution and abundance.[3,19,29,30]

Minimum and maximum temperatures[19,29,37] and measures of precipitation or humidity have been shown to be strong predictors of I. scapularis distributions.[19,29,37] Here, three measures of climate variability contributed significantly to model fit. Specifically, our model indicated a need for sufficient amounts of precipitation during the summer and moderate temperature fluctuation when I. scapularis nymphs are actively host-seeking. Lower diurnal temperature range, higher summer precipitation, and higher annual temperature variation were predicted to increase the number of nymphs. A negative association was found between the number of nymphs and the mean diurnal temperature range, suggesting that DON may be higher in areas where the daily temperature remains relatively stable. These areas are consistent with areas of higher humidity, which is important for maintaining water balance in ticks and preventing desiccation[55] and is essential for tick survival and host-seeking.[55] Furthermore, annual temperature range was positively correlated with predicted nymphal density. Areas with a broad range in annual temperature experience warmer summer temperatures and lower winter temperatures. Field studies have demonstrated that I. scapularis can survive subzero winter temperatures,[38,40] presumably because of an insulating effect of the leaf litter that is frequently associated with woodlands and the insulating effect of snow cover. Therefore, the positive association we observed between annual temperature range and predicted nymphal abundance likely reflects warmer summer temperatures favoring tick survival and reproduction. In general, temperatures above 30°C increase I. scapularis mortality, reduce oviposition, and inhibit host-seeking activity.[32,55,56] However, maximum summer temperatures rarely exceed this threshold in Minnesota, and ticks are capable of surviving warmer temperatures in areas of higher humidity. Warm temperatures experienced by ticks in Minnesota may accelerate the tick’s life cycle, thus supporting higher densities of ticks.[57]

In Minnesota, increased incidence of I. scapularis–borne human disease cases may be a result of changes in tick abundance and/or infection prevalence, range expansion of I. scapularis, changes in reporting and surveillance, or perhaps increased encounter rates between humans and ticks.[5,6,18,20,58,59] Although our model generally supports conclusions from previous models and epidemiologic surveillance activities, the areas of discordance that we highlighted are worthy of further investigation to determine if these are areas of model misclassification, previously underrecognized focal areas of risk, or areas of emergence. For this study, sites were located exclusively on public land; these sites may represent opportunities for recreational risk of acquiring I. scapularis–borne diseases and should motivate the use of preventative measures to avoid or mitigate the impacts of tick bites while recreating in tick habitat. Sampling on residential properties was not undertaken and, therefore, extrapolation of the models to residential areas should be viewed with caution and underscores the need to evaluate model predictions on residential properties. Understanding when and where ticks are active is fundamental to preventing tick-borne diseases. However, human activities that reduce risk, including avoiding tick habitat when ticks are active, using repellents containing 20–30% DEET, performing daily tick checks and removing ticks promptly, and showering after being outdoors, remain imperative to preventing tick-borne diseases.[60]

Supplemental information