Ecological Applications, 17(1), 2007, pp. 129–139 Ó 2007 by the Ecological Society of America

PERSISTENTLY HIGHEST RISK AREAS FOR HANTAVIRUS PULMONARY SYNDROME: POTENTIAL SITES FOR REFUGIA GREGORY E. GLASS,1,4 TIMOTHY SHIELDS,1 BIN CAI,1 TERRY L. YATES,2

AND

ROBERT PARMENTER2,3

1 The W. Harry Feinstone Department of Molecular Microbiology and Immunology, The Johns Hopkins Bloomberg School of Public Health, 615 North Wolfe Street, Baltimore, Maryland 21205 USA 2 The Department of Biology and the Museum of Southwestern Biology, The University of New Mexico, Albuquerque, New Mexico 87131 USA 3 The Valles Caldera National Trust, 18161 State Highway 4, P.O. Box 359, Jemez Springs, New Mexico 87025 USA

Abstract. Interannual variation in the number of cases of human disease caused by hantaviruses in North America has been hypothesized to reflect environmental changes that influence rodent reservoir populations. This hypothesis postulates that when cases are rare reservoir populations are geographically restricted in patches of suitable habitat. Identifying these sites, which is needed to test the hypothesis, has proven to be a challenge. Satellite imagery of the U.S. Southwest has shown associations among the likelihood of human hantaviral disease and increases in the rodent populations, as well as increased prevalence of Sin Nombre virus (SNV) in rodent populations. In this study we characterize local areas that had environmental signatures that persisted as predicted highest risk sites for human disease through much of the 1990s. These areas represent a small percentage (0.3%) of the region. Exploratory analyses indicate that these areas were not randomly distributed, but were associated with certain landscape characteristics. Characteristics of elevation, slope, aspect, and land cover were associated with persistent high risk. Using multivariate Poisson regression to control for confounding effects, sites with deciduous- or mixed-forest land cover on moderate to steep slopes (.58) above 2130 m elevation were associated with increasing numbers of years at highest risk. These are candidate locations for refugia. Sites associated with cleared ground or shrubland were less often associated with high risk compared to reference conditions. The seasonal patterns of vegetation growth in persistently high-risk areas were compared to matched locations using MODIS (moderate resolution imaging spectroradiometer) NDVI (normalized difference vegetation index) during a time of a severe drought in the region from 2002 to 2004. Despite the drought and regardless of land cover, the NDVI in persistently highest risk areas had an early onset, with significantly higher levels of green vegetation that lasted longer than at comparable sites. These observations identify locations that can be monitored for the abundance of P. maniculatus and presence of SNV. If these sites are refugia, we predict they will be occupied by infected deer mice when other monitored sites are unoccupied. Key words: deer mouse; GIS; hantavirus; hantavirus pulmonary syndrome; Peromyscus maniculatus; remote sensing; Sin Nombre virus; trophic cascade hypothesis.

INTRODUCTION Hantavirus Pulmonary Syndrome (HPS) is a human disease caused by infection with various hantaviruses. These viruses are zoonotic agents that naturally infect and are maintained in murine rodents (Yates et al. 2002). Extensive surveys, coupled with phylogenetic analyses of rodents and viruses, implicate single species of rodents as hosts of single hantaviruses with occasional spillover to other species, including humans (Childs et al. 1994. Mills et al. 1997, Yates et al. 2002). HPS was recognized following an outbreak of pulmonary disease among humans in the U.S. Southwest in 1993, and infection with Sin Nombre virus (SNV) was traced to its Manuscript received 25 April 2005; revised 7 March 2006; accepted 3 May 2006. Corresponding Editor: N. T. Hobbs. 4 E-mail:[email protected]

reservoir, Peromyscus maniculatus, the deer mouse (Nichol et al. 1993, Childs et al. 1994). The number of cases of HPS varies both geographically and among years (e.g. Hjelle and Glass 2000). The causes of these fluctuations are not clear but are presumed to be due, generally, to changes in contacts between infectious mice and humans (Childs et al. 1995, Hjelle and Glass 2000). Parmenter and colleagues (1993) postulated that one of the main drivers for the interannual variation in HPS cases was changes in environmental conditions that were thought to be due to fluctuations in the patterns of precipitation and temperature. These fluctuations affected the abundance of deer mice and subsequently the human risk of HPS (Parmenter et al. 1993, Yates et al. 2002). They hypothesized that as precipitation increased and temperatures were ameliorated net primary productivity (NPP) increased with attendant increases in ground-

129

130

