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

The warnings about the spread of Ixodes scapularis 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. Key words: climate, habitat suitability, Ixodes scapularis, Lyme disease, remote sensing.

**********

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 range-standardized, 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 10-day 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 host-finding 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-to-point 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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

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.sub.AB] is [R.sub.AB] = 1 - [d.sub.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.sub.A]) between candidate point A and the set of known record sites [T.sub.m] as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

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 bio-climatic 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 north-eastern 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 [degrees]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) along the northeastern Atlantic Coast from northern Maine through Maryland, along the southeastern Atlantic Coast and the eastern Gulf coastal areas, and in the north-central states of Michigan (upper peninsula), Wisconsin, Minnesota, and adjacent Iowa. I found a scattered distribution of positive counties in the south-central states of Missouri, Mississippi, Arkansas, Oklahoma, Louisiana, and eastern Texas. The model captured every county where the tick has reproductive populations and produced no false-positive counties in the HS output. These results suggest that the abiotic variables included in the modeling of potential distribution are important factors in the establishment of permanent populations of the tick. The model detected all of the counties where I. scapularis has been reported (Figure 1), indicating no false negatives in this category. The 593 counties in which the tick has been reported are generally clustered around the concentrated areas of established populations and have HS values ranging from 0.5 to 0.8. These results show a sensitivity of the model of 1 (no false negatives).

[FIGURE 1 OMITTED]

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), most of the suitable counties remained under this category, but new areas of positive suitability are detected in Kentucky and West Virginia. This nucleus is greatly expanded in the period 1990-1993 (Figure 4), when suitable habitat extends over areas of Virginia and North Carolina. Finally, in the period 1994-2000 (Figure 5), HS extends over areas of Kansas and Missouri, as well as in the northern range of distribution (Minnesota, Wisconsin, and northeastern counties on the United States-Canada border).

[FIGURES 2-5 OMITTED]

Two variables--the number of decadals with mean temperature within the interval 0 and 5[degrees]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 6 OMITTED]

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[degrees]C, and up to four decadals more of NDVI values over 0.5.

[FIGURE 7 OMITTED]

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 dearly 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 number of decadals with temperatures within 0-5[degrees]C (up to two decadals more in counties with predicted positive suitability) and the number of decadals with NDVI higher than 0.5. This temperature range is commonly found in winter; thus, the increase of this feature should be considered an expression of milder winters. These temperature conditions are the lower threshold limit for tick activity and metamorphosis (9). However, they may not have a direct influence on the tick life cycle, but on major changes in the functional relationships of the ecosystems, such as an increase in host animals or habitat modifications that may strongly influence tick survival and abundance. Data in Figure 7 show that both variables of concern changed together when HS of selected counties changed from unsuitable to suitable. The increase in the number of decadals with T within 0-5[degrees] must be complemented with an adequate increase in the number of decadals with NDVI higher than 0.5. The rise of temperature alone is not the only cause in habitat suitability changes.

Although climate warming and milder winters are well-documented features in the analysis of climate trends in recent years, the increase of NDVI values is a recently discovered trait (16). Most reports acknowledge the high humidity requirements for the survival of I. scapularis, which are reflected in the high NDVI values in suitable counties. Furthermore, it has been proposed (17) that harsh winter temperatures restrict this tick to the milder conditions of the Atlantic seaboard, which agrees with the data presented here. The figures provided in the present work show that the increase in geographic extension of I. scapularis is a feature driven in part by climate. Other features, such as changes in land use, habitat fragmentation, or variations in host distributions, have not been considered here because of the lack of an adequate, homogeneous database. It seems possible that some of these biotic variables may explain further the trend in habitat suitability for I. scapularis. The biologic meaning of these figures is that, leaving aside other bioclimatic features that define an area as suitable for tick development and survival, counties with adequate habitat have two main features: milder winters and higher NDVI (higher rainfall) in the vegetation growing season. Considering that milder winters and more humid conditions are dominant features in bioclimatic variables, I. scapularis populations have a wider area to colonize where hosts are already available, and pose a serious challenge to human health.

Table 1. Counties selected for the discriminant
analysis based on the habitat switch from unsuitable
(HS < 0.5) to suitable (HS > 0.5).

