III. Main text

Introduction

Terrestrial insects, and notably butterflies, are declining in abundance regionally and globally by 1-3% per year (Dirzo et al. 2014; Forister et al. 2021; Van Klink et al. 2020; Wepprichet al. 2019a). Habitat loss, pesticides, and climate change are implicated as interacting contributors to declines (Harvey et al.2023; Wagner et al. 2021; Warren et al. 2021). This complexity of stressors make predicting population trends and climate responses for any butterfly species particularly uncertain, as some species fare worse and some better under new climatic conditions (Breedet al. 2013; Crossley et al. 2021). Gradual warming and extreme events could affect butterfly abundances through different mechanisms, for example by shifting distributions, exceeding thermal limits of individuals, reducing host plant quality in drought (Harveyet al. 2023), or, the subject of this study, changing phenology and the number of generations they complete annually.
Even though insects have responded to climate change by shifting phenology and ranges to track temperature (Parmesan 2006), a key concern for conservation in a warming climate is whether insects will keep up with the rate of change given the constraints from local adaptations to the past (Forrest 2016). Voltinism, or the number of generations completed in a year, is one insect trait that is highly dependent on climate and a target for natural selection to keep in sync with changes in season length (Cohen 1970; Tauber & Tauber 1976). If season length changes too rapidly as a result of anthropogenic global warming, the locally-adapted environmental cues regulating voltinism may become maladaptive in an evolutionary trap (Schlaepfer et al. 2002). An evolutionary trap has been proposed as a potential cause for a local population collapse in a multivoltine butterfly attempting doomed late season generations, or “lost generations” (Van Dyck et al.2015). It is not known if such traps are widespread in multivoltine insects confronting climate warming and whether lost generations impact population viability. In this paper, we use 27 years of monitoring to test whether 30 multivoltine butterfly species are attempting lost generations, and whether the voltinism changes correlate with long-term population trends.
With thermal constraints on development, insect life history strategies vary spatially across climatic gradients and temporally with annual environmental variability to match the time available for seasonal activity (Danilevskii 1965; Masaki 1961; Tauber & Tauber 1976). In temperate climates, insects survive harsh overwinter conditions through diapause, a developmental pathway that provides a more hardy metabolic state (Koštál 2006). Multivoltine insects use environmental cues to induce pathways of attempting another generation with direct development or pausing normal development to prepare for diapause. Photoperiod is the predominant cue for initiating diapause development in temperate insects (Bradshaw & Holzapfel 2007; Tauber & Tauber 1976). Local adaptations in the photoperiodic response have evolved so that voltinism varies across species’ ranges corresponding to clines in season length and latitudinal gradients of day length (Masaki 1961). As photoperiod at any location is a consistent cue, it is a prime candidate for becoming less informative and mismatched with rapid environment change (McNamara et al. 2011), especially for insects where climate warming shifts the phenology of life stages sensitive to cues. However, photoperiodism is a labile trait, with rapid adaptive evolution in range-expanding insects (Bean et al. 2012; Bradshaw & Holzapfel 2001; Ittonen et al. 2022; Urbanski et al.2012). Open questions are whether locally-adapted multivoltine insects will be able to simultaneously adjust responses to diapause induction cues, and hence voltinism, through adaptation or phenotypic plasticity, and whether voltinism mismatches contribute to negative demographic consequences and, in turn, widespread insect declines (Forrest 2016).
Despite well-documented increases in annual number of generations in response to earlier emergence in animals, the effect of changes in voltinism on population dynamics is rarely tested (Altermatt 2010; Ileret al. 2021; Pöyry et al. 2011). Our approach quantifies both changes in voltinism and population demography to explore the potential consequences for matches and mismatches (Fig. 1). Given the plasticity and evolutionary potential in this trait, our expectation might be that it is largely adaptive for insects to increase voltinism in response to warming. Warmer autumn temperatures increase the costs of prolonged diapause, which is mitigated with direct development into extra generations (Gomi 2000; Nielsen et al.2022). In butterflies, multivoltine species have had more positive, long-term population trends than inflexible univoltine species (Colomet al. 2022; Macgregor et al. 2019; Michielini et al. 2021; Wepprich et al. 2019a). The potential for lost generations, given this context, may be limited to specific circumstances without widespread population impacts, akin to phenological mismatches in resource-consumer systems that have few long-term consequences for fitness or demography (Kharouba & Wolkovich 2023). However, mismatched voltinism is well known in insect introductions and in fact could create an alternative developmental trap – “missing generations” – if the photoperiod response causes diapause too early (Grevstad et al. 2022; Grevstad & Coop 2015). Mismatches, or either lost or missing generations, may be counterbalanced by the rewards of increased voltinism (i.e., “bonanzas” in Kerr et al. (2020)). We expect voltinism to affect population growth rates depending on winter conditions and that these matched or mismatched demographic responses will impact long-term population trends (Fig. 1).
Using nearly three decades of butterfly monitoring, we take a critical look at the lost generation hypothesis to test it more broadly. First, we determine how voltinism has changed and how responsive it is to phenological shifts across sites and years. Second, we pair observed voltinism differences with winter onset and severity to determine whether overwinter population growth rates show evidence of matches or mismatches. Third, we correlate species’ voltinism responses with long-term trends to determine whether matches or mismatches affect population viability. We find that increasing voltinism is the strongest predictor of higher long-term population trends across species. Even in species with potential voltinism mismatches in some circumstances, generally we find that overwinter demography would benefit from increased voltinism.

Material and methods

Study system