dwelling arthropod populations. Increased NPP acted both directly and indirectly on local deer mouse populations, resulting in larger populations. Increased rodent densities were associated with increased viral transmission among rodent population members and larger numbers of dispersing animals that, in turn, led to increased transmission to humans who came in contact with these animals (Yates et al. 2002). This ultimately became known as the Trophic Cascade Hypothesis (TCH). The TCH is being evaluated in several studies, though much work remains (Mills et al. 1997, Yates et al. 2002). According to the TCH, as environmental conditions become more severe, many local rodent populations decline and become extinct. However, in some locations the environment allows survival of sufficiently large rodent populations that the virus persists. It is from these sites that SNV later reappears and spreads. A key step in evaluating the TCH is to identify the conditions needed to maintain large reservoir populations that carry SNV and where they occur. Large reservoir populations are needed because hantaviruses are horizontally transmitted agents (agents that are spread among members of the population after birth). In the absence of vertical transmission (from mother to offspring, typically in utero), sufficient numbers of susceptible individuals must be recruited and acquire virus from infectious individuals so that the virus does not go extinct locally. Parmenter et al. (1993) proposed that locations where these environmental conditions existed served as ‘‘refugia’’ for the virus. Studies of rodent reservoir populations using traditional approaches have had difficulty testing the TCH because of challenges in locating refugia. Many of the local rodent populations monitored as part of the longitudinal studies have been subject to repeated extinction and only sporadic SNV infection (e.g., Calisher et al. 2001, Yates et al. 2002 and references therein), suggesting that these sites do not fit the conditions for refugia. A primary challenge with testing the TCH is that the environmental conditions needed to maintain P. maniculatus and SNV are not defined by the theory, but only by their existence. In earlier work (Glass et al. 2000), the risk of HPS was associated with environmental conditions characterized by satellite imagery nearly 10 months prior to disease outbreak. The ‘‘best’’ linear classifier, using Landsat Thematic Mapper (available online)5 (TM) imagery, combined TM bands 1, 5, and 7 and elevation (.2094 m). Subsequent studies of local rodent populations during an El Nin˜o/Southern Oscillation (ENSO) event in 1997–1998 examined areas where the environmental signature indicated, for humans, high risk, low risk, or areas that fluctuated between low and high risk during a two-year period. The prevalence of SNV infection 5

Ecological Applications Vol. 17, No. 1

GREGORY E. GLASS ET AL.

hhttp://eros.usgs.gov/products/satellite/tm.htmli

started from a common baseline in all cases, and the characteristics of local rodent populations in sites where human risk fluctuated between years were similar to low-risk sites. Sites that were predicted to be high risk for two consecutive years were characterized by portions of the mouse population known to be at highest risk for infection. The prevalence of SNV infection in deer mice from these sites exceeded 30% by the second year of study, even though they had started with a prevalence of infection comparable to low-risk sites (Glass et al. 2002). The apparent need for sites to persist as high risk for at least two years before infection prevalence increased was consistent with the interpretation that virus had been introduced into these local populations as suitable conditions persisted during the study. Here, we characterize physical conditions at sites where high risk of HPS persists for extended periods of time (persistently highest risk areas). The goals were to identify locations that could meet the criteria for refugia and to characterize their physical features with the intent of subsequently monitoring local rodent populations at these sites. We use the onset of a severe, long-term, regional drought to examine the seasonal dynamics of vegetative patterns, as measured by normalized difference vegetation index (NDVI), in persistently highest risk areas, compared to the remainder of the region. The rationale for characterizing the environment based on human risk is based on studies by Childs and colleagues (1995) who showed that the only demonstrable difference between households where HPS occurred and those where it didn’t were the number of infected P. maniculatus captured. Thus, sites that are at highest risk for HPS over many years also should be sites with environmental signatures consistently associated with large numbers of infected P. maniculatus. If the TCH is true, then environmental patterns in these areas could correspond to conditions that favor the survival of deer mice populations. METHODS Persistently highest risk (PHR) sites were identified in a region of the U.S. Southwest by extending previously described methods (Glass et al. 2000). Human risk was modeled from the case-control data set with logisticregression analyses, in which linear combinations of Landsat TM 5 bands and elevation .2094 meters were used. Annual models of risk were generated for five years, symmetrically spanning most of the 1990s (1992, 1993, 1995, 1997, and 1998). The annual HPS risk maps were coregistered and overlaid. The highest risk areas within a single year were located by examining the geographic extent for predicted values from the logistic-regression analysis. Thresholds for the predicted value corresponding to ;2.5% of the study area were selected. This choice was arbitrary but was used because it provided good concordance between the receiver operator characteristic (ROC) of NDVI and the ROC of predicted human risk in 1992

January 2007

HANTAVIRUS REFUGIA

