Ecology, 94(5), 2013, pp. 1036–1045 � 2013 by the Ecological Society of America Gradient-based habitat affinities predict species vulnerability to drought DIANE M. DEBINSKI,1,5 JENNET C. CARUTHERS,2 DIANNE COOK,3 JASON CROWLEY,3 AND HADLEY WICKHAM 4 1Department of Ecology, Evolution and Organismal Biology, 251 Bessey Hall, Iowa State University, Ames, Iowa 50011 USA 2Interdepartmental Program in Ecology and Evolutionary Biology, Iowa State University, Ames, Iowa 50011 USA 3Department of Statistics, Iowa State University, Ames, Iowa 50011 USA 4Department of Statistics, Rice University, Houston, Texas 77005 USA Abstract. Ecological fingerprints of climate change are becoming increasingly evident at broad geographical scales as measured by species range shifts and changes in phenology. However, finer-scale species-level responses to environmental fluctuations may also provide an important bellwether of impending future community responses. Here we examined changes in abundance of butterfly species along a hydrological gradient of six montane meadow habitat types in response to drought. Our data collection began prior to the drought, and we were able to track changes for 11 years, of which eight were considered mild to extreme drought conditions. We separated the species into those that had an affinity for hydric vs. xeric habitats. We suspected that drought would favor species with xeric habitat affinities, but that there could be variations in species-level responses along the hydrological gradient. We also suspected that mesic meadows would be most sensitive to drought conditions. Temporal trajectories were modeled for both species groups (hydric vs. xeric affinity) and individual species. Abundances of species with affinity for xeric habitats increased in virtually all meadow types. Conversely, abundances of species with affinity for hydric habitats decreased, particularly in mesic and xeric meadows. Mesic meadows showed the most striking temporal abundance trajectory: Increasing abundances of species with xeric habitat affinity were offset by decreasing or stable abundances of species with hydric habitat affinity. The one counterintuitive finding was that, in some hydric meadows, species with affinity for hydric habitats increased. In these cases, we suspect that decreasing moisture conditions in hydric meadows actually increased habitat suitability because sites near the limit of moisture extremes for some species became more acceptable. Thus, species responses were relatively predictable based upon habitat affinity and habitat location along the hydrological gradient, and mesic meadows showed the highest potential for changes in community composition. The implications of these results are that longer-term changes due to drought could simplify community composition, resulting in prevalence of species tolerant to drying conditions and a loss of species associated with wetter conditions. We contend that this application of gradient analysis could be valuable in assessing species vulnerability of other taxa and ecosystems. Key words: butterflies; climate change; drought; gradient analysis; Greater Yellowstone Ecosystem, USA; temporal trends. INTRODUCTION The science of detecting plant and animal responses to climate change has made major strides at the global scale, particularly by quantifying poleward and eleva- tional range shifts in species distribution patterns (Parmesan and Yohe 2003, Root et al. 2003, Parmesan 2006, Lenoir et al. 2008). However, it is much more difficult to move from global models and analyses of climate change to predicting biotic responses at fine geographic scales. Whereas thermal limitations such as minimum winter temperature are important determi- nants of species range distributions at the continental scale (Root 1988), and minimum or maximum temper- atures may determine range limits along an elevational gradient (e.g., Pounds et al. 2006), additional micro- habitat variables influence species distributions at finer geographic scales (e.g., Dennis and Sparks 2006). Because landscapes are composed of a complex matrix of other gradients, even within a constant elevation, early evidence of climate change may be found by quantifying shifts in abundance along these gradients. Here we describe how an analysis of changes in species distribution and abundance along a hydrological gradi- ent, before and after drought conditions, can inform scientists about potential fine-scale biotic responses to climate change in areas that are expected to become warmer and drier. Manuscript received 9 March 2012; revised 13 December 2012, accepted 18 December 2012. Corresponding Editor: J. T. Cronin. 5 E-mail: debinski@iastate.edu 1036 Debinski, Diane M., Jennet C. Caruthers, Dianne Cook, Jason Crowley and Hadley Wickham. 2013. Gradient- based habitat affinities predict species vulnerability to climate change. Ecology. 94(5): 1036-1045. Available at: http://dx.doi.org/ 10.1890/12-0359.1. Gradient analysis, pioneered by Whittaker (1956, 1972), has been used more recently by plant ecologists to assess niche utilization (Silvertown et al. 1999), to assess topographic gradients in species richness (Fleishman et al. 2000), and to predict local species richness (Grace et al. 2011). Landscapes with strong gradients, and particularly those of moisture and elevation, can have highly predictable distributions of plant and animal communities. Previously we used gradient analysis in these same sites to track changes in plant distribution in response to drought. Plant cover changes were tracked along a hydrological gradient of 55 montane meadows during an 11-year time period (1997–2007) that included eight years of mild to extreme drought (see Debinski et al. 2010 for drought data). Each of the six meadow types within the gradient is composed of a distinct plant community (Debinski et al. 2000, 2006). The beginning of our project (1997–1999) was characterized by above- normal to normal moisture conditions, but 2000–2007 was characterized by mild to extreme drought condi- tions. The plant community responses to drought included increases in bare ground, especially in mesic to xeric meadows, decreases in forb cover especially in mesic to xeric meadows, and increases in woody cover (Debinski et al. 2010). Given the number of significant changes in the plant community, we expected that there also could be corresponding changes within the insect community. Drought stress has been shown to be correlated with local extinctions of butterfly populations in California (Ehrlich et al. 1980) and Canada (Packer 1994), as well as changes in European butterfly distributions (Morecroft et al. 2002). At larger temporal and spatial scales, poleward (Parmesan et al. 1999) and elevational shifts (Forister et al. 2010) in distribution, as well as phenological changes (Forister and Shapiro 2003, Kearney et al. 2010, Diamond et al. 2011), have been documented in butterfly communities in response to climate change. Similarly, Breed et al. (2012) recently documented population size changes at species range margins in the eastern United States. Because butterflies have requirements for specific host plants as larvae and, to a slightly lesser degree, nectar species as adults, they have a tight association with plant communities (Ehrlich and Raven 1964, Gilbert and Smiley 1978, Debinski et al. 2000) and can serve as good indicators of environmental change. We suspected that in a complex landscape, the effects of drought would not be uniformly manifested, but would be predictable based upon landscape-level pat- terns of hydrological gradients and habitat affinities of the insect species associated with these meadows. Our goal was to test whether drought conditions could be correlated with temporal trajectories in butterfly abun- dance. Because the time series started out with wetter conditions followed by several years of drought condi- tions, we expected that there could be long-term (i.e., 11 year) trends in overall butterfly abundance patterns because butterflies could be responding both directly to the drought as well as indirectly to the reduction in host plant and nectar resources. In order to evaluate butterfly community responses, we grouped species into three major categories: (1) species with an affinity for hydric to mesic habitats (hereafter shortened to ‘‘hydric affinity’’), (2) species with an affinity for xeric habitats, and (3) species that could be found in either type of habitat. We expected that hydric-associated butterfly species would be most sensitive to drought conditions and would decrease in abundance during this time period. We expected that xeric meadow butterfly communities would be relatively well adapted to dry conditions, and might even increase in abundance during this time. Because the plant community had shown the most significant losses in forb cover within mesic meadows (Debinski et al. 2010), we also expected that butterflies in mesic meadows would be most vulnerable to change. MATERIALS AND METHODS Study area The relatively pristine nature and minimal human impacts of the study sites within the Greater Yellow- stone Ecosystem, USA, makes it an ideal location for studying how organisms such as butterflies respond to changes in the conditions of the environment. The Greater Yellowstone Ecosystem was divided up into two different study regions, referred to as the Gallatin (including 30 sites) and the Teton (includes 25 sites) regions. The two regions are separated by 192 km, yet both have similar plant and butterfly communities (Su et al. 2004). The meadows selected for the surveys in both regions have relatively homogenous topographic features. The average elevation of sites is 2098 m in the Gallatin region, and 2120 m in the Teton region, and the average meadow patch size was 360 ha (see Debinski et al. 2001 and 2006 for additional details on study sites). All sites were selected at approximately the same elevation to effectively hold elevation constant. Sites were considered suitable for sampling if they were at least 1003 100 m in size and a minimum distance of 500 m from other sites. Although there is some variance in meadow area between regions, meadow area did not have any significant effect on the butterfly counts in this data set (D. M. Debinski and D. Cook, unpublished data). Six meadow types (M1–M6) were characterized along a hydrologic gradient using satellite imagery (Debinski et al. 2000, 2001). The Gallatin region has five replicates of each meadow type from M1 and M2 (hydric; dominated by Salix spp. and sedge with sparse forbs), M3 and M4 (mesic; dominated by forbs and grasses), M5 and M6 (xeric; dominated by Artemisia spp. with some forbs), and the Teton region has five replicates of each meadow type except meadows characterized as M4. For a map of the study site locations see Saveraid et al. (2001). May 2013 1037GRADIENT-BASED VULNERABILITY ANALYSIS Butterfly surveys Field surveys were conducted primarily during June and July for two-week periods at each region, alternat- ing between the two regions, with two surveys for each region. Surveys were conducted in the Teton region in 1997–1998, 2000–2001, and 2003–2007. The Gallatin region was surveyed in 1997–1998, 2000–2001, and 2006–2007. Sites were located in the field using GPS coordinates. Within each meadow a direction from the center stake was randomly selected in 1997 for the placement of a 50 3 50 m plot, and the same quadrant was surveyed annually. Surveys were conducted on sunny days between 10:00–16:30 hours when the temperature was above 218C with low to moderate wind. Butterflies were surveyed by two observers walking through the plot and capturing as many butterflies as possible within the 20-min timeframe using aerial nets. Sampling did not include skippers. Butter- flies were stored in glassine envelopes until the end of the survey when they were tallied by species and then released. Species that were difficult to identify were taken back to the laboratory as vouchers or photo- graphed. Butterfly taxonomy was based upon the North American Butterfly Association taxonomy (NABA 2001). The total annual abundance for each species in each site was defined as the sum of the abundance across the two sampling times. In 2007, of the 25 Teton sites, bear closures prevented us from sampling three M1 and three M2 meadows. We statistically corrected for this smaller sampling effort in the analysis by scaling the individual species counts proportionally to estimate abundance as if all five sites of each meadow type had been surveyed based on the counts from the site of previous years with the final estimated count rounded to the nearest integer. Data analysis We used a combination of three different approaches to assess abundance changes in the butterfly community over time: ordination, habitat affinity-based trend analyses, and single-species-based trend analysis. But- terfly species with a total abundance ,10 individuals over all years for either region were eliminated from each of the analyses. In addition, species with a maximum abundance of four individuals or less (across all meadow types of the same region) were not included in Table 1, or in any statistical test, in order to minimize inclusion of species with such small abundances that a statistical trend might not also be associated with a biological trend. Although we would have ideally conducted an autocorrelation analysis, the data do not lend themselves to spatiotemporal autocorrelation mod- eling because they are subdivided into six different meadow types and two regions, and we did not have data consistently for each of the years. However, our criteria for site selection helped to minimize possible autocorrelation in the data. We had stringent rules about the selection of the sites when we started the experiment that minimized issues of spatial autocorre- lation. These rules are described in Debinski et al. (2000, 2001). Ordination Ordination was used to visually describe the relation- ship between butterfly species composition and meadow type across time from 1997 to 2007. The difference between these two points was characterized by a vector from time t ¼ 0 (1997) to time t ¼ 10 (2007), whose length was associated with the amount of change exhibited during this time period. We selected nonmetric multidimensional scaling (NMDS), an unconstrained ordination technique using the metaMDS function in TABLE 1. Summary of temporal trends by species based upon habitat-based trend analysis models (see Appendix C for more detailed description of results). Habitat affinity Species Increasing Decreasing Stable Hydric Lycaena editha (15), Lycaena helloides (60), Speyeria atlantis (14), Colias pelidne (5), Colias gigantia (20), Glaucopsyche piasus (5), Lycaeides idas (12) Coenonympha haydenii (100), Plebejus saepiolus (80), Boloria kriemheld (15), Glaucopsyche lygdamus (40), Lycaena nivalis (7), Erebia epipsodea (80), Pieris napi morph 1(14), Pieris napi morph 2 (10), Speyeria cybele (15), Lycaena hyllus (15), Colias interior (14) Coenonympha tulia inornata (100) Xeric Speyeria egleis (7), Speyeria callipe (15), Speyeria zerene (8), Euphydryas editha (8), Euphilotes enoptes ancilla (10), Euphydryas chalcedona (20), Cercyonis oetus (80), Plebejus lupinus (40), Lycaena heteronea (70), Plebejus shasta (6), Lycaeides melissa (12), Chlosyne palla (12), Parnassius clodius (12) Plebejus icariodes (120), Callophrys sheridanii (10), Parnassius smintheus phoebus (12), Callophrys dumetorum (10) Oeneis chryxus chryxus (8) Notes: Species included in this list exhibited significant trends. An increasing trend was defined where species showed increasing trends predominantly across all meadow types and both regions (11 separate graphs), and a decreasing trend showed an inverse pattern. The value in parentheses next to each species is the maximum abundance at which it occurred in any one year. DIANE M. DEBINSKI ET AL.1038 Ecology, Vol. 94, No. 5 the vegan package (Oksanen et al. 2008) of the statistical software R (R Development Core Team 2010). NMDS configures the points observed in the ordination plot based on species abundance data, and was plotted using the qplot function in the package ggplot2 (Wickham 2009). We used the Bray-Curtis distance metric to calculate the community similarities between meadow types for each year and region. Data plotted in the ordination represent the summation of butterfly abundance by species across all meadow replicates for a year within one meadow type and region. The Bray-Curtis index compares a pair of points (sites) to see how many species they share in common compared to all the species that differ between them (Jongman et al. 1995, Summerville and Crist 2003). Abundance for each species also weighs into this distance value. Points closer together on the ordination plot have more similar species compositions than those farther apart. Habitat affinity-based trend analysis In order to detect whether species with particular habitat affinities were responding similarly to environ- mental changes, we conducted analyses of groups of species that were associated with a particular habitat type. This habitat affinity-based trend analysis was conducted by grouping species into one of two major groups: a xeric-associated class or a hydric-associated class, and evaluating the trend for the entire group of species (see Appendix A: Table A1, Butterfly species categorized by habitat affinity). When conducting the habitat affinity-based trend analysis, we used the abundance of each single species as a data point and accounted for this pseudoreplication by having species as a random effect and evaluating the trend for the groups via an interaction term of year and group. A small number of species had no specific affinities for either type of meadow. These species were not included in the habitat affinity-based model, but were studied individually using single-species models. Species catego- rization was based upon host plant and habitat requirements as described in field guides (Ferris and Brown 1981, Scott 1986, Bird et al. 1995, Debinski and Pritchard 2002). There were a few species that occurred in one region but not in the other. These species were included in the habitat affinity-based trend analysis by inserting a zero for their abundance value in the region where they were not observed. To model habitat affinity-based trends, a mixed- effects Poisson regression was utilized. It is natural to assume that the response variable (count) follows the Poisson distribution, and a mixed-effects model allows the differences in counts between species to be accommodated by a random intercept. Year, region, and meadow type were treated as fixed effects, and species was included as a random effect. The fixed effects part of the model can be summarized as follows: logeðYijklÞ ¼ b0 þ b1ðYEARÞ þ aiðREGION ¼ GallatinÞ þ cjðMTYPE ¼ M1Þ þ dkðHABITAT AFFINITY 5 HydricÞ þ b1aiðYEAR 3 REGIONÞ þ b1cjðYEAR 3 MTYPEÞ þ b1dkðYEAR 3 HABITAT AFFINITYÞ þ cjdkðMTYPE 3 HABITAT AFFINITYÞ þ eijkl; i ¼ 1; 2 ðGallatin;TetonÞ; j ¼ 1; . . . ; 6 ðM1;M2;M3;M4;M5;M6Þ; k ¼ 1; 2 ðHydric;XericÞ; l ¼ 1; . . . ; nijk The R package lme4, and function lmer (Pinheiro and Bates 2009) was used to fit the mixed-effects model. Single-species-based trend analysis A Poisson regression model was fit separately for each species, which is what we call the single-species model. This has the benefit of allowing more freedom for the model for each species, so we can carefully examine the trends without the constraints of the trends in other species. The single-species models also allowed for test of significance at the species level, whereas the habitat- based affinity model did not. In both the habitat-based affinity and single-species models, we tested for linear as well as curvilinear patterns. We plotted trends as raw abundance values for the single-species models. We plotted trends on a log scale for the habitat-based affinity model because we were evaluating multiple species concurrently and many of the abundance values are small. This allowed us to see trends more clearly across species within a group. RESULTS Butterfly community A total of 87 butterfly species (including a few subspecies separations) were observed across the two regions, with 75–76 species observed in each region and 11–12 species unique to each region (Appendix A: Table A1). After accounting for minimal abundance values, 60 species were used in the ordination, 50 species were used in the habitat affinity-based (mixed-effects model), and 60 species were used in the single-species models. The sections below summarize results of the ordination, habitat affinity-based analyses, and single-species models. Additional details are provided in Appendix B (mixed- effects model) and Appendix C (single-species models). Ordination Two dimensions represented 16% of the variation between sites, showing the main differentiation, mois- May 2013 1039GRADIENT-BASED VULNERABILITY ANALYSIS ture, as the variable of interest driving butterfly community composition. The Teton region had more distinct clustering of sites within each meadow type, but less consistency in temporal trajectories among meadow types (Fig. 1a). The Gallatin region showed less distinction of butterfly communities among meadow types, but the temporal trajectory was consistent. In both regions, M4 and M5 meadows exhibited the least amount of change. Species locations in the ordination space are depicted in Fig. 1b, and species are color-coded by habitat affinity. Species clearly separate themselves by habitat affinity in axis 1 of Fig. 1b; species with hydric affinities assemble on the left and species with xeric affinities assemble on the right. Axis 2 is associated with regional affinity. Species with low scores on axis 2 (Boloria selene, Boloria frigga, Speyeria cybele, Colias interior, Cal- lophrys sheridanii, and Callophrys dumetorum) are much FIG. 1. Nonmetric multidimensional scaling (NMDS) ordination plot of butterfly survey sites by meadow type and region displayed via (a) arrows depicting temporal trend from 1997 to 2007. Meadow types are M1 and M2 (hydric; dominated by Salix spp. and sedge with sparse forbs), M3 and M4 (mesic; dominated by forbs and grasses), M5 and M6 (xeric; dominated by Artemisia spp. with some forbs), and regions are the Gallatin and Teton, in the Greater Yellowstone Ecosystem, USA. Color-coded arrows depict meadow-specific temporal trends by region from 1997 to 2007. The length of the arrow corresponds to the significance of the temporal trend. (b) The distribution of species by habitat affinity (hydric, xeric, and either) superimposed in the same NMDS space. The center of the text string represents the location of the point. Axis 1 of the ordinations represents the variation in butterfly community composition along the hydrologic gradient (M1 on the left to M6 on the right of zero), while Axis 2 represents the difference in butterfly species composition between the Gallatin (upper) and Teton (lower) regions. DIANE M. DEBINSKI ET AL.1040 Ecology, Vol. 94, No. 5 more numerous in the Teton region, and species with high scores on axis 2 (Speyeria atlantis hesperis, Phyciodes campestris, and Coenonympha haydenii [see Plate 1]) are much more numerous in the Gallatin region. Habitat affinity-based trend analysis Across both regions and all meadow types, xeric species increased over time (fixed-effects fit of region, meadow, and habitat affinity; Fig. 2). Hydric species decreased consistently across regions in the mesic and xeric meadows, but the M1 and M2 (hydric meadows) showed increases for hydric species, particularly in the Teton region. The Tetons also exhibited a larger increase in counts by year than the Gallatins. Interestingly, the temporal trend for M3 meadows in both regions involved an intersection of the abundance curves, caused by an increase in the xeric species that was offset by a decrease (Gallatin) or stable trend (Teton) in the hydric species. The habitat affinity-based trend model is a good fit to the data, explaining most of the variation in the butterfly counts. The model deviance was 10 717 reduced from a null deviance of 451 331. The model explains 97% of the variation in count, although the fixed effects part explains only 4% of the total variation (this represents a pseudo R2). This is to be expected because the main source of variation in the data is the species to species difference. Despite the small amount of variation explained, the fixed-effects model is statistically significant, and the evidence is strong that there are different trends for the two habitat affinity groups. Detailed results of the mixed-effects model, including estimates of factor effects, are included in Appendix B. Single-species models were plotted for the species with the top five increasing and decreasing trends by meadow type, region, and habitat affinity (Fig. 3). Notably, some of the most abundant species with affinities for xeric habitats showed the most striking increases over time, and this was particularly evident in xeric meadows. Table 1 summarizes species responses, as distilled from the single-species models, by habitat affinity and trend, and describes the diverging pattern of the two groups of butterflies. Of the species that showed a predominantly upward trend, 13 of these species had xeric affinities and only 7 had hydric affinities. Of the species that showed a predominantly downward trend, 12 of these species had hydric affinities and only 4 had xeric affinities. This pattern is statistically significant (v2 ¼ 5.09, P ¼ 0.024), giving further evidence of a relationship between trend and affinity (i.e., xeric species increased, while hydric species decreased). Plots of individual species trends are available in Appendices B and C for both the mixed- effects model and the single-species models. The single-species models and the ordination corrob- orated the habitat affinity-based model. The trends in species abundances were significant and were generally consistent with our hypotheses of xeric species increas- ing and hydric species decreasing (but there were some exceptions for hydric species in hydric meadows). We chose to emphasize the habitat affinity-based model for its simplicity and comprehensive incorporation of the information. However, it is important to emphasize that the single-species models allow for tests of significance at the species level. And for some species, the fit of the habitat affinity-based model was not as good as the single-species models, because in order to fit the full model, the region and meadow were constrained to have additive effects (e.g., see the single-species model for Speyeria mormonia in Appendix C). The single-species model should be used for making statements about the trend in these cases. DISCUSSION Our results indicate that across the community, butterflies varied in the way that their abundances changed during drought conditions, but that these responses were predictable based upon habitat affinity. Species with xeric affinities increased across all meadow FIG. 2. Mixed-effects model showing the temporal abundance trajectories for xeric and hydric meadow affinity butterflies shown separately by region and meadow type. Note that abundance was log-transformed and plotted on a log axis. May 2013 1041GRADIENT-BASED VULNERABILITY ANALYSIS types in both the Gallatin and Teton regions, and the increases were particularly striking at the xeric end of the gradient. Species with hydric affinities showed decreases in abundance in mesic and xeric meadows (Fig. 2). Mesic meadows, the meadow types we expected to be most vulnerable, showed a particularly intriguing temporal trajectory, starting with higher abundance of hydric species and ending with a higher abundance of xeric affinity species. These results concur with Aldridge et al. (2011), who found that the flowering-plant community in an alpine meadow shifted in floral resources during a drought and that the mesic meadows were most vulnerable to change. The implications of our data are that longer-term changes due to drought could simplify butterfly community composition, resulting in prevalence of a species tolerant to drying conditions and a loss of species associated with wetter conditions. Morecroft et al. (2002) tested similar hypotheses regarding butterfly responses to drought, but their response variable were slightly different than ours. They found that the most common species increased in abundance during drought, whereas species with low mobility showed decreases in abundance. Breed et al. (2012) studied population trends of butterflies in the northeastern United States from 1992 to 2010 and found that traits such as overwintering as eggs or unfed neonate larvae were strongly associated with declines. They suggest that such species would be more susceptible to dehydration if the climate becomes warmer and dryer. Assessment of additional life history traits on these patterns in our system may provide additional insight into these responses. The second community-level trend we observed from 1997 to 2007 was a differential shift in community composition between butterflies in the Teton region and those within the Gallatin region. Whereas Gallatin communities exhibited a consistent temporal trajectory, Teton region communities did not. Our results concur with results observed in British butterflies (Morecroft et al. 2002, González-Megı́as et al. 2008), where changes in community composition differed according to the habitat requirements of the species and their previous distributions. They also concur with Walther et al. (2002), who described how climate changes could be associated with altered butterfly community composi- FIG. 3. Individual species plots of abundance by year separately by meadow for xeric and hydric meadow affinity butterflies. Species with the top five increasing trends are shown in the left column, and species with the top five decreasing trends are shown in the right column. The random effects model for each species is overlaid on the data, and the models for both regions (solid lines and circles for Gallatin, dashed lines and triangles shows Teton) are displayed together. Note that scales vary in these plots based upon maximum abundance values for the species. DIANE M. DEBINSKI ET AL.1042 Ecology, Vol. 94, No. 5 tion due to the variability in rates at which species shift their ranges. Finally, they support our previous findings that Gallatin meadows were more consistent in their responses to environmental change than Teton meadows (Debinski et al. 2006). This difference could indicate potential indirect effects via changes in the plant communities in addition to climate. The one surprising result was found when we examined species responses at a finer scale along the hydrological gradient. Species with hydric affinities showed increases in hydric meadows, particularly in the Teton region. This result, however, may not be as surprising when we investigate the species driving the change. Of the seven hydric species showing significant increases (Table 1), none of them were restricted to one meadow type; rather, they could all be found in a range of meadow types. Thus, one potential explanation for this result might be that the hydric meadows are on the upper end of the moisture conditions butterfly species can stand, even for some of the species having hydric affinities. With a drought, forbs cover could increase in these willow and sedge-predom- inated hydric meadows (e.g., Debinski et al. 2010), providing additional nectar sources for butterflies. Decreasing moisture conditions because of droughts might thus increase habitat suitability of these meadows, particularly for the M2 meadows, leading to an increase in butterfly abundance. Moist meadows could then act as a kind of refuge for species with broader affinities. Such a change could also be envisioned as a sort of shift in habitats along the gradient where each habitat type shifts one step along the gradient to a more dry condition. However, it is important to point out that not all of the species with hydric affinities increased in abundance in these sites. Species that are more restricted in their hydric habitat selection (e.g., Boloria selene and Lycaena hyllus) showed decreasing trends. Finally, we consider species responses in the context of their range distribution for some of the species showing the most dramatic changes. The most abundant xeric species driving increasing trends (e.g., Cercyonis oetus and Lyceana heteronea) were broad-ranging western-U.S. species. It is interesting to note that they also showed some of the strongest increases in the most xeric meadows (Fig. 3). However, there was no PLATE 1. The butterfly Coenonympha haydenii is a species endemic to the Greater Yellowstone Ecosystem, and it showed particularly strong decreases in abundance from 1997 to 2007. Here it is pictured perched on Collomia linearis. Photo credit: J. C. Caruthers. May 2013 1043GRADIENT-BASED VULNERABILITY ANALYSIS consistent pattern observed in the size of geographic range for hydric species showing declines (e.g., Diamond et al. 2011). Both species with narrow ranges (e.g., Coenonympha haydenii ) and broad ranges (e.g., Plebejus saepiolus and Erebia epipsodea) showed declines. It may be noteworthy, however, that C. haydenii is endemic to the Greater Yellowstone region, and is the only species with such a narrow range. In summary, both our single-species, and habitat affinity-based trend analyses showed significant changes in the butterfly community over an 11-year time period. During this drought, species adapted to drier conditions predominantly exhibited increasing abundances, while species adapted to hydric conditions predominantly exhibited decreases in abundance. As expected, respons- es to drought differ drastically among the habitat types. More importantly, however, we found that drought might have turned suboptimal habitats into more favorable sites for some of the assumed losers of climate change. Those species with hydric habitat affinity, yet some flexibility in their habitat use, were able to take advantage of hydric meadows as meadows became drier. It is possible that a few years of wetter conditions may reverse some of the trends we have described (e.g., Morecroft et al. 2002). However, regional models of global climate change for the northern Rocky Moun- tains predict warmer temperatures (Reiners et al. 2003), and the western United States has generally been characterized by a hotter and drier climate, with an average of ;0.558C (1.08F) warmer during 2003–2007 as compared to the 20th-century average (Saunders et al. 2008). Thus, future drought conditions could have important consequences for the overall species diversity of the ecosystem. As such, species associated with wetter meadows may be at risk. Because butterflies are well studied, they provide some of the most comprehensive information about invertebrate community responses to climate change. This fine-scale gradient-level response of a diverse, well-studied group of insects may provide a window into understanding how broader biodiversity patterns could change with future climate change. ACKNOWLEDGMENTS Funding for this project was provided by EPA STAR Grant R825155, NSF LTREB Grant DEB 0518150, and Grand Teton Natural History Association. Thanks goes to Yellowstone National Park, Grand Teton National Park, University of Wyoming National Park Service Research Station, and the U.S. Forest Service for providing support and housing. We would like to acknowledge the many research technicians and field assistants, particularly D. Friedrick, D. Goshorn, T. Randall, and G. McAndrews. This paper benefitted from discussions with A. Toth, D. Warner, J. Sherwood, R. Moranz, J. Delaney, C. Ferris, and review by O. Schweiger and one anonymous reviewer. LITERATURE CITED Aldridge, G., D. W. Inouye, J. R. K. Forrest, W. A. Barr, and A. J. Miller-Rushing. 2011. Emergence of a mid-season period of low floral resources in a montane meadow ecosystem associated with climate change. Journal of Ecology 99:905–913. Bird, C. D., G. J. Hilchie, N. G. Kondla, E. M. Pike, and F. A. H. Sperling. 1995. Alberta butterflies. Provincial Museum of Alberta, Edmonton, Canada. Breed, G. A., S. Stichter, and E. E. Crone. 2012. Climate-driven changes in northeastern US butterfly communities. Nature Climate Change. http://dx.doi.org/10.1038/NCLIMATE1663 Debinski, D. M., M. E. Jakubauskas, and K. Kindscher. 2000. Montane meadows as indicators of environmental change. Environmental Monitoring and Assessment 64:213–225. Debinski, D. M., and J. A. Pritchard. 2002. A field guide to butterflies of the Greater Yellowstone Ecosystem. Roberts Rinehart Publishing, Maryland, USA. Debinski, D. M., C. Ray, and E. H. Saveraid. 2001. Species diversity and the scale of the landscape mosaic: do scales of movement and patch size affect diversity? Biological Con- servation 98:179–190. Debinski, D. M., R. E. Van Nimwegen, and M. E. Jakubaus- kas. 2006. Quantifying relationships between bird and butterfly community shifts and environmental change. Ecological Applications 16:380–393. Debinski, D. M., H. Wickham, K. Kindscher, J. C. Caruthers, and M. Germino. 2010. Montane meadow change during drought varies with background hydrologic regime and plant functional group. Ecology 91:1672–1681. Dennis, R. L. H., and T. H. Sparks. 2006. When is a habitat not a habitat? Dramatic resource use changes under differing weather conditions for the butterfly Plebejus argus. Biological Conservation 129:291–301. Diamond, S. E., A. M. Frame, R. A. Martin, and L. B. Buckley. 2011. Species’ traits predict phenological responses to climate change in butterflies. Ecology 92:1005–1012. Ehrlich, P. R., D. D. Murphy, M. C. Singer, C. B. Sherwood, R. R. White, and I. L. Brown. 1980. Extinction, reduction, stability and increase: the responses of checkerspot butterfly (Euphydryas) populations to the California drought. Oeco- logia 46:101–105. Ehrlich, P. R., and P. H. Raven. 1964. Butterflies and plants: a study in coevolution. Evolution 18:586–608. Ferris, C. D., and F. M. N. Brown. 1981. Butterflies of the Rocky Mountain states. University of Oklahoma Press, Norman, Oklahoma, USA. Fleishman, E., J. P. Fay, and D. D. Murphy. 2000. Upsides and downsides: contrasting topographic gradients in species richness and associated scenarios for climate change. Journal of Biogeography 27:1209–1219. Forister, M. L., A. C. McCall, N. J. Sanders, J. A. Fordyce, J. H. Thorne, J. O’Brien, D. P. Waetjen, and A. M. Shapiro. 2010. Compounded effects of climate change and habitat alteration shift patterns of butterfly diversity. Proceedings of the National Academy of Sciences USA 107:2088–2092. Forister, M. L., and A. M. Shapiro. 2003. Climatic trends and advancing spring flight of butterflies in lowland California. Global Change Biology 9:1130–1135. Gilbert, L. E., and J. T. Smiley. 1978. Determination of local diversity in phytophagous insects: host specialists in tropical environments. Symposium of the Royal Entomological Society of London 9:89–104. González-Megı́as, A., R. Menéndez, D. Roy, T. Brereton, and C. D. Thomas. 2008. Changes in the composition of British butterfly assemblages over two decades. Global Change Biology 14:1464–1474. Grace, J. B., S. Harrison, and E. Damschen. 2011. Local richness along gradients in the Siskiyou herb flora: R. H. Whittaker revisited. Ecology 9:108–120. Jongman, R. H. G., C. J. F. ter Braak, and O. F. R. Tongeren. 1995. Data analysis in community and landscape ecology. Cambridge University Press, Cambridge, UK. Kearney, M. R., N. J. Briscoe, D. J. Karoly, W. P. Porter, M. Norgate, and P. Sunnucks. 2010. Early emergence in a DIANE M. DEBINSKI ET AL.1044 Ecology, Vol. 94, No. 5 butterfly causally linked to anthropogenic warming. Biology Letters 6:674–677. Lenoir, J., J. C. Gégout P. A. Marquet, P. de Ruffray, and H. Brisse. 2008. A significant upward shift in plant species optimum elevation during the 20th century. Science 320: 1768–1771. Morecroft, M. D., C. E. Bealey, O. Howells, S. Rennie, and I. P. Woiwod. 2002. Effects of drought on contrasting insect and plant species in the UK in the mid-1990s. Global Ecology and Biogeography 11:7–22. NABA [North America Butterfly Association]. 2001. Checklist of North American butterflies occurring north of Mexico. Second edition. NABA, Morristown, New Jersey, USA. Oksanen, J., R. Kindt, P. Legendre, B. O’Hara L. G. Simpson, and H. H. Stevens. 2008. vegan: community ecology R package version 1.11-0. R Foundation for Statistical Computing, Vienna, Austria. Packer, L. 1994. The extirpation of the Karner blue butterfly in Ontario. Pages 143–151 in D. A. Andow, R. J. Baker, C. P. Lane, editors. Karner blue butterfly: a symbol of a vanishing landscape. Minnesota Agricultural Experiment Station, University of Minnesota, St. Paul, Minnesota, USA. Parmesan, C. 2006. Ecological and evolutionary responses to recent climate change. Annual Review of Ecology, Evolution, and Systematics 37:637–69. Parmesan, C., et al. 1999. Poleward shifts in geographical ranges of butterfly species associated with regional warming. Nature 399:579–583. Parmesan, C., and G. Yohe. 2003. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421:37–42. Pinheiro, J., and D. Bates. 2009. Mixed effects models in S and S-Plus. Springer, New York, New York, USA. Pounds, J. A., et al. 2006. Widespread amphibian extinctions from epidemic disease driven by global warming. Nature 439: 161–167. R Development Core Team. 2010. R: a language and environment for statistical computing, version 2.12.1. R Foundation for Statistical Computing, Vienna, Austria. Reiners, W. A., et al. 2003. Natural ecosystems 1. The Rocky Mountains. Pages 145–184 in F. Wagner, and T. Stohlgren, editors. Preparing for a changing climate: the potential consequences of climate variability and change. Rocky Moun- tain/Great Basin. A report of the RockyMountain/Great Basin Regional Assessment Team for the U.S. Global Change Research Program. Utah State Press, Logan, Utah, USA. Root, T. L. 1988. Environmental factors associated with avian distributional boundaries. Journal of Biogeography 15:489– 505. Root, T. L., J. T. Price, K. R. Hall, S. H. Schneider, C. Rosenzweig, and J. A. Pounds. 2003. Fingerprints of global warming on wild animals and plants. Nature 421:57–60. Saunders, S., C. Montgomery, T. Easley, and T. Spencer. 2008. Hotter and drier: the west’s changed climate. Rocky Mountain Climate Organization and the Natural Resources Defense Council, New York, New York, USA. Saveraid, E. H., D. M. Debinski, K. Kindscher, and M. E. Jakubauskas. 2001. A comparison of satellite data and landscape variables in predicting bird species occurrences in the Greater Yellowstone Ecosystem. Landscape Ecology 16: 71–83. Scott, J. A. 1986. The butterflies of North America. Stanford University Press, Palo Alto, California, USA. Silvertown, J., M. E. Dodd, D. Gowing DJG, and J. O. Mountford. 1999. Hydrologically defined niches reveal a basis for species richness in plant communities. Nature 400: 61–63. Su, J. C., D. M. Debinski, and K. Kindscher. 2004. Beyond species richness: Community similarity as a measure of cross taxon congruence for Coarse-Filter Conservation. Conser- vation biology 18:167–173. Summerville, K. S., and T. O. Crist. 2003. Determinants of lepidopteran community composition and species diversity in eastern deciduous forests: roles of season, eco-region and patch size. Oikos 100:134–148. Walther, G. R., E. Post, P. Convey, A. Menzel, C. Parmesan, T. J. C. Beebee, J.-M. Fromentin, O. Hoegh-Guldberg, and F. Bairlein. 2002. Ecological responses to recent climate change. Nature 416:389–395. Whittaker, R. H. 1956. Vegetation of the Great Smoky Mountains. Ecological Monographs 26:2–80. Whittaker, R. H. 1972. Evolution and measurement of species diversity. Taxon 21:213–251. Wickham, H. 2009. ggplot2: elegant graphics for data analysis. Springer, New York, New York, USA. SUPPLEMENTAL MATERIAL Appendix A Butterfly species categorized by habitat affinity (Ecological Archives E094-092-A1). Appendix B Mixed-effects model incorporating habitat affinity (Ecological Archives E094-092-A2). Appendix C Single-species model summary (Ecological Archives E094-092-A3). May 2013 1045GRADIENT-BASED VULNERABILITY ANALYSIS http://www.esapubs.org/archive/ecol/E094/092/ http://www.esapubs.org/archive/ecol/E094/092/ http://www.esapubs.org/archive/ecol/E094/092/ << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles true /AutoRotatePages /All /Binding /Left /CalGrayProfile (Gray Gamma 2.2) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Error /CompatibilityLevel 1.4 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.0000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams false /MaxSubsetPct 100 /Optimize true /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile (Color Management Off) /AlwaysEmbed [ true ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 300 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth 8 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.00333 /EncodeColorImages true /ColorImageFilter /FlateEncode /AutoFilterColorImages false /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 300 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth 8 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.00333 /EncodeGrayImages true /GrayImageFilter /FlateEncode /AutoFilterGrayImages false /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.00083 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName () /PDFXTrapped /False /CreateJDFFile false /SyntheticBoldness 1.000000 /Description << /CHS /CHT /DAN /DEU /ESP /FRA /ITA /JPN /KOR /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken voor kwaliteitsafdrukken op desktopprinters en proofers. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.) /NOR /PTB /SUO /SVE /ENU (Use these settings to create Adobe PDF documents for quality printing on desktop printers and proofers. Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.) >> /Namespace [ (Adobe) (Common) (1.0) ] /OtherNamespaces [ << /AsReaderSpreads false /CropImagesToFrames true /ErrorControl /WarnAndContinue /FlattenerIgnoreSpreadOverrides false /IncludeGuidesGrids false /IncludeNonPrinting false /IncludeSlug false /Namespace [ (Adobe) (InDesign) (4.0) ] /OmitPlacedBitmaps false /OmitPlacedEPS false /OmitPlacedPDF false /SimulateOverprint /Legacy >> << /AddBleedMarks false /AddColorBars false /AddCropMarks false /AddPageInfo false /AddRegMarks false /ConvertColors /NoConversion /DestinationProfileName () /DestinationProfileSelector /NA /Downsample16BitImages true /FlattenerPreset << /PresetSelector /MediumResolution >> /FormElements false /GenerateStructure true /IncludeBookmarks false /IncludeHyperlinks false /IncludeInteractive false /IncludeLayers false /IncludeProfiles true /MultimediaHandling /UseObjectSettings /Namespace [ (Adobe) (CreativeSuite) (2.0) ] /PDFXOutputIntentProfileSelector /NA /PreserveEditing true /UntaggedCMYKHandling /LeaveUntagged /UntaggedRGBHandling /LeaveUntagged /UseDocumentBleed false >> ] >> setdistillerparams << /HWResolution [1200 1200] /PageSize [612.000 792.000] >> setpagedevice