Impacts of a millennium drought on butterfly faunal dynamics

Background

Climate change is challenging plants and animals not only with increasing temperatures, but also with shortened intervals between extreme weather events. Relatively little is known about diverse assemblages of organisms responding to extreme weather, and even less is known about landscape and life history properties that might mitigate effects of extreme weather. Our aim was to address this knowledge gap using a multi-decadal dataset of 163 butterfly species that recently experienced a millennium-scale drought. To understand faunal dynamics in the context of the millennium drought, we investigated the behavior of phenology (including date of first flight), species richness and diversity indices through time at 10 study sites spanning an elevational gradient. Linear models were developed to understand the differential sensitivity of butterflies to climate at low and high elevations.

Results

Dates of first flight advanced across the elevational gradient during the drought, leading to an overall expansion of the flight window at low elevations and a compression of the flight window in the mountains. The number of species observed per year increased at lower elevations but decreased at higher elevations, apparently as a consequence of extreme sensitivity to hot and dry conditions.

Conclusion

Montane populations may be more sensitive to climatic extremes than expected based on availability of microclimates and spatial heterogeneity, while low-elevation populations (despite existing in degraded habitats) are buffered by life history plasticity.

Extreme weather events are occurring with increasing severity and frequency in recent decades, a trend linked to ongoing anthropogenic climate change and shifts in atmospheric circulation [1, 2]. While many studies of species and ecosystems have looked at climate impacts by using the average change in historical or projected climatic conditions [3, 4], we know far less about the ecological impacts of extreme weather events on wild organisms [5,6,7]. In part this is due to the regional geography in which extreme weather events occur, which requires baseline biological data as well as appropriately-scaled climatological data for a particular region prior to an event [8]. Droughts have been occurring more frequently, and across greater spatial scales, yet there are few studies addressing the effects of extreme drought on insects [9]. In general, warmer and drier conditions are known to have both positive and negative effects on insect populations. Positive effects can be associated with more rapid development and less time spent in vulnerable juvenile stages [10], while negative effects include reduced nectar availability and early senescence of host plants [11]. Here we take advantage of decades of data on 163 butterfly species across an elevational gradient in Northern California (Fig. 1a) to investigate the effects of a recent, extreme drought on butterfly populations. The years 2011 to 2015 were characterized by both record-setting high temperatures in California and low levels of precipitation, which combined to produce water deficits that were without precedent in more than 1000 years of reconstructed climatic records [12].

a Elevational profile of Northern California (left) and map of the same area (on the right) with ten study sites indicated on both; site names as follows, from west to east: Suisun Marsh (SM), Gates Canyon (GC), West Sacramento (WS), North Sacramento (NS), Rancho Cordova (RC), Washington (WA), Lang Crossing (LC), Donner Pass (DP), Castle Peak (CP) and Sierra Valley (SV). b Average dates of first flight and c last flight across species at each location and year. d Average flight days, which are the average fraction of days individual species are observed per year. In panels b, c, and d, color coding for individual lines corresponds to sites as in panel a, and the data are shown as z-standardized values. Grey rectangles in panels b, c, and d, indicate the major drought years from late 2011 into 2015

Full size image

Investigations of butterflies at our focal sites have reported that a majority of populations at the lowest elevations have been in decline since at least the mid-1990s [13], which has been attributed to changes in land use and warming temperatures [14]. Populations at higher elevations, in contrast, have shown relatively less directional change over time, with the exception of a decline in more dispersive, disturbance-associated species that rely on demographic connections with lower-elevation source populations [15]. Previous analyses of abiotic effects have noted responses to weather that were heterogeneous and idiosyncratic among sites and species [16, 17]. For example, associations with climatic variables often differ in sign among congeneric species [18], and even among populations of a single species [17].