(Glass et al. 2000). We also evaluated alternative thresholds of 5%, 10%, and 20%, which represented more extensive areas, but these did not alter the qualitative results. Exploratory analyses Each year’s highest risk areas were dichotomized above and below the 2.5% threshold, and persistently highest risk areas were identified by overlaying the annual risk maps. Sites that were in 2.5% highest risk categories throughout the study period were categorized as PHR. The spatial stability of PHR sites was examined by recording the number of pixels that were in the highest risk areas 0, 1, 2, 3, 4, or 5 years, and these were compared with the observed numbers expected under a binomial distribution with P ¼ 0.025. If highest risk sites persisted, then an excess of sites in the highest risk categories with four or five years recurrence and a paucity of sites with no or one year in the highest risk categories was expected. Failure to reject the hypothesis would indicate that pixels in highest risk categories varied geographically among years. As part of the exploratory analyses, the spatial relationships of PHR areas with physical and landcover characteristics of the region were examined to determine if potential ‘‘refugia’’ corresponded to recognizable landscape features. PHR sites were overlaid with digital raster graphics (DRGs) generated from U.S. Geological Survey 1:24 000 scale topographic maps. Images were visually examined to identify aspects of the landscape associated with PHR areas. Digital elevation models (DEMs), obtained from the U.S. Geological Survey (USGS) National Elevation Database in 7.5-min data, were used to estimate the elevation, slope, and aspect of sites. The DEM files were coregistered to the PHR areas. Elevation and aspect of PHR areas were compared with the remaining study area by Komolgorov-Smirnov two-sample test and Watson’s U2 test, respectively. The null hypotheses were that PHR areas were random samples of the underlying patterns of elevation and aspect within the study region. Land-cover characteristics were derived from U.S. EPA National Land Cover Data (NLCD) from the Multi-Resolution Land Characteristics (MRLC) Consortium. NLCD level II 1992 classifications were extracted for all pixels in PHR areas, as well as the entire study region. The proportions of the study area and PHR areas were compared in a 2 3 C table (two rows by the number of columns necessary for the comparison). Land-cover categories that individually represented ,0.5% of the study region were aggregated into an ‘‘other’’ class for the exploratory analysis. Multivariate analyses It was anticipated that many environmental variables, such as elevation and land cover, would be correlated among themselves, making it difficult to control for confounding effects. Multivariate Poisson regression

131

was performed to reduce confounding. Outcomes were the numbers of years (0–5) that an individual pixel was scored as highest risk. The predictor variables were derived from the environmental variables used in the exploratory analyses. Indicator variables were generated based on these variables and effects defined relative to reference classes. The large number of pixels in the study area (;98.5 million pixels), made it infeasible to perform a single analysis. Instead, a Monte Carlo analysis was conducted to obtain point and interval estimates of the regression coefficients. Approximately 11 000 pixels were randomly selected without replacement in the study area and their environmental variables extracted. Poisson regression was performed and the coefficient for each variable was generated while including the effects of all other variables. The process was repeated 200 times, with replacement between samples. The median value of each coefficient derived from the 200 repetitions was chosen as the point estimate for the environmental parameters. The empirical interval estimate excluded the extreme 2.5% of values for each parameter. Parameters with estimates that overlapped 0 were not considered to differ significantly from the reference category but were retained in the model to control for their impact on the remaining predictor variables. Environmental dynamics in PHR areas Seasonal dynamics in vegetation greenness of PHR areas were characterized by the normalized difference vegetation index (NDVI). NDVI values were generated as a data product (MOD13) from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor for 2002, 2003, and the first quarter 2004. Spatial resolution of the MODIS NDVI was 250 m. Temporal resolution was 16 days. NDVI, a commonly used indicator of vegetation growth/productivity (Lillesand et al. 2004), was previously compared with the HPS risk algorithm (Glass et al. 2000). It compared favorably with the risk algorithm in highest and lowest risk areas, though it was not as sensitive in the intermediate range of risk. However, intermediate risk areas were not included in the PHR sites. The PHR image that was overlaid on the NLCD file was used to stratify PHR areas and the remaining portion of the study site by land-cover type. We selected a series of locations from PHR sites and randomly selected lower risk sites stratified into the four most common land-cover categories: evergreen forest, shrubland, grasslands/herbaceous, and pasture/hay accounted for more than 95% of the region. Site selection was restricted by requiring that the location be at least 10 pixels from boundaries between PHR areas and the remaining image in each of the land-cover classifications. This was done to reduce the potential misclassification of the selected locations. The changes in NDVI for PHR and remaining sites were compared, within land-cover categories, using time series analyses to determine if the onset, intensity, or

132