Butterfly communities have been monitored from 1996 to present across the Midwestern state of Ohio, USA (38.4 – 42.3°N, -84.8 – -80.5°W, 116,100 km2 land area). Fixed transects are walked weekly from April through October using the Pollard method where butterflies seen within a 5-meter box centered on the observer are counted (Pollard & Yates 1993). We used 32,794 transects from 1996-2022 from 161 sites (median 23 surveys per year and 8 years of monitoring) including 21 sites contributing 20 or more years. The monitoring program uses species naming conventions from Iftner et al. (1992) and Opler & Krizek (1984).
Previous analysis showed overall abundance declines of 2% per year across all species. Univoltine species had lower abundance trends than bivoltine and multivoltine species (Wepprich et al. 2019a). For this analysis, we selected non-migratory species with two or more generations with limited overlap and an adequate sample size for the models. The initial set of 42 species included species with partial second generations observed at multiple sites even though they are often classified as univoltine (Iftner et al. 1992; Wepprich et al. 2019a).

Statistical analysis and reproducibility

All analyses were performed in R 4.3.1 (R Core Team 2023; Wepprichet al. 2024). Software packages were used for temperature data (‘daymetr’), phenology models (‘mgcv’), mixture models (‘Mclust’), population models (‘lme4’, ‘lmerTest’, ‘merTools’), phylogenetic models (‘ape’, ‘caper’), and data visualization (‘corrplot’, ‘ggeffects’, ‘tidyverse) (Bates et al. 2015; Hufkens et al. 2018; Knowles & Frederick 2023; Kuznetsova et al. 2017; Lüdecke 2018; Orme et al. 2023; Paradis & Schliep 2019; Scrucca et al.2016; Wei & Simko 2021; Wickham et al. 2019; Wood 2017).

Environmental covariates

Environmental covariates were obtained for each site from Daymet daily temperature estimates from 1996-2022 (Thornton et al. 2022). Degree-days were calculated with the single-sine method (Buckleyet al. 2023). While 10℃ base thresholds are commonly used for temperate insects when physiological constraints are unknown (Caytonet al. 2015), we used 5℃ base thresholds so that early spring observations were not assigned zero accumulated degree-days.
For measuring winter severity, we used the ordinal date of the first hard frost in the fall (daily minimum < -2℃) as winter onset and the 5-month mean of the daily minimum temperature starting October 31st, the end of the monitoring season, aswinter temperature . Other potential variables, such as snow accumulation or autumn degree-days, were not used because they were highly correlated with winter onset and temperature. We calculated environmental trends with linear mixed models with a year covariate and site random intercepts.

Abundance and phenology of separate generations

We modified existing methods for estimating abundance and phenology to smooth weekly counts, interpolate missing surveys, estimate annual abundance indices, and separate generations (Dennis et al. 2016, 2022). The key difference in our approach was using degree-days (i.e., physiological time) as the seasonal timescale, as it standardizes phenological responses across spatial and temporal temperature variation (Barker et al. 2020).
For each species, weekly counts were smoothed with generalized additive models to predict phenological variation caused by temperature variation and latitudinal gradients in voltinism (Hodgson et al. 2011; Wood 2017). We used these models to select multivoltine species for analysis and impute missing weekly counts (Dennis et al. 2013). We modeled weekly counts by degree-days and geographic coordinates to account for voltinism variation and spatial correlation. Random effects accounted for differences in relative abundance at sites and years (Wepprichet al. 2019a).
Mixture models can classify observations in weekly surveys to different generations and estimate generation size and phenology across sites (Dennis et al. 2022; Matechou et al. 2014). For each species, we fit univariate Gaussian mixture models to all counts with both ordinal date and degree-day timescales, selecting the timescale that minimized classification uncertainty. Model fitting required adjustments, such as limiting the number of clusters and requiring equal variance (Matechou et al. 2014). Species with too much generation overlap were removed from analysis.
Weekly counts for each species, site, and year were multiplied by their classification probability for each generation from the mixture model, which permitted a single survey’s counts to be apportioned to different generations. For each generation’s apportioned counts, we calculated a population index with the trapezoidal rule and the peak phenology with a weighted mean.

Last generation size variability

Voltinism was quantified by the last generation size accounting for the size of the penultimate generation. Although the diapause switch for an individual insect is binary, at a population-level both the diapause proportion over time and the last generation size vary as continuous traits (Gomi 2000; Kerr et al. 2020). Last generation size was calculated as the natural logarithm of population growth rates between the penultimate and last generations. We modeled last generation size in response to penultimate generation phenology, density-dependence, and year with variables scaled as follows.
Phenology of the penultimate generation was the predicted driver of last generation size. We used the ordinal date of penultimate generation peak as a proxy for photoperiod exposure. Phenology varies across the climatic gradient of sites and in response to annual weather. We split the phenology variable to account for this hierarchy, including bothsite mean phenology and annual variation in phenologyscaled by site, using within-group centering (van de Pol & Wright 2009). An interaction between these two variables allowed models to estimate site-dependent responses to annual phenological variation. A comparison of these two variables quantified how voltinism responses were driven by local adaptation versus phenotypic plasticity. The degree of local adaptation was estimated as the sensitivity of the last generation size to spatial differences in phenology, while the degree of phenotypic plasticity was estimated as the sensitivity to annual variation in phenology (Phillimore et al. 2012; Roy et al.2015).
Density-dependence was the annual variation in the natural logarithm of the penultimate generation size mean-centered by site. Spurious density-dependence may be detected in population time-series due to the uncertainty in counts (Freckleton et al. 2006), so we assume that its effect includes both observation bias and real density-dependence (Wenger et al. 2022).
Year , scaled by subtracting the monitoring year from the midpoint, was used to estimate the trend in last generation size beyond voltinism changes accounted for by spatial and temporal variation in penultimate generation phenology.
We fit linear mixed effects models for each species with sites as random intercepts. All covariates described above were used as fixed effects. The random effect was removed via likelihood ratio tests if not significant (Kuznetsova et al. 2017).
A community model using all species was fit with the variables in the species model, but with an added variable for the species’ mean penultimate generation peak . We assumed that there would be fundamental differences between earlier- and later-flying species. This species trait interacted with site mean and annual variation in last generation phenology. Species and sites were random intercepts. Model performance was assessed by the variation explained when compared to a simpler model containing only density-dependence, year, and random intercepts. Fixed and random effects were removed from the full model via likelihood ratio tests.