County (state)       Switch  HS (B)  HS (A)

Saint Lawrence (NY)   1986    0.01    0.58
Hamilton (NY)         1990    0.04    0.59
Genesee (NY)          1986    0.07    0.65
Chautauqua (NY)       1986    0.01    0.63
Allegany (NY)         1986    0.02    0.63
Cortland (NY)         1986    0.01    0.64
Bradford (PA)         1986    0.01    0.68
Tioga (PA)            1986    0.03    0.59
Susquehanna (PA)      1986    0.01    0.58
Venango (PA)          1986    0.01    0.61
Fayette (PA)          1986    0.03    0.61
Green (WI)            1994    0.02    0.64
Dunn (WI)             1994    0.03    0.61
Shawano (WI)          1994    0.01    0.59
Menominee (WI)        1990    0.03    0.62
Ashland (WI)          1986    0.02    0.64
Juneau (WI)           1986    0.04    0.61
Door (WI)             1990    0.01    0.69
Polk (WI)             1994    0.03    0.65
Saint Croix (WI)      1994    0.01    0.64
Waupaca (WI)          1994    0.01    0.59
Langlade (WI)         1994    0.01    0.69
Walworth (WI)         1990    0.03    0.58
Pepin (WI)            1990    0.01    0.61
Iron (WI)             1990    0.02    0.61
Prince William (VA)   1994    0.03    0.59
Arlington (VA)        1994    0.01    0.61
Rappahannock (VA)     1990    0.05    0.61
Fairfax (VA)          1986    0.04    0.62
Harrisonburg (VA)     1990    0.11    0.57
Orange (VA)           1990    0.06    0.53
Fredericksburg (VA)   1990    0.04    0.58
Westmoreland (VA)     1994    0.02    0.54
Northumberland (VA)   1994    0.06    0.56
Amherst (VA)          1990    0.05    0.52
Botetourt (VA)        1990    0.08    0.59
Powhatan (VA)         1990    0.04    0.57
Montgomery (VA)       1990    0.03    0.54
Dickenson (VA)        1990    0.02    0.57
New Kent (VA)         1990    0.08    0.52
Bedford (VA)          1990    0.09    0.61
Alleghany (NC)        1990    0.01    0.58
Surry (NC)            1990    0.02    0.62
Stokes (NC)           1990    0.01    0.61
Mecklenburg (NC)      1990    0.04    0.57
Graham (NC)           1986    0.04    0.59
Gaston (NC)           1994    0.05    0.61
Hyde (NC)             1994    0.08    0.62
Polk (NC)             1986    0.04    0.59
Macon (NC)            1986    0.01    0.59
Durham (NC)           1990    0.03    0.52
Madison (NC)          1990    0.05    0.59
Alexander (NC)        1990    0.06    0.61
Pearson (NC)          1990    0.05    0.64
Ponder (NC)           1994    0.05    0.68
Cleveland (NC)        1990    0.01    0.64

Data include the name of the county, the year in which
habitat switched as from model data (switch) and the HS
obtained for that county before (B) and after (A) the suitability
swap. After suitability change, the HS of these
countries remained within the limits of positive suitability.
Table 2. Time evolution of the number of counties with HS within the
range 0.5-0.8 (medium HS), as detected with the range-standardized,
point-to-point similarity metric and the remotely sensed bioclimatic
variables, and the number of counties that remained within the
category of suitable habitat in subsequent time periods.

Time period  No. counties  1982-1985  1986-1989  1990-1993

1982-1985         604
1986-1989         952         542
1990-1993       1,301         569        881
1994-2000       1,552         586        927       1,213
Table 3. Eigenvalues of the discriminant analysis
performed with the counties included in Table 1,
and the set of raw variables converted to accumulated
values.

Variable                        Percent variance  Cumulative %

0[degrees] < T < 5[degrees]          51.76           53.76
NDVI > 0.5                           20.03           71.79
-12[degrees] < T < -7[degrees]       14.12           85.91
0.25 < NDVI < 0.35                    2.05           87.96
5[degrees] < T < 8[degrees]           0.85           88.81
8[degrees] < T < 10[degrees]          0.12           88.93
NDVl < 0                              0.08           89.01
T > 10[degrees]                       0.02           89.03
-7[degrees] < T< 0[degrees]         < 0.01           89.03
0 < NDVI < 0.1                      < 0.01           89.03
0.1 < NDVl < 0.25                   < 0.01           89.03
T < -12[degrees]                    < 0.01           89.03