Ecological Applications Vol. 17, No. 1

GREGORY E. GLASS ET AL.

sites. We also expected that SNV should be introduced into sites closest to PHR areas than sites farther away. RESULTS

PLATE 1. Persistently high risk area associated with Hantavirus Pulmonary Syndrome. Photo credit: Douglas E. Norris.

duration of vegetation greenness differed. Seasonal Autoregressive Integrated Moving Average (SARIMA) models were fit to the data. Differences between PHR and lower risk sites within strata were identified by the need to postulate alternative models to characterize the comparable time series, such as different model structures or significant differences in the coefficients and lag structures for individual models. Yearly changes in abundance of local deer mouse (P. maniculatus) populations near PHR sites were examined using data from the 38 locations which had been previously sampled for small rodent populations and tested for SNV during 1998 and 1999. Approximately half the sites were in high-risk areas during two years of the study, while the remaining locations were either in low-risk areas or sites that varied in risk during the two year study. Counts of deer mice trapped using standard protocols and the presence of mice with antibodies to SNV, as tested by ELISA, were used as measures of population abundance and the presence of virus. Methods of sampling and testing are described in detail elsewhere (Glass et al. 2002). We expected that if PHR sites were refugia, that populations nearer PHR sites were more likely to retain resident populations of deer mice and were more likely to be recolonized than sites in lower quality habitat that were further away from PHR

Persistently high risk (PHR) sites represented a small fraction of the study area, accounting for 0.3% (311 km2 out of 105 200 km2) of the region. Despite their limited geographic extent, PHR areas had a coherent structure related to physical features, elevation, land cover, and seasonal vegetative growth patterns (see Plate 1). Sites that had previously been sampled for deer mice near PHR locations were more likely, over the two year period, to acquire and retain deer mouse populations that were infected with SNV. PHR sites were spatially stable across years. Sites that were at highest risk during three, four, or five years of the study were approximately 138, 1300, and 310 000 times more likely, respectively, than expected than if the highest risk pixels were randomly distributed (v2 ¼ 210.08, df ¼ 5, P , 0.001) among years (Fig. 1). These sites frequently corresponded to well-defined physical features on the digital raster graphics (DRGs). Often, they occurred along arroyos or canyons associated with mountains and mesas (Fig. 2). However, even within the same canyon, only portions were PHR areas. These sites occurred adjacent to streams, rivers, or other water bodies. Near these features, they occurred at slightly higher elevations along the drainages rather than immediately adjacent to the watercourses themselves. PHR sites were predominantly restricted to a narrow range of elevations compared to the study region as a whole (Komolgorov-Smirnov two-sample test ¼ 0.76; P , 0.0001). PHR sites did not occur below 1600 m elevation, even though more than one third of the study region was lower (Fig. 3). Approximately 80% of the PHR sites were within a narrow 400 m elevation band between 2600 m and 3000 m, although only 6% of the study region occurred within this range. PHR sites also were not found above 3200 m, even though parts of the study region reached 3600 m. In exploratory analyses, PHR sites were more frequent in areas with east

FIG. 1. Temporal persistence of highest risk sites. The histogram shows observed number of pixels rated as the highest risk category/expected number of highest risk pixels if they were under random distribution. Sites classified in the highest risk category tended to persist across years.

January 2007

HANTAVIRUS REFUGIA

133

FIG. 2. PHR sites (red circles) were frequently associated with arroyos and canyons on mesas or mountains (panel A; DRG NJ 12-11, Series V502 quadrangle) or along streams and rivers (panel B; DRG NJ 12-12, Series V502). Locations are indicated by Digital Raster Graphics (DRG) identifiers assigned by the U.S. Geological Survey.

northeast–east southeast aspects than the study region as a whole (16.4% vs. 12.2%) (Watson’s U2 ¼ 5.78, P , 0.001) and they were less common on flat areas than expected (27.8% vs. 40.9%). This, in part, reflected their absence in lower elevation areas, where much of the flat ground occurred. PHR areas were disproportionately associated with evergreen forest classification in the NLCD and negatively associated with shrubland; these two classes

accounted for ;80% of the study area (v2 ¼ 82.45, df ¼ 1, P
0.0001). Most of the remaining PHR areas were grassland or pasture, though in proportions consistent with their occurrence in the region (Table 1). Although more than two-thirds of PHR areas were associated with evergreen forest, this still represented only 0.65% of all the land with this classification (see Plate 1). In contrast, some land-cover classes that occupied vanishingly small portions of the region, such as herbaceous wetlands

134

Ecological Applications Vol. 17, No. 1

GREGORY E. GLASS ET AL.

