Increasing habitat suitability in the United States for the tick that transmits Lyme disease: a remote sensing approach.

The warnings about the spread of (italic)Ixodes scapularis(/italic), one of the vectors of Lyme disease, into the United States are based on reports about regional distribution and increasing local abundance. In a modeling approach, I used the recorded, current distribution of this tick and remotely sensed bioclimatic factors over the United States to establish the changes of habitat for this tick since 1982 and to detect the areas with factors adequate to support tick colonization. Results indicate the geographic expansion of areas with adequate habitat suitability in the period 1982-2000. A discriminant analysis of counties with different degrees of habitat suitability shows that the increase in winter temperatures and in vegetation vitality (as a direct consequence of higher rainfall) is key to habitat switch from unsuitable to suitable.

Arthropods are extremely sensitive to climate. Throughout the 20th century, public health researchers have understood that climate circumscribes the distribution of mosquito-borne diseases, whereas weather affects the timing and intensity of outbreaks (1). A growing number of investigators propose that vector-borne diseases shift their range in response to climate change (2). Lyme disease has become the most prevalent vector-borne illnesses in the United States, and its distribution is expanding (3). One of the tick species involved in Lyme disease transmission is Ixodes scapularis (Say), which is currently well established in several areas of the northeastern, southeastern, and midwestern United States. Knowledge of the environmental factors that influence the survival of this vector has acquired special importance, because the regional distribution and local abundance of I. scapularis in North America have increased over the past two decades, apparently at an unprecedented rate (4).
Considerable effort has been devoted to studying the ecologic and social factors that influence where and in what abundance this tick is found. However, the integration of these findings into a paradigm that explains tick colonization of new regions and changes in local abundance is still being developed. Added to this is concern that global climatic changes may create widespread conditions that permit expansion into areas that are not currently infested. We now have evidence that I. ricinus, an important tick species of the western Palearctic, has increased its density and geographic range (latitude and altitude shifts) in Sweden because of changes in the medium-term climate (10-15-year period) in that region (5). The aim of the present study was to detect changes in the habitat for I. scapularis in the United States in the period 1982-2000. I used a modeling approach together with remotely sensed climate and vegetation features at a medium resolution of 8 km to elucidate the temporal drift of these abiotic variables and how they influence the forecasted habitat suitability for the tick.