Overwinter population growth rates

If voltinism mismatches were a concern, we would expect an interaction between last generation size and winter conditions that lowered population growth rates (Fig. 1). We estimated overwinter population growth rates as natural logarithms of population growth rates between the penultimate and first springtime generations. We modeled how last generation size, winter onset, winter severity, year, and density-dependence predicted overwinter population growth rates. We describe variable selection and scaling below.
Last generation size varied by species, site, and year. Context-dependent differences in annual voltinism were the main interest in this study. Mean last generation size for each species was used as a voltinism trait in a later analysis. Site mean last generation size and annual variation of last generation size , mean-centered within site, were included in the overwinter model.
Winter onset and winter temperature were mean-centered within sites to create site mean and annual variation terms. Including both mean and annual terms permits an analysis of context-dependent demographic responses to annual environmental conditions (Amburgeyet al. 2018). However, the voltinism gradient follows the climatic gradient across sites, so site mean winter variables were highly correlated with the site mean last generation size. We only included uncorrelated annual variation in winter onset and temperature.
For each species, we fit linear mixed effects models of overwinter population growth rates with site mean last generation size, annual variation in last generation size, annual variation in winter onset, and annual variation in winter temperature and their 3-way interactions. We included density-dependence and year terms without interactions as fixed terms and sites as random intercepts. Variables were removed with likelihood ratio tests.
A community model using all species was fit with all variables in the species model, but with an added variable for the species’ mean penultimate generation peak and its interactions with voltinism and weather terms. Species and sites were included as random intercepts. We fit models as stated above.
We simulated the potential population consequences of a larger last generation, incorporating the variable interactions for each species. It was unclear from the regression models alone if higher population growth in years with well-matched voltinism would outweigh mismatched lost or missing generations. We estimated overwinter population growth rates for each observation with one standard deviation added or subtracted from the last generation size for 100 simulations that accounted for uncertainty in model residuals, fixed effects, and random effects (Knowles & Frederick 2023). We calculated the geometric mean of predicted overwinter population growth rates for each site across years, and then averaged them for a statewide simulated effect of larger last generations.

Voltinism traits and long-term population trends

We correlated species-level traits, derived from the above analysis, with species’ population trends. Statewide trends were calculated with generalized linear mixed models following Wepprich et al. (2019), with 6 additional years of monitoring in this study compared to the previous study. To avoid confounding population trends with voltinism changes, we only used first generation counts. Trends were correlated with the following traits: species phenology (mean peak date of penultimate generation), species mean last generation size, trend in last generation size, sensitivity of last generation size to spatial variation in phenology (i.e., local adaptation), sensitivity of last generation size to annual variation in phenology (i.e., phenotypical plasticity), and the simulated effect of larger last generations on overwinter population growth.
Population trends and traits were estimated with uncertainty from their respective models. To account for this uncertainty, we used a parametric bootstrap of 1000 draws from normal distributions with their estimates and standard errors to perform the correlation analysis and report the median and 95% range of bootstrap results. To account for non-independence of traits due to phylogenetic relatedness between species, we constructed phylogeny from publicly available sequences (Wepprich et al. 2019a, b). We fit models to test for effects of traits on population trends with phylogenetic generalized least squares (Orme et al. 2023) and tested the robustness of the results by refitting models using the parametric bootstrap values.
Following the conceptual model of voltinism matches (Fig. 1), the interaction of two traits—trend in last generation size and simulated effect of a larger last generation on overwinter growth rates—would test for the lost generation hypothesis most directly. We tested if these two traits interacted using the phylogenetic models. To highlight species that fit scenarios proposed in Fig. 1, we visualized these two traits and their population trends.

Results

Environmental covariates

The annual thermal window for insect development, measured by accumulated degree-days, ranges from 2,157 to 3,687 with the combined effects of spatial and temporal variation. Over the study period (1996-2022), monitoring sites warmed and growing seasons lengthened. Annual degree-days increased at a rate of 98.5 degree-days per decade (SE = 2.63, t = 37.5, P < 0.001). First fall frosts trended 1.26 days later per decade (SE = 0.165, t = 7.60,P < 0.001). Mean daily fall (September, October, November) and winter (December, January, February) temperatures increased at 0.478℃ per decade (SE = 0.0168, t = 28.4, P< 0.001) and 0.202℃ per decade (SE = 0.0371, t = 5.45,P < 0.001), respectively.

Abundance and phenology of separate generations

After filtering out species poor discrimination of separate generations, we included thirty species out of forty-two candidates. Generalized additive models demonstrated spatial and temporal variation in voltinism, with the proportion of null deviance explained by species models ranging from 0.56 to 0.94 (median 0.75). Mixture models for 28 of 30 species had lower uncertainty assigning observations to generations when performed on a degree-day timescale rather than an ordinal day timescale. An example species, Phyciodes tharos , has voltinism variation in warmer and cooler sites on degree-day and calendar date timescales (Fig. 2 A-C).
By species, the proportion of counts classified as the last generation ranged from 4% to 64% of the total counts across generations. Species had two to five maximum generations per year. The mean penultimate generation peak phenology ranged from May 21st to September 5th, with earlier phenology generally associated with larger last generations (Fig S1).

Last generation size variability