FIG. 3. Comparison of elevations at PHR sites (persistently highest risk) and those within the study region (background elevation). PHR sites occurred within a narrow range of elevations. Approximately 80% of the PHR areas occurred between 2600 and 3000 m, while only 6% of the study area was at these elevations.

(accounting for 0.004% of the region), were frequently (18.3%) classified as PHR. In multivariate analyses, the number of years a site was in the highest risk level was related to a combination of land cover, elevation, and slope (Table 2). Referent sites were agricultural or pasture land, below 1801 m elevation on flat ground. Bare ground and shrubland were associated with significantly fewer numbers of years at highest risk, compared to referent sites, while mixed and deciduous forest classified areas were associated with significantly more years as highest risk. Evergreen forest was no more likely to contribute additional years at highest risk than the referent landcover class. Increasing elevation was associated with increasing years at highest risk, with the effect markedly increasing above 2130 m (Table 2). Similarly, compared to flat ground, moderate and steep slopes (.48) were associated with more years at highest risk. After controlling for land cover, elevation, and slope, aspect did not contribute significantly to the numbers of years at highest risk (Table 2). Eight to twelve sites, in each of the four major landcover classes (evergreen forest, shrubland, grassland, or pasture) had NDVI values extracted from MODIS NDVI for 848 days from January 2002 through April 2004 (Fig. 4A–D). Regardless of land cover, PHR sites showed a marked seasonal change in NDVI values with a rapid increase beginning about day 64 of each year (Fig. 4, times 4, 29, and 53). NDVI values in PHR sites peaked over a broad time period, between day 160 and 320 of each year, and then declined. In contrast, lower risk sites showed a slower greening through the growing season and only occasionally reached NDVI values approaching those in PHR areas (Fig. 4; Appendix A: Table A1).

Autocorrelation and partial autocorrelation plots were used to identify appropriate models for NDVI each category of sites (Appendix B: Table B1). The time series were differenced to achieve stationarity. The most common class of models was seasonal-mixed-autoregressive moving-average models in which the autoregressive and moving-average components were first order and the seasonal-moving-average component also was first order, with a seasonal component between 352 and 384 days. The best models fitting the seasonal changes in NDVI at PHR sites in evergeen forests, shrublands, and grasslands differed substantially from corresponding sites that were not highest risk (Appendix B: Table B1). Human-modified pastureland sites, whether PHR or not, (Fig. 4D; Appendix B: Table B1), were most similar to one another with similar AR(1) structure, though they differed in their seasonal MA coefficients. Sites sampled for small mammals and tested for SNV in 1998 and 1999 were examined relative to PHR areas. None of the previously trapped sites were within any PHR areas, although 13 of the trapped locations were within a kilometer, and three were within 200 m (Fig. 5). On average, trapped high-risk sites were nearer (2.0 km 6 2.9 km; mean 6 SD) than trapped low-risk sites (13.5 km 6 8.0 km) to PHR areas. No deer mice were trapped at seven sites in 1998. Four sites were low-risk and three were high-risk areas (Fig. 5). During resampling in 1999, only a single deer mouse was captured at one of these four low-risk sites and it was uninfected with SNV, indicating colonization by mice and introduction of SNV was rare at these locations. In addition, seven of the 11 low-risk sites that had deer mice in 1998 lost their mice by 1999 (Fig. 5), indicative of high local extinction rates in these habitats. By 1999, 10 out of 18 low-risk sites were unoccupied. Infected deer mice were found at only one of the low-risk sites. Each of the three unoccupied high-risk sites in 1998 yielded deer mice in 1999, and at least one mouse at each TABLE 1. Comparison of persistently highest risk sites and the remaining portion of the study area by major land-cover classes.

Land cover

Persistently highest risk (%)

Not highest risk (%)

Evergreen forest Shrubland Grassland Pasture Other 

69.04 18.22 6.74 1.00 5.00

27.71 51.95 17.75 0.83 1.76

Notes: Values represent the percentage of the study area that fell into a risk category, stratified by land cover. All classes with .0.5% coverage were categorized separately. Land-cover data were based on the National Land Cover Data (NLCD) 1992 database.   Other: Open water, low-density residential, commercial, bare rock, quarries, deciduous forest, mixed forest, orchards, row crops, small grains, parks, woody wetlands, herbaceous wetlands.

January 2007

HANTAVIRUS REFUGIA

135

TABLE 2. Effect of selected environmental variables on the number of years pixels were identified as persistently highest risk.

Variable

Referent class

Regression coefficient

95% interval estimate