Variables refer to the number of decadals in a year
within the given interval of temperature (T in degrees
Celsius) and NDVI.

REFERENCES AND NOTES

(1.) Dobson A, Carper R. Biodiversity. Lancet 342:1096-1099 (1993).

(2.) McMichael A, Haines A, Sloff R. Kovats S, eds. Climate change and human health. Geneva:World Health Organization, 1996.

(3.) Centers for Disease Control and Prevention. Lyme disease-United States. Morb Modal Wkly Rep 45:481-484 (1996).

(4.) Wilson ML. Distribution and abundance of Ixodes scapularis (Acari: Ixodidae) in North America: ecological processes and spatial analysis. J Med Entomol 35:446-457 (1998).

(5.) Lindgren E, Talleklint L, Polfeldt T. Impact of climate change on the northern latitude limit and population density of the disease-transmitting European tick Ixodes ricinus. Environ Health Perspect 108:119-123 (2000).

(6.) Dennis DT, Nekamoto TS, Victor JC, Paul WS, Piesman J. Reported distribution of Ixodes scapularis and Ixodes pacificus (Acari: Ixodidae) in the United States. J Med Entomol 35:629-638 (1998).

(7.) Carpenter G, Gillison AN, Winter J. DOMAIN: a flexible modelling procedure for mapping potential distribution of plants and animals. Biodiv and Cons 2:667-680 (1993).

(8.) GSFC Earth Sciences. Atmospheric Dynamics Data and Resources. Available: http://daac.gsfc.nasa.gov/[cited 11 April 2002].

(9.) Mount GA, Haile DG, Daniels E. Simulation of black-legged tick (Acari: Ixodidae) population dynamics and transmission of Borrelia burgdorferi. J Med Entomol 34:461-484 (1997).

(10.) Hugh-Jones M. Applications of remote sensing to the identification of habitats of parasites and disease vectors. Parasitol Today 5:244-251 (1989).

(11.) French JB, Schell WL, Kazmierczak JJ, Davis JP. Changes in population density and distribution of Ixodes dammini(Acari: Ixodidae)in Wisconsin during the 1980s. J Med Entomol 29:723-728 (1992).

(12.) Fish D. Environmental risk and prevention of Lyme disease. Am J Med 28(suppl 4A):2S-8S (1995).

(13.) Smith RP Jr, Rand PW, Lacombe EH, Morris SR, Holmes DW, Caporale DA. Role of bird migration in the long-distance dispersal of Ixodes dammini, vector of Lyme disease. J infect Dis 174:221-224 (1996).

(14.) Kitron U, Kazmierczak JJ. Spatial analysis of the distribution of Lyme disease in Wisconsin. Am J Epidemiol, 145:558-566 (1997).

(15.) National Atlas of the United States. Available: http://www.nationalatlas.gov [cited 10 August 2001].

(16.) Climate change impacts on the United States: The potential consequences of climate variability and change. Available: http://www.usgcrp.gov and http://sedac.ciesinorg/NationalAssesment [cited 10 August 2001].

(17.) McEnroe WD. The restriction of the species range of Ixodes scapularis Say, in Massachusetts by fall and winter temperature. Acarologia 18:618-625 (1977).

Address correspondence to A. Estrada-Pena, Department of Parasitology, Faculty of Veterinary Medicine, Miguel Servet, 177, 50013-Zaragoza, Spain. Telephone: [phone]. Fax: [phone]. E-mail: [e-mail]

The data used here include data produced through funding from the Earth Observing System Pathfinder program of NASA's Mission to Planet Earth in cooperation with the National Oceanic and Atmospheric Administration. The data were provided by the Earth Observing System Data and Information System (EOS-DIS), Distributed Active Archive Center at Goddard Space Flight Center, which archives, manages, and distributes this data set.

Received 4 September 2001; accepted 10 December 2001.

Agustin Estrada-Pena
Department of Parasitology, Faculty of Veterinary Medicine,
Miguel Servet, Zaragoza, Spain