In the community model across 30 species, last generation size increased over time (β = 0.014, SE = 0.002, P < 0.001, Table S1). Penultimate generation phenology affects species more strongly later in the season. For mid- and late-season species, years with earlier phenology leads to larger last generation sizes and sites with earlier phenology (i.e., warmer sites) also have larger last generation sizes (Figure S2, Table S1). Last generation size varies more strongly across sites (β = -1.546, SE = 0.103, P < 0.001) than across years (β = -0.259, SE = 0.035, P < 0.001), consistent with stronger effects of local adaptation than phenotypic plasticity on voltinism (Table S1).
From species models, 12 species had increasing last generation size over time, while one had decreasing voltinism. In the test for local adaptation, 14 species had smaller last generation sizes when site mean phenology was later, while two species had the opposite response. In the test for phenotypic plasticity, six species had larger last generations with earlier annual phenology and five species had the opposite response (Fig. 3). Of the seven species with sensitivity to both site and annual variation in penultimate generation phenology, five had co-gradient negative responses (including our example species in Fig. 1D), two had co-gradient positive responses, and one species had counter-gradient responses (Fig. 3). All species had strong negative density dependence, where larger penultimate generations lowered the population growth rate producing the last generation.

Overwinter population growth rates

In the community model across 30 species, overwinter population growth rates have declined over time (β = -0.011, SE = 0.002, P < 0.001, Table S2). Larger last generations at warmer sites increased overwinter population growth rates regardless of winter temperature or species phenology. However, overwinter population growth declined with larger last generations at cooler sites with cold winters or early frosts (Fig. 4, Table S2). Winter onset and winter temperature had similar effect sizes in their interactions with other variables, but their interaction with each other shows that negative effects of early winter onset (for species with late season phenology only) are counteracted by positive effects of warmer winter temperatures (Fig. S3).
Nearly all species’ overwinter population growth rates were influenced by both last generation size and winter variables, often with interactions. Species models included site mean last generation size for 24 species, annual variation in last generation size for 25 species, annual variation in winter temperature for 23 species, and annual variation in winter onset for 19 species. While warmer winter temperature was evenly split between having a positive or negative effect on growth rates, earlier frosts had a negative effect in nearly all species. All but one species had strong negative density dependence, where larger penultimate generations lowered the overwinter population growth rate producing the springtime generation. Eleven species had declines over time in overwinter growth rates with no species having an increase.
Simulations for individual species show that larger last generations generally benefit geometric mean growth rates, with a median increase in overwinter growth rate of 0.60 (-0.57 to 0.56 range). Fifteen species had positive consequences with larger last generation sizes, while only two species had negative responses to larger last generation sizes (Fig. 3).

Voltinism traits and long-term population trends

Long-term population trends have a median -0.01 slope per year (-0.08 to 0.11 range) with 10 positive trends (3 significant) and 20 negative trends (6 significant) (Fig. 3). These trends are less negative compared to those across all species, including univoltine, and all generations analyzed in Wepprich et al. (2019).
The species traits positively correlated with population trends were trends in last generation size (r = 0.42, CI -0.03 – 0.67), simulated effects of larger last generations (r = 0.29, CI 0.01 – 0.54), and mean size of the last generation (r = 0.26, CI 0.12 – 0.40) (Fig. 5A). Species with generations later in the season were more likely to have last generation sizes that were sensitive to site phenology (r = -0.42, CI -0.52 – -0.24) and annual variation in phenology (r = -0.63, CI -0.73 – -0.52) of the penultimate generation, matching the community model of last generation size variability (Table S1). These two traits showing the last generation sensitivity to spatial and temporal phenological variation (i.e., local adaptation and phenotypic plasticity, respectively) were positively correlated with each other (r = 0.57, CI 0.27 – 0.76) (Fig. 5A).
The phylogenetic regression models of traits versus population trends demonstrated that trends in last generation size predicted population trends (β = 0.5021, SE = 0.176, P = 0.0083, Table S3, Fig. S4). The interaction of the trend in last generation size and the simulated consequences of larger last generations was not robust when trait uncertainty was included (Table S3, Fig. S5). However, grouping species based on this interaction visualizes the proposed scenarios of voltinism matches and mismatches even if there is no consistent relationship across all species (Fig. 5B).

Discussion