Residential/industrial Bare ground/quarries  Deciduous/mixed forest  Evergreen Shrubland  Wetlands Elevation (1801–1900 m) Elevation (1901–2130 m) Elevation (2131–2800 m)  Elevation (.2800 m)  Slope (18–48) Slope (58–258)  Slope (.258)  Aspect (3168–3608) Aspect (18–448) Aspect (458–898) Aspect (908–1348) Aspect (1358–1798) Aspect (1808–2248) Aspect (2258–2698) Aspect (2708–3158)

agriculture/pasture agriculture/pasture agriculture/pasture agriculture/pasture agriculture/pasture agriculture/pasture elevation  1800 m elevation  1800 m elevation  1800 m elevation  1800 m flat ground flat ground flat ground flat ground flat ground flat ground flat ground flat ground flat ground flat ground flat ground

1.09 0.72 0.48 0.02 0.45 0.76 0.19 0.10 3.43 4.21 0.04 0.38 0.51 0.04 0.10 0.14 0.10 0.05 0.01 ;0.03 0.00

15.56, 0.93 1.41, 0.09 0.26, 0.68 0.13, 0.12 0.61, 0.30 14.27, 1.09 0.60, 0.10 0.21, 0.44 3.27, 3.65 4.04, 4.43 0.023, 0.37 0.16, 0.68 0.22, 0.84 0.039, 0.24 0.27, 0.34 0.22, 0.38 0.26, 0.38 0.44, 0.25 0.38, 0.25 0.36, 0.25 0.29, 0.24

Notes: Multivariate Poisson regression was used to control for confounding effects by estimating coefficients from a randomly selected 11 000 pixels in the region. The models were repeated 200 times, and empirical point and interval estimates were derived as the median value for each variable and the 95% interval. Variables that differ significantly from zero (as derived by Monte Carlo estimates; see Methods: Multivariate analysis) are indicated by daggers ( ). All changes are in comparison to a selected category (referent class).

of these sites had antibodies to SNV (Fig. 5). All of the high-risk sites had deer mice by 1999 and 18 out of 20 had at least one deer mouse infected with SNV, suggesting successful population persistence and that viral introduction was common in suitable habitats near PHR sites. Sites with infected deer mice were, with one exception, within 8 km of PHR areas. Sites at farther distances either no longer had deer mice by the second year, or those that were present were not infected. DISCUSSION Sites with persistent environmental signatures similar to where people were at highest risk for HPS occurred in very limited locations within a geographic region that is viewed as a major epicenter for disease. If PHR sites represent refugia, their limited extent could explain previous workers’ difficulty locating them. PHR sites had well-defined physical characteristics and were focused in a narrow elevational range with moderate to steep slopes. Land-cover categories were not good descriptors of increasing risk, though this may reflect the habitat generalist nature of the rodent reservoir or the coarse categorization of land classes in the data set. Deer mice occur in diverse habitats throughout much of North America, and in the U.S. Southwest, they range from desert to alpine habitats (Mills et al. 1997). Despite this wide range of tolerances, Mills and colleagues found that the abundance of P. maniculatus varied with altitude and vegetation community, as did the prevalence of Sin Nombre Virus (SNV) in the local

rodent populations. They found deer mice less commonly at the lowest and highest elevations and prevalence of SNV decreased at the elevation extremes. Nearly half (20/41 sites) of their locations where deer mice occurred yielded no mice with antibodies to SNV. Occasionally this was due to small sample sizes, but many of their sites had either a disproportionate overabundance or paucity of mice with SNV. The focal nature of infection was unrelated to population density and the bases for the patterns remained unexplained. However, their results mirror the altitudinal and land cover distributions of PHR sites, with a disproportionate fraction of sites above 2100 m associated with a variety of land-cover categories. PHR and adjacent areas also may account for the focal distribution of SNV infection in mouse populations, as areas adjacent to PHR sites were more likely to acquire and retain infected deer mice than sites at greater distances (Fig. 5). Similarly, Calisher and colleagues (2001) have hypothesized that local edaphic and vegetation conditions drive reservoir population dynamics to favor either an excess or paucity of infected population members by affecting rates of reproduction and altering age structure so that in some sites long-lived individuals serve as the sources for SNV infection and maintain transmission. The primary advantage of the current analysis is that it provides the detailed characterization of these sites, as well as their geographic locations. A key factor distinguishing PHR sites from other sites was the seasonal pattern of vegetative growth (Fig. 4).

136

GREGORY E. GLASS ET AL.

Ecological Applications Vol. 17, No. 1

FIG. 4. Seasonal changes in normalized difference vegetation index (NDVI) derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) between 2002 and the first quarter of 2004 at selected PHR sites (dotted lines) and sites that were not PHR (solid lines). Colors of lines were assigned to individual locations arbitrarily. Sites are stratified by land cover: (A) evergreen forest, (B) shrubland, (C) grassland, and (D) pasture. These four land cover categories accounted for .98% of the study area. PHR sites tended to show an earlier, more rapid and more extensive vegetative growth than did comparable lower risk sites.