Materials and Methods
I took the tick distribution data for this study from the distribution of I. scapularis published by U.S. counties (6). This compilation covers the known range of the tick from published records and collections as of 1998. In this compilation, I considered I. scapularis ticks to be "reported" in a county if at least one tick at any life stage had been identified. I also considered tick populations "established" in a county if at least six ticks or two of the three tick life stages had been identified in a single collection period. Both definitions have a biologic meaning. More than one stage represents a reproductive population, whereas the collection of just a single tick implies the initial stages in the event sequence from invasion to colonization. In both cases, the county is the geographic unit, because it provides the ecologic background of collections (collecting point and its surroundings). In the present study, I used this published distribution as the baseline for a procedure based on a rangestandardized, point-to-point similarity metric (the DOMAIN procedure) that provides a robust method for modeling potential distribution according to habitat suitability (HS), as defined by remotely sensed temperature and vegetation data (7).
Bioclimatic factors have been observed through the period 1982-2000 by remote sensing. I obtained the series of the advanced very high-resolution radiometer sensor of the National Oceanic and Atmospheric Administration (NOAA) satellite, already georeferenced and processed to eliminate artifacts, for a region over the Nearctic covering Canada, the United States, and most of Central America (8). I obtained data as 10day intervals (decadals) with a resolution of 8 km/pixel. I obtained imagery for both channels 4 and 5, as well as for the already derived Normalized Derived Vegetation Index (NDVI). The NDVI is a measurement of vegetation stress and is directly related to relative humidity and rainfall, as well as to vegetation growth. I calculated the surface temperature (T) with data from channels 4 and 5 using a window split algorithm, according to instructions of the data set. With decadal data about temperature and NDVI, I calculated each of the four yearly variables: mean, absolute maximum, monthly mean, and monthly mean maximum for both T and NDVI. Temperatures are important in defining the development rates of the tick population and the hostfinding rates, because tick activity is temperature mediated. Incorporating NDVI into computations provides an indirect evaluation of relative humidity in the habitat, as one of the most important variables (alone or together temperature) in the modeling of tick survival. I computed these eight variables both for each year in the period 1982-2000 and for the entire period of study.
The modeling procedure uses a point-topoint similarity metric to assign a classification value to a candidate site based on the proximity in environmental space of the most similar record site. I defined the distance d between two points A and B in the Euclidean p-dimensional space as This metric uses range standardization to equalize the contribution from each bioclimatic attribute. This method of standardization is preferred over variance standardization in this application because it is less susceptible to bias arising from dense clusters of sample points. The complementary similarity measure R AB is R AB = 1 -d AB . R is constrained between 0 and 1 for points within the ranges of the equation but may yield negative values for points outside this range. I defined the maximum similarity (S A ) between candidate point A and the set of known record sites T m as Evaluating S for all grid points in a target area generates a matrix of similarity values.
I performed two kinds of HS calculations. The first calculation used the mean abiotic variables for the entire 1982-2000 period to check the robustness of the model by comparing the computed HS for this period and the known distribution of the tick. The second calculation used the abiotic factors for defined multiyear time steps to identify how the HS changed with time. A 1-year time period is too short to reveal the trend of abiotic variables and thence the effect on habitat suitability. Therefore, I decided to use yearly data to detect the cycles in abiotic variables, and then I recomputed HS with the abiotic variables of the multiyear period.
It is of interest to know which bioclimatic variables are the most involved in habitat changes, as well as the average trends of these variables in years predicted to be unsuitable versus those modeled as suitable. I performed a discriminant analysis using bioclimatic variables from counties of five states that displayed different habitat values at given time periods. Because of results   obtained in the preliminary modeling, I decided that a county switched into positive suitability when HS values were higher than 0.5. To increase the performance of the discriminant analysis, and to convert raw abiotic variables into values with biologic meaning, I transformed the variables from the original 10-day recording intervals to interval cumulative relatives. The groups of counties selected for this analysis were in Pennsylvania and New York (collectively called the northeastern group), Wisconsin, Virginia, and North Carolina. I selected these counties because they represent a relatively wide range of abiotic conditions, and because they maintain positive suitability after switching from an unsuitable habitat. Table 1 lists data for every selected county, with notes about date of habitat switching from unsuitable to suitable together with the modeled value of HS before and after habitat switch. I undertook the discriminant analysis with two objectives: to determine the variables that are most important in the switch from nonsuitability to suitability, and to understand whether both single counties and groups of counties can be easily separated according to simple abiotic variables recorded before and after the habitat became suitable. Variables for the discriminant analysis consisted of temperatures and NDVI data for each county, ranked according to the number of decadals with a variable within certain intervals. For temperature, the intervals (in°C ) were < -12, -12 to -7, -7 to 0, 0 to 5, 5 to 8, 8 to 10, and > 10. For NDVI, the intervals were < 0, 0 to 0.1, 0.1 to 0.25, 0.25 to 0.35, and > 0.5. I selected temperature ranges according to reports on the critical factors for the activity and survival of the ticks, obtained from summarized data (9). I selected NDVI intervals to define vegetation categories (10). I used these values as predictors for the four groups of counties, before and after habitat became suitable, in a discriminant analysis based on the Mahalanobis distance.