Using long-term butterfly monitoring data from Ohio, we find that increasing voltinism generally benefits butterflies, thus rejecting the generality of the lost-generation hypothesis. If the lost generation hypothesis were commonplace, years with larger attempted last generations would have fared worse when confronted with early frosts or colder winters and these mismatches would be reflected in population trends. Instead, we found that larger last generations generally increase species overwinter growth rates. We also found that long-term population trends are positively correlated with increasing voltinism over time. These lines of evidence suggest that lost generations are rare, and that increasing voltinism may allow many butterflies to track, and even adapt to warming.
The degree to which increasing voltinism has positive effects on population growth varies across space, time, and species. Even without warming, spatiotemporal variation in phenology and winter conditions can flip voltinism mismatches to demographic bonanzas in different years or regions (Kerr et al. 2020). Here, lost generations where larger last generation size leads to lower overwinter growth rates occur in a limited set of circumstances across species. The species and populations most susceptible to lost generations are ones that inhabit cooler sites within this sample. Even in cooler sites, the consequences of increased voltinism are only negative in cold years when larger last generations are confronted with cold winters or early frosts (Fig. 4, Table S2). If observed trends in season length, first fall frost, and winter temperature continue in Ohio, the benefits of adding generations will likely increase. This is one mechanism that may facilitate expansion at poleward range limits (Ittonen et al. 2022; Macgregor et al. 2019), for which Cyllopsis gemma provides an example as it expands in southern Ohio and has the largest magnitude population and voltinism trends (Fig. 3).
Our classification of species traits finds few likely candidates for lost generation mismatches that have both increasing voltinism over time and a simulated negative impact of larger generation sizes (lower right quadrant, Fig. 5B). There are more species that have not increased voltinism despite overwinter models predicting higher growth rates with larger last generations (e.g., Lycaena phlaeus and Papilio polyxenes , upper left quadrant of Fig. 5B). These species with “missing generations” may have voltinism limited by late-season resource availability or might have an unknown evolutionary trap that limits voltinism flexibility. We acknowledge that our search for lost or missing generations in these species could be hindered by analyzing a snapshot in time without knowing how natural selection has already shaped voltinism patterns.
Our results suggest that multivoltine species have the capacity to adapt to longer seasons even though their voltinism response is locally adapted across sites. Local adaptation in butterfly springtime phenology may constrain responses to climate warming compared to phenotypic plasticity (Roy et al. 2015). Even though we found smaller contributions of phenotypic plasticity in the voltinism response, larger last generations seem to match the annual variation in overwinter conditions. There are other factors that influence voltinism beyond the photoperiod exposure that we assume determines diapause. When diapause and voltinism are examined in more detail for any species, nuances of interacting cues and co-evolved traits modulate development to fit within the growing season. For example, faster development rates that allow for more generations can vary by genotype in sawtooth latitudinal clines (Levy et al. 2015) or plastically with exposure to shorter day lengths near season end (Lindestad et al. 2019). Host plant quality, which we presume also varies with annual weather conditions or butterfly density, changes diapause rates and voltinism (Abarca 2019; Henry et al. 2022). Testing for evolution in critical photoperiods for diapause could determine whether Ohio butterflies are adapting to warming by using multiple cues or evolving photoperiod reaction norms (Nielsen & Kingsolver 2020).
Despite the possible evidence for adaptation to climate change, the benefits of an extra generation are not enough to reverse declines of most species. Consistent with flexible voltinism helping species keep pace with a changing climate, we found that species with increasing voltinism have more positive state-wide population trends (Fig. 5A, Fig. S4). However, this subset of multivoltine species were already known to have higher population trends compared to obligate univoltine species (Wepprich et al. 2019a) and even here there are more declines than bonanzas (Fig. 3). Climate change is only one driver of insect declines and species are forced to respond to warming in a highly fragmented and degraded landscape (Hodgson et al. 2022; Warrenet al. 2001). Pesticide use in Ohio, specifically neonicotinoids, is associated with lower abundance in common butterflies (Van Deynzeet al. 2022). To increase insect adaptability to climate change, conservation actions should prioritize actions that increase adaptive capacity. Large, connected populations can maintain a high level of genetic diversity that allows for in-situ adaptation (Habel & Schmitt 2018). The locally-adapted variability in phenological response speaks to the need to maintain variation in microclimates across, and even within, sites to allow for variation in phenology on which natural selection can act (Suggitt et al. 2018). When considering the entire butterfly community, univoltine species do not have the adaptability of multivoltine species and should be a greater conservation concern than voltinism mismatches.
Long-term monitoring allowed us to characterize late-season generations as potential adaptations to a warming climate and not just random or rare events. Insect populations fluctuate notoriously, making it especially challenging to separate signal from noise without long-term data. Further complicating this are the complex ways we have demonstrated that sites and weather interact to shape local population dynamics. With monitoring across sites spanning a 7℃ mean annual temperature gradient confronting 27 years of winter variation, we have a dataset that includes 12,377 observations of voltinism variation and 8,676 observations of overwinter demography. It is only with this extensive record that we can test, and mostly reject, the lost generation hypothesis. Most monitoring studies are designed for surveillance for conservation, however, these long-term datasets are key to understanding how, when, and where species are adapting to climate change or not. With this long-term, large-scale data, we can conclude that it is other factors besides the lost generation hypothesis that explain butterfly declines in Ohio.

Acknowledgements

We thank the volunteers organized by the Ohio Lepidopterists for their monitoring efforts and Jerome Wiedmann for data curation. TW was supported by a graduate fellowship awarded by the U.S. Geological Survey (USGS) Southeast Climate Adaptation Science Center at North Carolina State University.

References