Even during an extensive drought, PHR sites had an early season onset of vegetation growth (early March), with an extensive growing season (until late October to early November), before a rapid decline with the onset of winter conditions, typically in mid- to late November. The ability of vegetation in PHR sites to show such vigorous and long-lasting growing seasons often corresponded to locations along arroyos, drainage systems where moisture is likely to accumulate and persist. The association with increased moisture also is consistent with the increase in years at risk with deciduous and mixed deciduous land cover (Table 2) which, in this

region, is predominantly associated with oaks (Quercus spp). Thus, although the region is generally viewed as high risk for HPS, the real area of consistent, long-term, high risk is quite limited. The specific characteristics of these sites may explain why during most years there are few human cases, while during some years larger numbers of HPS are observed. The increased number of years at risk associated with sites .2130 m (Table 2) occurs near the maximum elevation occupied by most people in the region. The altitudinal distribution of the human control sites observed by Glass et al. (2000) was 2370 m (1962 6

January 2007

HANTAVIRUS REFUGIA

137

FIG. 4. Continued.

198 m; mean 6 SD). PHR sites, in particular, were predominantly within the 2600–3000 m band (Fig. 3), above the elevation of human residence. Because few people reside at elevations with PHR sites, relatively few are likely to become infected during most years. If PHR areas are refugia, then during most years mice would need to disperse from these locations to lower elevations before reaching human habitations. Distance from PHR sites might be a surrogate of human risk and the timing of pathogen acquisition. The current analysis does not establish that refugia for SNV exist or that the PHR areas have identified them. However, it shows that areas with environmental signatures similar to those that distinguish high-risk

locations for HPS from lower risk sites may persist for many years. These sites may represent refugia because the physical and ecological conditions are consistent with environmental features that favor survival and reproduction in deer mouse populations (Yates et al. 2002). These sites also have characteristics similar to those favoring SNV persistence (Childs et al. 1995, Mills et al. 1997). The observation that many of the highest risk sites remained as such for four or five of the years studied (Fig. 1) suggests that they could serve as sites where SNV reaches high levels in local reservoir populations. A single year as a high-risk area is not sufficient for SNV to reach high levels in reservoir populations, but within

138

GREGORY E. GLASS ET AL.

Ecological Applications Vol. 17, No. 1

FIG. 5. Number of deer mice captured in 1998 (dark bars) and 1999 (light bars) with increasing distance (km) from PHR sites. Classified high-risk sites (indicated by daggers above bars) tended to be closer to PHR sites than low-risk sites. Low-risk sites were more likely to lose deer mice between 1998 and 1999.

two years as a high-risk site, the average prevalence of infection in local populations can reach more than 30% (Glass et al. 2002), especially near PHR areas. This is consistent with the original epidemiologic investigations of HPS cases (Childs et al. 1995) where, in case households, the prevalence of SNV infection among mice was more than 30%. By implication, the Landsat image classifier identifies environmental conditions that are associated with large deer mouse populations that have high prevalence of SNV. If PHR sites are refugia, then we expect that sampling in locations with these environmental signatures will yield deer mice at times when sampling in other locations finds few if any of this species. In addition, local populations of deer mice sampled in these sites should be more likely infected with SNV than samples of mice collected in sites that are not persistently highest risk locations. The loss of deer mice from low-risk sites in 1999 (Fig. 5) and the colonization of high-risk sites near PHR areas, with the introduction of SNV, is consistent with this hypothesis. If the trophic cascade hypothesis (TCH) is true, then sampling with increasing distance from PHR should be able to document the movement of individual mice from PHR to lower elevations as environmental conditions become more favorable and net primary productivity (NPP) increases. According to the TCH, the movement of mice from refugia should conform to metapopulation theory (Hanski 1999), and we expect that locations away from refugia should be colonized as a decreasing function of distance from the refugia. Local extinction among these colonized sites should be inversely related to local

