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.