Abarca, M. (2019). Herbivore seasonality responds to conflicting cues: Untangling the effects of host, temperature, and photoperiod. PLOS ONE , 14, e0222227.
Altermatt, F. (2010). Climatic warming increases voltinism in European butterflies and moths. Proc. R. Soc. B. , 277, 1281–1287.
Amburgey, S.M., Miller, D.A.W., Campbell Grant, E.H., Rittenhouse, T.A.G., Benard, M.F., Richardson, J.L., et al. (2018). Range position and climate sensitivity: The structure of among‐population demographic responses to climatic variation. Glob Change Biol , 24, 439–454.
Barker, B.S., Coop, L., Wepprich, T., Grevstad, F. & Cook, G. (2020). DDRP: Real-time phenology and climatic suitability modeling of invasive insects. PLOS ONE , 15, e0244005.
Bates, D., Mächler, M., Bolker, B. & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. J Stat Softw , 67, 1–48.
Bean, D.W., Dalin, P. & Dudley, T.L. (2012). Evolution of critical day length for diapause induction enables range expansion of Diorhabda carinulata , a biological control agent against tamarisk (Tamarix spp.): Evolution and range expansion of a biocontrol agent. Evol Appl , 5, 511–523.
Bradshaw, W.E. & Holzapfel, C.M. (2001). Genetic shift in photoperiodic response correlated with global warming. Proc Natl Acad Sci USA , 98, 14509–14511.
Bradshaw, W.E. & Holzapfel, C.M. (2007). Evolution of Animal Photoperiodism. Annu Rev Ecol Evol S , 38, 1–25.
Breed, G.A., Stichter, S. & Crone, E.E. (2013). Climate-driven changes in northeastern US butterfly communities. Nature Clim Change , 3, 142–145.
Buckley, L.B., Ortiz, B.A.B., Caruso, I., John, A., Levy, O., Meyer, A.V., et al. (2023). TrenchR: An R package for modular and accessible microclimate and biophysical ecology. PLOS Climate , 2, e0000139.
Cayton, H.L., Haddad, N.M., Gross, K., Diamond, S.E. & Ries, L. (2015). Do growing degree days predict phenology across butterfly species?Ecology , 96, 1473–1479.
Cohen, D. (1970). A Theoretical Model for the Optimal Timing of Diapause. Am Nat , 104, 389–400.
Colom, P., Ninyerola, M., Pons, X., Traveset, A. & Stefanescu, C. (2022). Phenological sensitivity and seasonal variability explain climate-driven trends in Mediterranean butterflies. Proc. R. Soc. B. , 289, 20220251.
Crossley, M.S., Smith, O.M., Berry, L.L., Phillips‐Cosio, R., Glassberg, J., Holman, K.M., et al. (2021). Recent climate change is creating hotspots of butterfly increase and decline across North America. Glob Change Biol , 27, 2702–2714.
Danilevskii, A.S. (1965). Photoperiodism and seasonal development of insects. Oliver & Boyd, Edinburgh and London, UK.
Dennis, E.B., Fagard-Jenkin, C. & Morgan, B.J.T. (2022). rGAI: An R package for fitting the generalized abundance index to seasonal count data. Ecol Evol , 12, e9200.
Dennis, E.B., Freeman, S.N., Brereton, T. & Roy, D.B. (2013). Indexing butterfly abundance whilst accounting for missing counts and variability in seasonal pattern. Methods Ecol Evol , 4, 637–645.
Dennis, E.B., Morgan, B.J.T., Freeman, S.N., Brereton, T.M. & Roy, D.B. (2016). A generalized abundance index for seasonal invertebrates.Biometrics , 72, 1305–1314.
Dirzo, R., Young, H.S., Galetti, M., Ceballos, G., Isaac, N.J.B. & Collen, B. (2014). Defaunation in the Anthropocene. Science , 345, 401–406.
Forister, M.L., Halsch, C.A., Nice, C.C., Fordyce, J.A., Dilts, T.E., Oliver, J.C., et al. (2021). Fewer butterflies seen by community scientists across the warming and drying landscapes of the American West. Science , 371, 1042–1045.
Forrest, J.R. (2016). Complex responses of insect phenology to climate change. Curr Opin Insect Sci , 17, 49–54.
Freckleton, R.P., Watkinson, A.R., Green, R.E. & Sutherland, W.J. (2006). Census error and the detection of density dependence: Census error and density dependence. J Anim Ecol , 75, 837–851.
Gomi, T. (2000). Effects of timing of diapause induction on winter survival and reproductive success in Hyphantria cunea in a transition area of voltinism. Entomol Sci , 3, 433–438.
Grevstad, F.S. & Coop, L.B. (2015). The consequences of photoperiodism for organisms in new climates. Ecol Appl , 25, 1506–1517.
Grevstad, F.S., Wepprich, T., Barker, B., Coop, L.B., Shaw, R. & Bourchier, R.S. (2022). Combining photoperiod and thermal responses to predict phenological mismatch for introduced insects. Ecol Appl , 32, e2557.
Harvey, J.A., Tougeron, K., Gols, R., Heinen, R., Abarca, M., Abram, P.K., et al. (2023). Scientists’ warning on climate change and insects. Ecol Monogr , 93, e1553.
Henry, E.H., Terando, A.J., Morris, W.F., Daniels, J.C. & Haddad, N.M. (2022). Shifting precipitation regimes alter the phenology and population dynamics of low latitude ectotherms. Clim Change Ecol , 3, 100051.
Hodgson, J.A., Randle, Z., Shortall, C.R. & Oliver, T.H. (2022). Where and why are species’ range shifts hampered by unsuitable landscapes?Glob Change Biol , 28, 4765–4774.
Hodgson, J.A., Thomas, C.D., Oliver, T.H., Anderson, B.J., Brereton, T.M. & Crone, E.E. (2011). Predicting insect phenology across space and time. Glob Change Biol , 17, 1289–1300.
Hufkens, K., Basler, D., Milliman, T., Melaas, E.K. & Richardson, A.D. (2018). An integrated phenology modelling framework in R. Methods Ecol Evol , 9, 1–10.
Iftner, D.C., Shuey, J.A. & Calhoun, J.V. (1992). Butterflies and skippers of Ohio . Ohio Biol. Surv. Bull. New Series. College of Biological Sciences, Ohio State University, Columbus, Ohio.
Iler, A.M., CaraDonna, P.J., Forrest, J.R.K. & Post, E. (2021). Demographic Consequences of Phenological Shifts in Response to Climate Change. Annu Rev Ecol Evol S , 52, 221–245.
Ittonen, M., Hagelin, A., Wiklund, C. & Gotthard, K. (2022). Local adaptation to seasonal cues at the fronts of two parallel, climate-induced butterfly range expansions. Ecol Lett , 25, 2022–2033.
Kerr, N.Z., Wepprich, T., Grevstad, F.S., Dopman, E.B., Chew, F.S. & Crone, E.E. (2020). Developmental trap or demographic bonanza? Opposing consequences of earlier phenology in a changing climate for a multivoltine butterfly. Glob Change Biol , 26, 2014–2027.
Kharouba, H.M. & Wolkovich, E.M. (2023). Lack of evidence for the match-mismatch hypothesis across terrestrial trophic interactions.Ecol Lett , 26, 955–964.
Knowles, J.E. & Frederick, C. (2023). merTools: Tools for Analyzing Mixed Effect Regression Models, R package version 0.6.1.
Koštál, V. (2006). Eco-physiological phases of insect diapause. J Insect Physiol , 52, 113–127.
Kuznetsova, A., Brockhoff, P.B. & Christensen, R.H.B. (2017). lmerTest Package: Tests in Linear Mixed Effects Models. J Stat Softw , 82, 1–26.
Levy, R.C., Kozak, G.M., Wadsworth, C.B., Coates, B.S. & Dopman, E.B. (2015). Explaining the sawtooth: latitudinal periodicity in a circadian gene correlates with shifts in generation number. J Evol Biol , 28, 40–53.
Lindestad, O., Wheat, C.W., Nylin, S. & Gotthard, K. (2019). Local adaptation of photoperiodic plasticity maintains life cycle variation within latitudes in a butterfly. Ecology , 100, e02550.
Lüdecke, D. (2018). ggeffects: Tidy Data Frames of Marginal Effects from Regression Models. J Open Source Softw , 3, 772.
Macgregor, C.J., Thomas, C.D., Roy, D.B., Beaumont, M.A., Bell, J.R., Brereton, T., et al. (2019). Climate-induced phenology shifts linked to range expansions in species with multiple reproductive cycles per year. Nat Commun , 10, 4455.
Masaki, S. (1961). Geographic variation of diapause in insects. 弘前大学農学部学術報告, 66–98.
Matechou, E., Dennis, E.B., Freeman, S.N. & Brereton, T. (2014). Monitoring abundance and phenology in (multivoltine) butterfly species: a novel mixture model. J Appl Ecol , 51, 766–775.
McNamara, J.M., Barta, Z., Klaassen, M. & Bauer, S. (2011). Cues and the optimal timing of activities under environmental changes. Ecol Lett , 14, 1183–1190.
Michielini, J.P., Dopman, E.B. & Crone, E.E. (2021). Changes in flight period predict trends in abundance of Massachusetts butterflies.Ecol Lett , 24, 249–257.
Nielsen, M.E. & Kingsolver, J.G. (2020). Compensating for climate change–induced cue-environment mismatches: evidence for contemporary evolution of a photoperiodic reaction norm in Colias butterflies.Ecol Lett , 23, 1129–1136.
Nielsen, M.E., Lehmann, P. & Gotthard, K. (2022). Longer and warmer prewinter periods reduce post-winter fitness in a diapausing insect.Funct Ecol , 36, 1151–1162.
Opler, P.A. & Krizek, G.O. (1984). Butterflies east of the Great Plains: an illustrated natural history . Johns Hopkins Univ. Press, Baltimore, Maryland.
Orme, D., Freckleton, R., Thomas, G., Petzoldt, T., Fritz, S., Isaac, N., et al. (2023). caper: Comparative Analyses of Phylogenetics and Evolution in R, R package version 1.0.3.
Paradis, E. & Schliep, K. (2019). ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics , 35, 526–528.
Parmesan, C. (2006). Ecological and Evolutionary Responses to Recent Climate Change. Annu Rev Ecol Evol S , 37, 637–669.
Phillimore, A.B., Stålhandske, S., Smithers, R.J. & Bernard, R. (2012). Dissecting the Contributions of Plasticity and Local Adaptation to the Phenology of a Butterfly and Its Host Plants. Am Nat , 180, 655–670.
van de Pol, M. & Wright, J. (2009). A simple method for distinguishing within- versus between-subject effects using mixed models. Anim Behav , 77, 753–758.
Pollard, E. & Yates, T.J. (1993). Monitoring butterflies for ecology and conservation: the British butterfly monitoring scheme . Chapman & Hall, London.
Pöyry, J., Leinonen, R., Söderman, G., Nieminen, M., Heikkinen, R.K. & Carter, T.R. (2011). Climate-induced increase of moth multivoltinism in boreal regions: Climate-induced increase in moth multivoltinism.Global Ecol Biogeogr , 20, 289–298.
R Core Team. (2023). R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing, Vienna, Austria.
Roy, D.B., Oliver, T.H., Botham, M.S., Beckmann, B., Brereton, T., Dennis, R.L.H., et al. (2015). Similarities in butterfly emergence dates among populations suggest local adaptation to climate.Glob Change Biol , 21, 3313–3322.
Schlaepfer, M.A., Runge, M.C. & Sherman, P.W. (2002). Ecological and evolutionary traps. Trends Ecol Evol , 17, 474–480.
Scrucca, L., Fop, M., Murphy, T.B. & Raftery, A.E. (2016). mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models. R J , 8, 289–317.
Suggitt, A.J., Wilson, R.J., Isaac, N.J.B., Beale, C.M., Auffret, A.G., August, T., et al. (2018). Extinction risk from climate change is reduced by microclimatic buffering. Nature Clim Change , 8, 713–717.
Tauber, M.J. & Tauber, C.A. (1976). Insect Seasonality: Diapause Maintenance, Termination, and Postdiapause Development. Annu Rev Entomol , 21, 81–107.
Thornton, M.M., Shrestha, R., Wei, Y., Thornton, P.E., Kao, S.-C. & Wilson, B.E. (2022). Daymet: Daily Surface Weather Data on a 1-km Grid for North America, Version 4 R1.
Urbanski, J., Mogi, M., O’Donnell, D., DeCotiis, M., Toma, T. & Armbruster, P. (2012). Rapid Adaptive Evolution of Photoperiodic Response during Invasion and Range Expansion across a Climatic Gradient.Am Nat , 179, 490–500.
Van Deynze, B., Swinton, S.M., Hennessy, D.A., Haddad, N.M. & Ries, L. (2022). Neonicotinoids, more than herbicides, land use, and climate, drive recent butterfly declines in the American Midwest. bioRxiv .
Van Dyck, H., Bonte, D., Puls, R., Gotthard, K. & Maes, D. (2015). The lost generation hypothesis: could climate change drive ectotherms into a developmental trap? Oikos , 124, 54–61.
Van Klink, R., Bowler, D.E., Gongalsky, K.B., Swengel, A.B., Gentile, A. & Chase, J.M. (2020). Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances. Science , 368, 417–420.
Wagner, D.L., Grames, E.M., Forister, M.L., Berenbaum, M.R. & Stopak, D. (2021). Insect decline in the Anthropocene: Death by a thousand cuts.Proc Natl Acad Sci USA , 118, e2023989118.
Warren, M., Hill, J., Thomas, J., Asher, J., Fox, R., Huntley, B.,et al. (2001). Rapid responses of British butterflies to opposing forces of climate and habitat change. Nature , 414, 65–69.
Warren, M.S., Maes, D., van Swaay, C.A.M., Goffart, P., Van Dyck, H., Bourn, N.A.D., et al. (2021). The decline of butterflies in Europe: Problems, significance, and possible solutions. Proc Natl Acad Sci USA , 118, e2002551117.
Wei, T. & Simko, V. (2021). R package “corrplot”: Visualization of a Correlation Matrix, R package version 0.92.
Wenger, S.J., Stowe, E.S., Gido, K.B., Freeman, M.C., Kanno, Y., Franssen, N.R., et al. (2022). Simple statistical models can be sufficient for testing hypotheses with population time-series data.Ecol Evol , 12, e9339.
Wepprich, T., Adrion, J.R., Ries, L., Wiedmann, J. & Haddad, N.M. (2019a). Butterfly abundance declines over 20 years of systematic monitoring in Ohio, USA. PLOS ONE , 14, e0216270.
Wepprich, T., Adrion, J.R., Ries, L., Wiedmann, J. & Haddad, N.M. (2019b). Data from: Butterfly abundance declines over 20 years of systematic monitoring in Ohio, USA. Dryad Digital Repository. Available at: https://doi.org/10.5061/dryad.cf78420.
Wepprich, T., Henry, E. & Haddad, N.M. (2024). Data from: Decades of butterfly monitoring reveal adaptation of multivoltine species to climate warming. Zenodo. Available at: https://zenodo.org/doi/10.5281/zenodo.10607626.
Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L.D., François, R., et al. (2019). Welcome to the tidyverse. J Open Source Softw , 4, 1686.
Wood, S.N. (2017). Generalized additive models: an introduction with R . Chapman & Hall/CRC, Boca Raton, FL.