The extremity of the 2011–2015 drought provides a unique opportunity to study abiotic effects on butterfly populations that have already been well characterized in terms of climatic relations and long-term demographic trends. Specifically, we asked the following: 1) does an extraordinary, multi-year drought elicit a faunal response that is extreme relative to faunal behavior in previous dry years? 2) Are impacts on the fauna consistent or divergent across elevations? A theoretical expectation is that organisms living in more heterogeneous environments should be more resilient to extremes of temporal variation [19, 20]. We predicted that butterflies at montane sites would be robust relative to populations at lower elevations in landscapes that are both less spatially variable and already impacted by a history of human activity. 3) Finally, we asked if population-level responses to drought are mediated by phenological shifts. Species that are able to begin activity earlier in the spring might reach higher population densities [21, 22], potentially offsetting detrimental drought effects. Another possibility is that ectotherms exposed to longer growing seasons could fall into a developmental trap by which extra generations fail because of insufficient time [23].

Butterfly data

Ten study sites (Fig. 1a) were visited biweekly (every 2 weeks) by a one of us (AMS) for between 45 and 29 years, depending on the site, and only during good “butterfly weather” when conditions were suitable for insect flight (nearly year round at low elevations, and during a more narrow period at higher elevations). At each site, a fixed route was walked and the presence and absence of all butterfly species noted. Maps of survey routes and site-specific details, as well as publically-archived data can be found at http://butterfly.ucdavis.edu/.

For most analyses, we restricted data to a common set of years, from 1988 to 2016, for which we have data from all sites (the plots in Fig. 1 that go back to 1985 are an exception, and do not include all sites in the first few years). Plots and analyses (described below) primarily involve species richness or phenological data, specifically dates of first flight (DFF) or dates of last flight (DLF). The latter two variables (DFF and DLF) involved filtering to avoid biases associated with variation in the intensity or timing of site visits. Specifically, DFF values were only used if they were proceeded by an absence; in other words, there must be a negative observation before a positive observation is taken as a DFF record. Similarly, DLF values were not used if they were not followed by an absence; so any species observed on the last visit to a site in a particular year did not have a record of DLF for that year. If a species was only observed on a single day in a particular year, then that date was used as a DFF (and only if proceeded by an absence) but not as a DLF, in order to not use the same data point twice.

For a subset of years and sites, absolute counts of individual butterflies (by species) were taken in addition to the presence/absence data; this was done at the 5 lowest elevation sites starting in 1999. These data are used here to investigate the dynamics of the low elevation butterflies during the drought years, specifically the relationship between phenological shifts, changing abundance and dependence on irrigation. For the latter (dependence on irrigation), one of us (AMS) ranked species a priori (without knowing the results of analyses) based on natural history observations. Dependency on irrigated areas was categorized as follows: 1), butterfly species that are essentially independent of irrigation; 2), species that use irrigated, non-native hosts in some areas as well as native, non-irrigated hosts in other areas; 3), species that use irrigated, non-native hosts in at least one of multiple generations; and 4), species that are completely dependent on irrigated, non-native hosts.

Weather variables