Results
Counties recorded as "established" were suitable in the model results, with HS values ranging from 0.8 to 1. These counties were concentrated (Figure 1   However, the modeling also detects suitable habitat in the same range (0.5-0.8) in 1,066 additional counties for the period 1982-2000 where I. scapularis has never been noticed (specificity, 0.78). These counties cover an extensive area from eastern Kansas and Oklahoma through parts of Missouri, Arkansas, Kentucky, Tennessee, Ohio, West Virginia, Virginia, both Carolinas and Georgia, together with isolated patches in northeastern counties, Wisconsin, Minnesota, and Michigan.
When I performed HS modeling using bioclimatic variables along time steps in the period 1982-2000, I saw a trend toward increasing HS. Obviously, faint short-term changes of the abiotic conditions in a short time period (i.e., a year) did not account for the presence or absence of the tick in the same period. Only when the model considers conditions from 3-6 years can the trend in HS be deduced with confidence. The search for these medium-term cycles in abiotic variables reveals four well-defined "waves" that correlate with drastic changes in HS. In the period 1982-1985, the model detected a total of 604 counties with HS ranging from 0.5 to 0.8 (Table 2, Figure 2). The distribution of these counties almost coincides with that of counties where the tick has been reported. In the period 1986-1989 (Figure 3) Two variables-the number of decadals with mean temperature within the interval 0 and 5°C, and the number of decadals with NDVI over 0.5-explained more than 71% of total variability and provided the better separation of the clusters of counties before   and after habitat suitability switching ( Figure  6). Ellipses of 95% confidence show that the model achieves good separation within and between the groups of counties with these variables for both the Wisconsin and northeastern counties. The separation between Virginia and North Carolina is optimal. Table 3 includes data about the eigenvalues in the discriminant analysis, showing the strength of associations between variables and separation of counties. Figure 7 shows that the values of these two variables from the averaged data for periods of positive HS tends to be greater than the values for periods of negative HS. Comparing the averaged conditions, counties with positive suitability may have up to two decadals more of temperatures within the range 0-5°C, and up to four decadals more of NDVI values over 0.5.

Discussion
Because only a limited sample of the actual distribution of a parasite is usually achieved in practice, epidemiologists rely increasingly on spatial models to indicate the potential distribution range and the likely consequences of environmental impacts on a given species. The findings presented in this study support the hypothesis that the habitat availability for I. scapularis in the United States is undergoing geographic changes, allowing the development and survival of this tick in zones where 20 years ago the habitat was unsuitable. However, the number of counties with highly adequate habitat remains constant throughout the period of study. In whole, these data correspond well with reports of an increase in the regional distribution and local abundance of this tick in North America. Reports concerning regional changes in the distribution of I. scapularis in parts of Maine (11), Connecticut (12), New York (13), and Wisconsin (14) fit well with the changes in HS calculated herein. We have no available evidence of such an expansion in the southern United States (5), and the data calculated in the present study reveal no change in habitat availability for I. scapularis over this region. Bearing in mind the sensitivity and specificity of the modeling procedure, as evidenced by the accuracy with which it predicts counties with established populations, these findings seem to represent potential areas of suitable habitat for first events in the colonization of new zones, rather than false positives.
The geographic extension of suitable habitat may not be a direct reflection of the existence of I. scapularis inside a particular range, although HS measurement provides an adequate estimation of the areas susceptible to tick colonization. The expansion of I. scapularis over time may not occur at the same rate as do changes in HS, because expansion depends on the invasion of newly suitable areas with specimens from already established and active foci. This process could easily come about through the movement of migratory birds or mammals from well-established areas to contiguous regions, and seems to be the rule in the northeastern and midwestern United States (11,12). The ability of passerine birds to introduce ticks into a nonendemic area has recently been confirmed (13). Invasive movements of this tick are slow (14), and areas such as Virginia, North Carolina, and Missouri are far enough from active foci to prevent introduction by ranging mammals. In addition, they are distant from the main north-south bird migratory routes, which seem to be more common than east-west ones (15). This paucity in the introduction of specimens to newly available areas would explain the small number of I. scapularis collections in the range from the Allegheny Mountains to the Mississippi Valley. Analysis of abiotic factors suggests that this area became suitable in the period 1990-1993, and then remained almost unchanged through 1994-2000. Evaluating the sampling efforts in the many different reports that comprise the data source used for this study is not possible. Therefore, making inferences about the lack of identifications in individual counties is difficult. Finally, year-to-year variations must be taken into account. Habitat is evaluated with remotely sensed bioclimatic variables detected along homogeneous periods, according to variance minimization down a temporal window of variable length. Small year-to-year changes in HS can not taken into account by this method. It remains possible that subtle climate variations from year to year within the same period render the habitat unsuitable at intermittent periods of sufficient duration to prevent tick colonization. However, because HS computations in this study for the period 1982-1985 clearly reflect the situation as currently known, it is wise to consider that HS calculations as performed for 1986-2000 should be regarded as an early warning of future colonization processes, if trends of abiotic variables remain.
Abiotic variables changed in the groups of counties involved into discriminant analysis. Most highly involved variables are the Articles • Climate and trend in habitat suitability for I. scapularis Environmental Health Perspectives • VOLUME 110 | NUMBER 7 | July 2002 639 Figure 6. Discriminant analysis separation of the groups of counties according to the variables explained in "Materials and Methods." The first axis is loaded with the number of decadals with an average temperature within the range 0-5°C. The second axis is loaded with the number of decadals with NDVI over 0.5. The analysis discriminates the counties included into two groups: before (-) and after (+) the habitat switched to suitability. Table 1 lists counties included.    Table 1 lists counties included. Data are the means of the periods when the habitat was unsuitable (predicted HS < 0.5) and the means for the years when the habitat was adequate (predicted HS > 0.5).