Figure legends

Fig. 1: Conceptual diagram of population consequences from voltinism changes.

Lost generations are predicted when increased voltinism is mismatched with winter (onset or severity) and would potentially decrease fitness and population growth rates. There are three other potential scenarios for species with voltinism change confronting various winter conditions, including matches that would increase population growth rates. A mismatch when species do not increase voltinism enough for warming winters is coined “missing” generations as a counterpoint to lost generations.

Fig. 2: Phenological variation for an example species (Phyciodes tharos).

Seasonal phenology was modeled from weekly counts to predict relative abundance on a degree-day timescale, with example model predictions from warmer or cooler sites during warmer or cooler years (A). When these same model predictions are visualized on a calendar date timescale, phenological shifts in response to site and annual temperature misalign generations (B). Mixture models assign counts to different generations on the degree-day timescale, with example shown from all sites during warm and cool years, split by warm and cool sites (C). The log population growth rate between the penultimate and last generations is the measure of last generation size, which varies in part due to photoperiod exposure that varies by latitude and phenology of the sensitive life stage (D). Years in this example (2008-2017) are split in half into warm or cool years. Sites in this example are split into warm or cool by southern and northern regions of Ohio.

Fig. 3: Voltinism trait estimates and long-term population trends.

Species’ estimates and uncertainty (95% confidence intervals) for four voltinism traits and their statewide population trend. Solid black coloration indicates that the estimate’s uncertainty does not overlap zero, while gray coloration indicates otherwise. Species are ordered by the population trend estimate. LG = “last generation”.