Analyses included the following weather variables: maximum and minimum daily temperatures, total precipitation, and a sea surface temperature variable associated with regional conditions [17]. Following previous analyses [15], maximum and minimum temperatures were averaged and precipitation was totaled from the start of September of the previous year through August of the current year. Previous fall through current summer is a useful climatological time period in a Mediterranean climate and captures precipitation and overwintering conditions that potentially affect butterflies through both direct effects on juvenile and adult stages, and indirect effects through host and nectar plants. Data were generated as monthly values using the PRISM system (Parameter-elevation Relationships on Independent Slopes Model, PRISM Climate Group; http://prism.oregonstate.edu) for latitude and longitude coordinates at the center of each survey route.

As a complement to the site-specific, PRISM-derived weather variables, we used an index of sea-surface temperature associated with the El Niño Southern Oscillation (ENSO). Specifically, we used the ONI (Oceanic Niño Index) values for December, January and February (a single value is reported for those winter months; http://cpc.ncep.noaa.gov) from the winter preceding butterfly observations for a given year. Higher values of this index correspond to regionally warmer and wetter conditions. We also downloaded snowfall data from the Central Sierra Snow Lab located near our Donner Pass site (station number 428, https://wcc.sc.egov.usda.gov/reportGenerator/), but preliminary investigations found that annual snowfall totals were highly correlated with annual precipitation totals. Correlation coefficients between snowfall and precipitation were between 0.80 and 0.88, and the inclusion of snowfall caused variance inflation factors from linear models (described below) to often exceed 10; thus snowfall was not included in final models. In contrast, correlations among other weather variables (maximum temperatures, minimum temperatures, precipitation, and ENSO values) were generally lower: across sites and weather variables, the mean of the absolute value of correlation coefficients was 0.31 (standard deviation = 0.23).

Weather variables that were included in models were z-standardized within sites to be in units of standard deviations. This allows variables from sites with different average conditions (e.g., mountain and valley sites) to be readily compared and, more important, it allows for slopes from multiple regression models to be compared among weather variables that are themselves measured on different scales (as is the case with precipitation and temperature).

Analyses

Statistical investigations involved two phases. First, we used plots of z-standardized data to visualize patterns in phenology (DFF and DLF) and flight days over time; the latter variable, the number of days flying, was expressed as the fraction of days that a species is observed divided by the number of visits to a site per year (this has been referred to as “day positives” in other publications using these data [17]). DFF, DLF, and flight days were z-standardized within species at individual sites and then averaged across species to facilitate comparisons of patterns across sites and years. We also used plots of species richness to explore change over time at each site. Plots of richness were accompanied by splines (with 5 degrees of freedom) and predicted values from random forest analyses [24], which both allow for visualization of non-linear relationships. The spline analysis has the advantage of producing smoothed relationships (between richness and years), while the random forest analysis, performed with the randomForest [25] package in R [26], has the complementary advantage of being able to incorporate covariates (in this case the number of visits per year) as well as the advantage of making no assumptions about the shape of the relationship (between richness and years, while controlling for sampling effort). Patterns in species richness over time were also explored using diversity indices and Hill numbers that weight rare and common species differently (at different levels of q, which determines the sensitivity of the analysis to rare species) [27,28,29], using the vegetarian (v1.2) package [30] in R [26]. In addition, we used the combination of spline and random forest analyses to investigate changes in abundance (numbers of individuals observed per year) at the low elevation sites where abundance data were available. As with other variables, abundance values were z-standardized within species and sites, and z-scores were averaged across species.

Following the visualization phase of investigation, we developed simple linear models that were focused on prediction of dates of first flight (DFF) and species richness, summarized as described above (as z-scores averaged across species within each site and year). Independent variables for both sets of models included average daily minimum temperatures, average daily maximum temperatures, total precipitation, ENSO (ONI sea surface temperatures), sampling effort (number of visits), and year. Models of species richness included date of first flight as an additional variable because we were interested in the possibility that the timing of species emergence affects butterfly populations and consequently observed species richness. These models (for both DFF and richness) were estimated for each site individually and for high and low sites as groups of sites (5 sites in each model). Additional model complexities were explored that included interactions between weather variables, time lags (effects of previous years on current year dynamics) and cumulative effects (sliding windows of averaged precipitation values). Interactions were rare, but we report interactions between weather variables that were significant at P < 0.05. Lagged and cumulative weather variables did not add to the explanatory power of models and the individual lagged and cumulative effects were rarely significant and not discussed further.

As we have found elsewhere for analyses of phenology and richness with these data [13, 31], linear models satisfied assumptions of normality and homogeneity of variance, and autocorrelation plots were examined to confirm independence of residuals. To address potential collinearity among predictor variables, variance inflation factors were investigated and found generally to be between 0 and 5, and in a few cases between 5 and 10. For instances where inflation factors approached 10, quality control was conducted by including and excluding correlated variables to verify that estimated β coefficients were not affected. Linear models were also used to test the hypothesis that phenological shifts at low elevations have demographic consequences for individual species. For each species at the lowest sites (SM, WS, NS, and RC; see Fig. 1 legend for site abbreviations), we separately regressed dates of first flight against years, and annual abundance against years. Slopes from those regressions were then compared using correlation to ask if species that shifted to an earlier flight (negative slopes for DFF versus years) were also species that became more abundant (positive slopes of abundance versus years). This was done for the years 2008–2016 to capture the transition into the millennium drought years, and only included species that were present in at least 6 of those years. As with other analyses, linear models were performed and assumptions investigated using R [26].

A prolonged and consistent shift towards earlier spring flight during the recent drought years can be seen in Fig. 1b. While the shift in phenology is evident across elevations, the dynamics of the flight window diverge later in the season: at higher elevations, the date of last flight shifted to an earlier time during the drought, while at lower elevations the last flight dates from 2011 to 2015 are closer to the long-term average (Fig. 1c; also see Additional file 1: Figure S1 for the same results without sampling filters as described above). In other words, the total flight window expanded at lower elevations, while in the mountains the flight window shifted and compressed towards the start of the season, a change that is reflected in fewer overall flight days at higher sites (Fig. 1d).

Along with the recent reduction in the average number of days that butterflies were observed flying at higher elevations during the drought years, there have been fewer butterfly species observed per year at the same sites (Fig. 2a-e). In some cases, the millennium

Species richness (a through j) for all sites, and abundance (k through o) for low elevation sites (abundance data were only collected at the low sites). Two letter site names and colors follow Fig. 1a. In all panels, patterns are visualized with both a spline fit with five degrees of freedom (thin black line) and predicted values from random forest analysis (thick gray line) incorporating variation in sampling effort. In panels k through o, values plotted are z-standardized values of total abundance (number of individuals) per year averaged across species

Full size image

drought was associated with a discrete downturn (e.g., Fig. 2b and c), while at other montane sites the recent drought years contributed to ongoing, negative trends (Fig. 2d and e). A downward trend in species richness is less evident at the highest site (CP, Fig. 2a), which previous analyses have found to be receiving immigrants from lower elevations as populations shift upslope in warmer years [15]. In a dramatic reversal of fortunes, the lowest elevation sites during the millennium drought experienced some of their most productive years in nearly two decades, as reflected both in numbers of species (Fig. 2g-j) and numbers of individuals (Fig. 2l-o). Results shown in Fig. 2 are for simple richness (the number of species observed per year). We repeated analyses using alpha diversity Hill numbers that down weight the importance of rare species (Additional file 1: Figures S2, S3), and found similar results for all sites except for GC, where a long-term decline in the number of species becomes evident when rare, or transient species have less influence.

Why did the low elevation sites apparently rebound during the drought years? Using the lowest sites (RC, WS, NS and SM) and a span of years starting just before the millennium drought, we discovered a potential effect of phenological plasticity. Specifically, species whose first flight shifted to an earlier day were the species that became more abundant (r = − 0.50, P < 0.001; Fig. 3). Butterflies at the lowest elevations are almost entirely multivoltine, and an earlier start for those species could have allowed populations more time to build to greater densities which might have in turn translated to overall higher levels of observed abundance and species richness (Fig. 2). However, we note that causality in the other direction is also possible (higher abundance in favorable years leading to earlier detection of first flight dates), and the two possibilities cannot be separated at this time. Nevertheless, the association between change in abundance and change in first flight date (Fig. 3) raises the potential importance of phenological flexibility, which we sought to further understand by modeling the dates of first flight as a function of maximum and minimum temperatures, precipitation, and El Niño (ENSO) conditions. Models explained 60% of the inter-annual variation at low elevations (F _6,138  = 35.17, P < 0.001), and 72% at the higher elevation sites (F _6,138  = 59.68, P < 0.001) (Additional file 1: Table S1). Minimum and maximum temperatures had negative effects on first flight dates (warmer temperatures lead to earlier flights), and the effect of the former was most noticeable at higher elevations (Fig. 4a and b). Precipitation, as reflected by local weather and ENSO conditions, had a delaying effect on phenology (positive β coefficients in Fig. 4c and d), which is expected as wetter conditions are associated with cooler, cloudy days and delayed spring emergence.

Change in abundance as a function of change in date of first flight: species that shifted to an earlier date of first flight (negative values on the x axis) during the years 2008–2016 were species that increased their overall abundance (positive values on the y axis) (Pearson correlation = − 0.50, P < 0.001). Symbols represent dependency on irrigated areas, as follows: circles, essentially independent of irrigation; squares, use irrigated, non-native hosts in some areas as well as native, non-irrigated hosts; triangles, use irrigated, non-native hosts in at least one of multiple generations; diamond, completely dependent on irrigated, non-native hosts (there is only one species coded with a diamond: Agraulis vanillae, not native to the area). Categories of irrigation dependence are not significant predictors of changing abundance, either considered in a model with change in date of first flight or alone (F_3,39 = 1.60, P = 0.21)

Full size image

Results from a model of phenology (a through d; DFF = date of first flight) and a model of species richness (e through h; SR = species richness), as well as plots of weather variables through time (i through l; ENSO = El Niño Southern Oscillation). In the model results (a through h), the values shown are β coefficients (with standard errors) that summarize the effect of a particular weather variable (while controlling for other variables) on either phenology (DFF) or species richnes (SR). Trend lines are only shown in plots where the effect of a particular weather variable changes with elevation (* P < 0.05; ** P < 0.01; *** P < 0.001; see Additional file 1: Table S4 for additional details). In panels i through k, weather patterns are visualized using splines with five degrees of freedom; panel l is the El Niño index (ONI) for each year. Weather variables are shown as z-standardized values in panels i through k, and site specific colors in all plots are the same as in Fig. 1a

Full size image

Models of species richness revealed even more pronounced variation in weather effects across elevations, including an increased importance of minimum temperatures (Fig. 4e), maximum temperatures (Fig. 4f), and precipitation at higher elevations (Fig. 4g) (see Additional file 1: Table S2 for full details). The highest elevations are most negatively affected by dry years with warmer nighttime temperatures. While daily maximum temperatures have risen everywhere (Fig. 4j) and patterns of precipitation have fluctuated in concert across sites (Fig. 4k), minimum temperatures have risen most steeply at the mountain sites (Fig. 4i). Models of species richness also included phenology (date of first flight) as an explanatory variable, and we found an overall negative association (Fig. 4h), such that earlier emergence is associated with elevated richness (consistent with the effect of phenology on abundance at low elevations; Fig. 3). However, the beneficial effects of earlier emergence at higher elevations might not be as consequential because of a lack of multivoltine species, or they may simply be outweighed by the negative effects of minimum temperatures. Negative effects of minimum temperatures at the higher elevation sites range from 0.48 fewer species seen for every degree Celsius of warming at WA, to 6.46 fewer species seen for every degree at DP (Additional file 1: Table S3). Effects of ENSO on species richness did not differ across elevations (not shown in Fig. 4, see Additional file 1: Table S4), and were generally weak at individual sites (Additional file 1: Table S2).

In contrast to both the previously documented heterogeneity in population response at our focal sites [16] and the observation from other systems that extreme weather elicits species-specific responses [32], the millennium drought in California produced a consistent response across many sites that included advanced dates of first flight with elevation-specific changes in flight windows and species richness. The resilience exhibited by the lowest elevations appears to be associated with phenological flexibility combined with multivoltine life histories and climatic associations that are less detrimental (in the context of current climatic trends) than biotic-abiotic associations at higher elevations. These results are interesting in the light of the finding that a major drought in 1995 in Great Britain shifted species assemblages towards more wide ranging, multivoltine butterfly species with generalist host associations [33]. The drought response in the British butterflies also included an increase in abundance for a subset of sites, although (in contrast to our findings), the populations that increased during the drought tended to be associated with wetter sites, while we saw positive responses at the hotter and drier low elevations.

Many researchers have hypothesized an impending mismatch between trophic levels as a result of climate change [34]. The results from low elevation butterflies in California perhaps challenge that hypothesis, or at least suggest that a shift in phenology at the consumer trophic level need not always have negative consequences. In addition to having multiple generations each year, populations at the lowest elevations have access to agricultural lands. Although association with irrigation does not appear to predict population dynamics during the drought (Fig. 3), we cannot rule out the possibility that low elevation populations were buffered during the drought by irrigated crops or agricultural margins. If true, a reliance on agriculture would be interesting in the light of a recent hypothesis that long-term declines in low elevations butterfly populations are associated with intensification of pesticide applications [35]. It is possible that the rebound of the drought years could be followed by a more severe decline following concentrated agriculture dependency and toxin exposure.

It has been known for some time that high latitude environments are warming faster than the rest of the planet, with negative consequences for many high latitude species [36] but positive or neutral effects for others [37]. It is only recently that climatologists have become aware that higher elevations may also be experiencing a disproportionate share of warming [38], which raises the issue of how cold-adapted, montane ecosystems will respond. Contrary to the expectation that mountains offer microclimatic refugia and preadapt species for climatic variation [39], we found high elevation butterfly communities to be declining and especially sensitive to dry years with warmer minimum temperatures. Warmer and drier years are associated with lower productivity of mesic-adapted plant communities [40], and shorter windows during which montane plants are optimal for nectar and herbivory [41]. We did not model snowfall because it is highly correlated with annual precipitation at our sites (see Methods), but reduced snowpack in dry, warm years would have additional negative effects including higher overwinter mortality for life history stages that typically spend the winter under a blanket of snow [42]. The daily temperature range (the difference between maximum and minimum temperatures) has been shrinking around the globe [43], but the ecological consequences of this thermal homogenization are poorly understood and not yet incorporated into theoretical expectations of global change biology [44]. The results reported here suggest that we have much yet to learn about organismal responses to extreme weather, and the extent to which different habitat types might or might not buffer populations against climate change [45, 46]. However, there is hope for progress because powerful and simple statistical models predicting faunal dynamics with annual climatic values are possible for ectotherms even in the face of unprecedented climatic variation.

Thanks to the Plant-Insect group at UNR for discussion and ideas that led to the analyses reported here.

Funding

Research was funded by the National Science Foundation (DEB-1638773 to CCN, DEB-1638922 to JAF, and DEB-1638793 to MLF), and MLF was supported by a Trevor James McMinn professorship.

Availability of data and materials

Data analyzed in this paper is available at http://butterfly.ucdavis.edu/.

Affiliations

1. Department of Biology, Program in Ecology, Evolution and Conservation Biology, University of Nevada, Reno, NV, 89557, USA
   • Matthew L. Forister
2. Department of Ecology and Evolutionary Biology, University of Tennessee, Knoxville, TN, 37996, USA
   • James A. Fordyce
3. Department of Biology, Population and Conservation Biology Program, Texas State University, San Marcos, TX, 78666, USA
   • Chris C. Nice
4. Department of Environmental Science and Policy, University of California, Davis, CA, 95616, USA
   • James H. Thorne
   •  & David P. Waetjen
5. Center for Population Biology, University of California, Davis, CA, 95616, USA
   • Arthur M. Shapiro

Contributions

AMS designed and carried out data collection. JHT and DPW managed data entry and curation. MLF, JAF and CCN analyzed the data. All authors contributed to writing. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Matthew L. Forister.

Ethics approval and consent to participate

Not applicable.

Competing interests

The authors declare they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Reprints and Permissions