habitat quality and patch size and should occur as environmental conditions become more severe. Further detailed studies are needed to characterize the ecological dynamics of PHR areas, compared to similar areas that are not PHR, to better understand the dynamics that influence deer mouse populations. Additional studies also are needed to evaluate the quality of the environmental data layers used in these analyses. A primary advantage of using these environmental data is the well-characterized metadata. However, groundtruthing is needed to clarify some issues. For example, the PHR areas in NLCD-classified (National Land Cover Data) pasture and grassland, though representing small geographic regions, were unexpected because these are not usually considered suitable habitats for deer mice. This could indicate misclassification in the landcover database. However, Mills and colleagues (1997) reported deer mice from some of their sites in these habitats and Parmenter (R. Parmenter, unpublished data) has observed relatively large local populations during the autumn in these habitats in northern New Mexico. Future ground-based studies could target these apparently anomalous locations to determine the reasons for the observations and evaluate if these locations truly are suitable for reservoir species. Regardless, combining remotely sensed environmental data with ground studies provides a rational strategy for targeted field studies of local populations of deer mice that are responsible for transmission of SNV. By identifying the physical structure and spatial context of PHR locations, we can use the geographic extent and contiguity of local populations to create testable

January 2007

HANTAVIRUS REFUGIA

hypotheses concerning the persistence and movement of individuals responsible for viral transmission. ACKNOWLEDGMENTS We thank the many individuals who assisted with the field aspects of the study, especially those individuals who conducted the sampling efforts in 1998 and 1999. Financial support for this study was provided by funding from the NSF and NIH program in Ecology of Infectious Diseases DEB-0326757, NASA NCC5-305, and NOAA NA96GP0419. Comments from M. L. Farnsworth and an anonymous reviewer were especially helpful and improved the quality of the final manuscript. LITERATURE CITED Calisher, C. H., J. N. Mills, W. P. Sweeney, J. R. Choate, D. E. Sharp, K. M. Canestorp, and B. J. Beaty. 2001. Do unusual site-specific population dynamics of rodent reservoirs provide clues to the natural history of hantaviruses? Journal of Wildlife Diseases 37:280–288. Childs, J. E., et al. 1994. Serological and genetic identification of Peromyscus maniculatus as the primary rodent reservoir for a new hantavirus in the southwestern United States. Journal of Infectious Diseases 169:1271–1280. Childs, J. E., et al. 1995. A household-based, case-control study of environmental factors associated with hantavirus pulmonary syndrome in the southwestern United States. American Journal of Tropical Medicine and Hygiene 52:393–397. Glass, G. E., et al. 2000. Anticipating risk areas for hantavirus pulmonary syndrome with remotely sensed data: re-examination of the 1993 outbreak. Emerging Infectious Diseases 6: 238–247. Glass, G. E., et al. 2002. Satellite imagery characterizes local animal reservoir populations of Sin Nombre virus in

139

Southwestern United States. Proceedings of the National Academy of Sciences (USA) 99:16817–16822. Hanski, I. 1999. Metapopulation ecology. Oxford University Press, Oxford, UK. Hjelle, B., and G. E. Glass. 2000. Outbreak of hantavirus infection in the Four Corners Region of the US in the wake of the 1997–98 El Nin˜o Southern Oscillation. Journal of Infectious Diseases 181:1568–1573. Lillesand, T. M., R. W. Kieffer, and J. W. Chipman. 2004. Remote sensing and image interpretation. Fifth edition. John Wiley and Sons, Hoboken, New Jersey, USA. Mills, J. N., et al. 1997. Patterns of association with host and habitat: antibody reactive with Sin Nombre virus in small mammals in the major biotic communities of the southwestern United States. American Journal of Tropical Medicine and Hygiene 56:273–284. Nichol, S. T., C. F. Spiropoulou, S. Morzunov, P. E. Rollin, T. G. Ksiazek, H. Feldmann, A. Sanchez, S. Zaki, J. Childs, and C. J. Peters. 1993. Genetic identification of a novel hantavirus associated with an outbreak of acute respiratory illness in the southwestern United States. Science 262:914– 917. Parmenter, R. R., J. W. Brunt, D. I. Moore, and M. S. Ernest. 1993. The Hantavirus epidemic in the Southwest: rodent population dynamics and the implications for transmission of hantavirus-associated adult respiratory distress syndrome (HARDS) in the Four Corners Region. Pages 1–45 in Sevilleta Publication No. 41. Sevilleta Long Term Ecological Research (LTER) Project, University of New Mexico, Albuquerque, New Mexico, USA Yates, T. L., et al. 2002. The ecology and evolutionary history of an emergent disease: Hantavirus Pulmonary Syndrome. BioScience 52:989–998.

APPENDIX A Descriptive statistics for NDVI in persistently highest risk (refugia) and other (not refugia) areas classified as evergreen forest, shrubland, grassland, and pasture (Ecological Archives A017-007-A1).

APPENDIX B Coefficients from SARIMA models for average NDVI values in four major land-cover classes (evergreen forest, grassland, shrubland, pasture) in highest risk areas (refugia) and lower risk areas (not refugia) (Ecological Archives A017-007-A2).