Fig. 4: Overwinter population growth rate across all 30 species.

Overwinter population growth rate (y-axis) is influenced by the interaction of last generation size (annual variation and site mean), species phenology (earlier-, mid-, or later-season), and winter temperature. Annual variation in last generation size on the x-axis shows increased voltinism at larger numbers. The three rows of species phenology shows the model prediction if the penultimate generation peaks on July 2 (top), August 1 (middle), or August 31 (bottom). The two columns are divided by cooler winter temperatures (left) and warmer winter temperatures (right), scaled by +/- 1 SD of the annual variation in winter temperature. Site variation is shown by separate lines in each panel for warmer sites in red and cooler sites in blue (+/- 1 SD of site mean last generation size).

Fig. 5: Correlation of voltinism traits and population trends.

A) A correlation matrix of six voltinism traits and population trends are based on Pearson correlations between variables performed with a parametric bootstrap. Values are the median bootstrap correlation with 95% bootstrap confidence intervals in parentheses. Colors indicate the direction (positive or negative) and magnitude (circle size) of the correlation.
B) Grouping species by the trend in last generation size and the simulated consequence of larger generations shows 4 zones of the potential scenarios in Fig. 1. The lower right quadrant would fit the lost generation hypothesis. The upper right quadrant would fit a “bonanza” scenario of matching voltinism to a warming climate. The upper left quadrant would also be a mismatch scenario with species not taking advantage of a warming climate by increasing voltinism. The lower left quadrant would be species not increasing voltinism, an appropriate match to the detrimental effect of attempting larger last generations.

Figures

Fig. 1: Conceptual diagram of population consequences from voltinism changes.