\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Changes in Flight Period Predict Trends in Abundance of Massachusetts
Butterflies}



\author[1]{James Michielini}%
\author[1]{Erik Dopman}%
\author[1]{Elizabeth Crone}%
\affil[1]{Tufts University}%


\vspace{-1em}



  
  \date{May 19, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Phenological shifts are well-documented in the ecological literature.
However, their significance for changes in demography and abundance is
less clear. We used 27 years of citizen science monitoring to quantify
trends in phenology and relative abundance across 89 butterfly species.
We calculated shifts in phenology using quantile regression and shifts
in relative abundance using list length analysis and counts from club
trips. Elongated activity periods within a year were the strongest
predictor of increases in relative abundance. These changes may be
driven in part by changes in voltinism, as this association was stronger
in multivoltine species. Some species appear to be adding a late-season
generation while other species appear to be adding a spring generation,
revealing a possible shift from vagrant to resident. Our results
emphasize the importance of evaluating phenological changes throughout
species' activity periods and understanding the consequences for such
climate-related changes on viability or population dynamics.%
\end{abstract}%



\sloppy


\textbf{Introduction}

Phenological shifts, changes the timing of species life history events,
have become one of the best-documented signatures of anthropogenic
climate change (Forrest \& Miller-Rushing, 2010; Parmesan \& Yohe,
2003). Many prominent studies of phenology have focused on the timing of
reproduction, especially in plants (Parmesan \& Yohe, 2003). Some have
documented shifts in seasonal migration (M\selectlanguage{ngerman}øller et al., 2008). Others
have documented changes in emergence from or entry into dormant life
stages, including hibernation in vertebrate animals (Inouye, 2000),
leaf-out periods for trees (Polgar \& Primack, 2011), and diapause in
insects (Bale \& Howard, 2010). Although the majority of past studies
have demonstrated advances in phenology associated with climate warming,
there is considerable variation, with some species delaying phenology
(Fric et al., 2020; Forister \& Shapiro, 2003; Roy \& Sparks, 2000).

Despite widespread documentation of changes in phenology, few empirical
studies have tested whether phenological changes are associated with
long-term trends in population size (Ramula et al., 2015). In
terrestrial warming experiments, plants with earlier phenology in warmed
plots tended to have higher biomass, growth, and/or reproduction,
suggesting that phenological advances represent adaptive responses to
environmental change (Cleland et al., 2012). Advanced phenology was also
positively associated with long-term population trends of plants in
Concord, Massachusetts (Willis et al., 2008) and may be beneficial for
insect populations that can avoid competition for resources or feed on a
wider range of vegetation (Rathcke \& Lacey, 1985). Nevertheless,
changing phenology has also been shown to harm species, if changes lead
to phenological mismatches with resources or interacting species (Both
et al., 2006).

For insects, changes in phenology often arise through changes in the
timing of entry into diapause, as well as timing of spring emergence
(Bale et al., 2002). Many insect species exhibit geographically variable
patterns of voltinism (number of broods within a year) based on location
within their range. Under a warming climate, populations which were
formerly thermally restricted in parts of their range may be capable of
producing an additional generation within the extended growing season
(Kozak et al., 2019; Grevstad \& Coop, 2015; Mitton \& Ferrenberg, 2012;
Altermatt, 2010; Tobin et al., 2008) which could be beneficial or
detrimental for populations. Species that can successfully add a
generation in the extended growing season may benefit from another bout
of reproduction, leading to higher overall population growth rates (Kerr
et al., 2020; Kerr et al., 2019), and the potential for more rapid
evolutionary responses to climate change (Chevin et al., 2010). In spite
of its potential benefits, an increase in voltinism can also be
detrimental to insects. Mismatched phenological cues could cause a
population to start an additional generation that fails to reach its
diapausing life stage before frost (Levy et al., 2015; VanDyck et al.,
2014). Such developmental traps can potentially cause populations to
decline, possibly dramatically, as the flight period increases (VanDyck
et al., 2014). Knowing how these different responses affect population
dynamics is necessary to understand the consequences of phenological
shifts, as well as the longer-term effect of increased growing season
length on insect populations.

Butterflies are known to be sensitive to changes in temperature and have
been widely used as models to study phenological change (Bale et al.,
2002). Phenological studies of butterflies have demonstrated advances in
adult emergence, which have been attributed to warming climates
(Forister \& Shapiro, 2003; Roy \& Sparks, 2000). In contrast to
numerous examples of phenological advances in the onset of butterfly
flight, few studies have investigated potential changes in the end of
flight. In the two studies that have investigated empirical patterns of
late-season phenology in butterfly communities, Zipf et al. (2017) and
Westwood \& Blair (2010) both observed later end dates of flight
activity, correlated with increasing temperatures. These studies are
noteworthy in phenological research because much less is known about
late-season phenology than early-season phenology (Gallinat et al.,
2015; Karlsson, 2014). Nonetheless, in spite of the large body of work
on butterfly phenology, we do not know whether phenological changes are
beneficial or deleterious responses to changing environments. One recent
study has demonstrated that advanced emergence in British Lepidoptera
was associated with significantly higher rates of demographic abundance
within multivoltine species (MacGregor et al., 2019). Some studies of
birds have also generally demonstrated associations between delayed
phenology and decreasing abundance trends (Saino et al., 2011; Møller et
al., 2008; Both et al., 2006). Understanding the population dynamics
associated with phenology is imperative to translating climate-related
phenological changes into their impacts for long-term viability.

One common feature of past research on phenology change is analysis
using simple metrics such as first, average, or last observation date.
Inferences from such metrics may be limited because first and last
observation dates have known biases (Miller-Rushing, 2008) and averages
do not capture changes throughout the activity period (Inouye et al.,
2019). In this study, we use quantile regression (Cade and Noon, 2003)
to evaluate phenological change throughout the activity period. Quantile
regression is a statistical modelling technique that enables us to
robustly estimate changes in the onset and end of adult flight. Like
linear regression, quantile regression generates a slope-intercept line
through a specific part of the distribution. Linear regression minimizes
squared deviations in the response variable around a trend, whereas
quantile regression minimizes the absolute deviations from a trend,
subject to the constraint that some proportion of the data be below the
line (e.g. the 0.1 quantile, or 10\textsuperscript{th} percentile, fits
a line that minimizes absolute deviations with 10\% of observations
below the line). Quantile regression uses the complete data set to fit
this line, and so it differs from fitting a line to the first or last
x\% of data points, an \emph{ad hoc} technique that has occasionally
been used in ecological studies of phenology (Zipf, et al., 2017; Brooks
et al., 2014; Polgar et al., 2013). To date, a handful of studies have
used quantile regression to study trends in phenology, mostly in the
context of bird migration., These studies have revealed that changes in
phenology are not uniform (e.g., Barton \& Sandercock, 2018; Gimesi et
al., 2012; Gordo et al., 2013). Although under-used, quantile regression
represents a formal and robust technique for evaluating changes in
phenology throughout the activity period.

Here, we quantify long-term trends in phenology and abundance of
butterflies using 27 years of citizen science records from the
Massachusetts Butterfly Club (hereafter MBC). Using these observational
data for Massachusetts butterflies, we test the relationship between
year-to-year changes in flight phenology and abundance. For each
species, we estimate abundance trends through time using list length
analysis, updating previous analyses of an earlier subset of the same
data (Breed et al., 2013), and analyses of counts individuals. We test
whether changes in the onset of activity, the end of activity, the
average date, and the flight period are associated with increases or
declines in abundance. We also evaluate relationships between life
history traits, trends in phenology, and trends in abundance using
structural equation modelling to elucidate the potential mechanistic
pathways. We compare our results to a recent study which documented
similar associations between phenology and abundance for British
(MacGregor et al., 2019). Finally, we discuss the implications of our
findings for insect population viability in their phenological response
to global change.

\textbf{Methods}

\emph{Study system} : The Massachusetts Butterfly Club (MBC) is a group
of citizen scientists which documents sightings of butterflies and
maintains a database of observation records. The MBC compiles records of
butterfly sightings from organized group trips and as reported
opportunistically by individuals. The current club database of records
from 1992-2018 includes 30,707 species lists from these reports
(hereafter referred to as trips) encompassing 162,670 sightings of 99
different butterfly species. We thinned this to 10,884 trips
encompassing 126,087 sightings of 89 species suitable for analysis using
multiple criteria. First, although the database includes historic
records used to build the Atlas of Massachusetts Butterflies, we
excluded records prior to 1992, the year in which the club started to
record large numbers of lists (i.e., \textgreater{} 1200
sightings/year). Next, following Breed et al. (2013), trips with
\textless{} 5 species seen were excluded, to minimize the effect of
anecdotal sightings (see \emph{Appendix S1 in Supporting Information} ).
Unlike Breed et al. (2013), we combined and retained records for two
species groups. Canadian Tiger Swallowtail \emph{Paplio canadensis} and
Eastern Tiger Swallowtail \emph{Paplio glaucus} records were combined
(rather than excluded) because they are known to hybridize and can
sometimes be difficult to visually distinguish. Spring
Azure~\emph{Celastrina ladon} and Summer Azure \emph{Celastrina
neglecta} were also combined because of a change in species reporting
within the MBC during the study period. With these two exceptions,
species were excluded if during the study period they underwent a
taxonomic shift or realignment, or if there was an outbreak of the
population uncharacteristic of conventional population dynamics
(following Breed et al., 2013; see \emph{Appendix S2} ). Species were
also excluded if there were fewer than 15 observations before the year
2000, to prevent outliers or a lack of available data from driving
trends (See~\emph{Appendix S4} ). MBC records were processed into a
database, generating 89 species-specific lists of all trips. Each
species file included presence/absence of the species on each trip, the
date of each trip, the number of individuals of that species seen on
each trip, and the total number of species seen on each trip (see
\emph{Appendix S1} ).

To interpret patterns of phenology and abundance change we compiled data
on two life history traits: voltinism and range type. Following Breed et
al. (2013), voltinism was categorized as either obligately univoltine (a
single generation per year) or multivoltine (two or more generations per
year). Categorization agreed with Breed et al. (2013), with a few
exceptions in which scatter plots of observation day of year supported
an update (see \emph{Appendix S4} ). Following Breed et al. (2013), we
categorized species' range type based on the center of their geographic
range distribution, relative to Boston, Massachusetts (the approximate
mid-latitude for the state of Massachusetts). Species were ordinally
identified as having ranges centered north, central, or south of the
latitude of Boston (species classified as northern would therefore be
closer to the southern limit of their range).

\emph{Trends in abundance and phenology:~}We estimated relative
abundance using `list length' analysis (Breed et al., 2013; Szabo et
al., 2010; Franklin, 1999). List length analysis uses the total number
of species observed in an outing, i.e. the list's length, as a proxy for
search effort. List length may be more reliable than observer hours,
observer miles, or number of observers because the probability of seeing
a species does not scale directly with these metrics for observation
(Szabo et al., 2010; Franklin, 1991). Relative abundance of each species
was estimated using a binomial generalized linear mixed model (GLMM)
with fixed effects of year and list length, and a presence/absence
response (see \emph{Appendix S2} ). We also included year as a random
effect to account for interannual variation in relative abundance of
species. The yearly slope coefficients of these models represent changes
in relative abundance, measured as the change in species presence on
lists over time. Analyses were implemented in R using the lme4 package
(R Core Team, 2018; Bates et al., 2014). To supplement inference from
list length analysis, we analyzed trends in the number of butterflies of
each species seen per trip. These counts are a very coarse metric of
abundance because they were collected idiosyncratically (not with a
uniform survey protocol). However, congruent results from analyzing
these counts and list length analysis would suggest that both are
reasonable metrics of abundance. Count models used Poisson GLMMs with
fixed effects of year and random effects of year. Models were fit using
the glmer() function from the package `lme4' (Bates et al., 2014).

We used quantile regression (Cade \& Noon, 2003) to estimate trends in
the onset, end, and total length of observations of each species across
years. Changes in the onset of flight were estimated as the slope of the
0.1 quantile of observation dates for each species as a function of
year. Therefore, negative values indicate advances in phenology and
positive values indicate delays in phenology. Changes in the end of
flight were estimated as the slope of the 0.9 quantile of observation
dates of each species as a function of year. Therefore, negative values
indicate advances and positive values indicate delays. The change in the
total flight period was estimated as the difference between the 0.9 and
0.1 quantile slopes, where negative values indicate shortening flight
periods and positive values indicate lengthening flight periods.
Quantile regression was carried out using the rq() function from the
package `quantreg' (Koenker, 2019). To make an explicit comparison
between the changes in phenology of Massachusetts and British
butterflies from Macgregor et al. (2019), we also use linear regression
to measure the annual change in mean observation date as a function of
time with year included as both a fixed and random effect. Negative
values again indicate advances in the mean observation date and positive
values indicate delays. We used the glmer() function from the package
`lme4' (Bates et al., 2014).

As a simple metric of the association between changes in phenology and
changes in abundance, we correlated metrics of phenological change
(slope for each species with respect to time of the 0.1 quantile, 0.9
quantile, flight period, and mean date) with changes in abundance (slope
for each species from list length analysis). Because all slopes were
measured with error, we verified the relationship with parametric
bootstrapping, i.e., resampling estimates for each species from a random
normal distribution using the true mean and standard deviation for each
species (10,000 bootstrap samples). We evaluated taxonomic effects as
predictors of all response variables. Trends in abundance and various
phenology metrics were not associated with taxonomic groups (at the
family or genus level), and these effects are not discussed further
(see~\emph{Appendix S3} ).

\emph{Structural Equation Models:~}We observed a strong association
between changes in flight period and changes in abundance (see Results).
This correlation could be the result of various causal factors. We
explored two of these factors using structural equation models (SEMs).
SEMs can be used to analyze the effect size and directionality of causal
pathways among the variables (Grace et al., 2015). Breed et al. (2013)
had previously found that species were more likely to be increasing if
they were multivoltine or if our study area (Massachusetts) was north of
the center of their range. Our \emph{a priori} model had a direct effect
of range type (northern, central, or southern) on voltinism and
abundance trends; voltinism had a direct effect on flight period trends;
and range type had a direct effect on flight period (See \emph{Appendix
S5}). We scaled each variable in the SEM to x\selectlanguage{english} = 0 and sd = 1 such that
each per unit change in the response is with respect to 1 change
standard deviation of the predictor. This scaling facilitates
comparisons between variables with different units.

SEM pathways were built and evaluated using the `lavaan' package in R
(Rosseel et al., 2017). After fitting our \emph{a priori} model, we
explored model fit using the modindices() function to look for
correlations that were not included in our initial hypothesis (following
Grace et al., 2015). We used this exploratory tool to add additional
ecologically plausible relationships among voltinism, range type, and
trends in flight period and abundance. We also interpreted the original
and updated models to understand the causal relationships between trends
in flight period and trends in abundance.

\textbf{Results}

\emph{Trends in Abundance and Phenology:~}Because list length analysis
is a metric of relative abundance, average trends in abundance are
expected to be centered around 0 across species; this expectation held
(mean = 0.003 , sd = 0.045, greatest decline = -0.137 (\emph{Nymphalis
vaualbum} Compton tortoiseshell; Fig 1a), greatest increase = 0.190
(\emph{Poanes zabulon} Zabulon skipper; Fig 1b)). Trends estimated from
count models were also broadly centered around 0 (mean = 0.012, sd =
0.057, greatest decline = -0.159 (\emph{Nymphalis vaualbum} Compton
tortoiseshell), greatest increase = 0.213 (\emph{Poanes zabulon} Zabulon
skipper)). Trends estimated from list length were strongly congruent
with trends estimated from counts (Pearson's r = 0.89, P \textless{}
0.001). Therefore, we focus the rest of our analyses on trends estimated
from list length analyses (see~\emph{Appendix S2} ).

Across 26 years of continuous observation, species were, on average,
advancing their onset of flight approximately 0.20 days per year (sd =
0.42, range = (advancing) -1.75 days per year for \emph{Battus
philenor~}(Pipevine swallowtail; Fig 1c), to (delaying) +1.22 days per
year for~\emph{Vanessa cardui} (Painted lady; Fig 1d)) with an estimated
75\% of species advancing their 0.1 quantile observations. Species were,
on average, delaying their 0.9 quantile day of year observation
approximately 0.04 days per year (sd = 0.50, range = -1.76 days per year
(\emph{Nymphalis vaualbum} Compton Tortoiseshell), +1.29 days per year
(\emph{Erynnis horatius} Horace Duskywing)). Approximately 58\% of
species delayed their 0.9 quantile observations. Species were, on
average, increasing their flight period (measured as the difference
between the 0.9 quantile trend and 0.1 quantile trend) by 0.16 days per
year (sd = 0.55 days per year, range = -1.96 days per year
for~\emph{Nymphalis vaualbum} (Compton tortoiseshell), +2.05 days per
year~\emph{Battus philenor} (Pipevine swallowtail). Approximately an
estimated 63\% of species increasing their observed flight period.

Trends in the onset of flight activity (0.1 quantile) were negatively
correlated with changes in relative abundance, i.e., advancing the start
of the flight period was associated with increasing abundance (r =
-0.27, t = -2.7 , df = 87, p = 0.009, 98.03\% of bootstrapped
correlations \textless{} 0; Fig 2a). Trends in the end of flight
activity (0.9 quantile) were positively correlated with changes in
relative abundance, i.e., species that were extending their activity
later into the year also tended to be increasing in abundance (r = 0.36,
t = 3.59, df = 87, p \textless{} 0.001; 99.07\% of bootstrapped
correlations \textgreater{} 0; Fig 2b). The strongest association in our
data was a positive correlation between the annual change in the flight
period and the annual change in the relative abundance (r = 0.54, t =
5.9, df = 87, p \textless{} 0.001; 99.88\% of bootstrapped correlations
\textgreater{} 0; Fig 2c) such that species elongating their total
flight time were increasing in relative abundance.

Mean flight dates were significantly changing in univoltine species
(slope of mean day of year observation vs time \selectlanguage{ngerman}± SE: -0.27 ± 0.16) but
not for multivoltine species (slope of mean day of year observation vs.
time ± SE: -0.03 ± 0.26. Overall, there was no significant association
between trend in mean date and trend in abundance for all species
combined (r = 0.061, t = 0.574, df = 87, p = 0.568; 93.5\% of
bootstrapped correlations \textgreater{} 0). Grouping by voltinism and
using abundance trends estimated using counts for direct comparison with
MacGregor et al. (2019), mean day of year trends were not significantly
associated with abundance trends for multivoltine species (r = -0.06, t
= -0.43, df = 50, p = 0.672; 74.18\% of bootstrapped correlations
\textless{} 0 ), but were marginally significantly associated with
abundance trends for univoltine species (r = 0.32, t = 2.00, df = 35, p
= 0.053; 88.75\% of bootstrapped correlations \textgreater{} 0).

\emph{Structural equation models:} Our \emph{a priori~}structural
equation model (Fig. 3) showed significant direct effects of flight
period trends (scaled regression coefficient ± SE, \selectlanguage{greek}β \selectlanguage{english}= 0.472 \selectlanguage{ngerman}± 0.087)
and species' range type (\selectlanguage{greek}β \selectlanguage{english}= 0.308 \selectlanguage{ngerman}± 0.101) on trends in abundance,
which presumably reflects other traits associated with northern vs.
southern species. The effect of flight period trends on abundance trends
was larger than the direct effects of range type. Differences in flight
period were significantly associated with differences in voltinism (\selectlanguage{greek}β \selectlanguage{english}=
0.430 \selectlanguage{ngerman}± 096). Compared to univoltine species, multivoltine species more
often increased their flight period (days/year ± SE: -0.13 ± 0.09 and
0.35 ± 0.0 for univoltine and multivoltine, respectively) and increased
in relative abundance (trend ± SE: -0.01 ± 0.01 and 0.02 ± 0.01 for
univoltine and multivoltine, respectively). Differences in voltinism
were weaker but significantly associated with range type (proportion
multivoltine ± SE: 0.69 ± 0.06 and 0.32 ± .10 for southern and northern
species, respectively). This \emph{a priori} model demonstrated
adequate, but not exceptional fit (\selectlanguage{greek}χ\selectlanguage{english}\textsuperscript{2} = 3.076, df = 2,
p = 0.215; recalling that, in SEMs p \textless{} 0.05 represents
significant inconsistency with data).

\emph{Post hoc} inspection of modification indices revealed one
additional factor that was not present in our \emph{a priori} model.
This missing factor was a relationship between range type and flight
period trends; species with southerly ranges were increasing flight
periods more than species with northerly ranges (days/year \selectlanguage{ngerman}± SE: 0.35 ±
0.07 and -0.13 ± 0.09, for southern and northern species, respectively),
regardless of voltinism. Including this relationship improved model fit
(\selectlanguage{greek}χ\selectlanguage{english}\textsuperscript{2} = 0.614, df = 1, p = 0.433), and slightly
decreased the strength of the association between voltinism and flight
period (updated \selectlanguage{greek}β \selectlanguage{english}= 0.381 \selectlanguage{ngerman}± 0.099).

\textbf{Discussion}

For butterflies in Massachusetts, changes in abundance were most
strongly associated with changes in flight period. Species that have
extended their flight period have increased in abundance relative to
those that have not. Past studies that have evaluated the effect of
warming temperatures have generally demonstrated that species advancing
phenological events under warming environmental conditions are generally
increasing in abundance or have positive demographic growth (Cleland et
al., 2014; Hulme, 2011; Ozgul et al. 2010; Willis et al., 2008; Møller
et al., 2008). Species which fail to change their phenology, however,
appear to be declining (MacGregor et al. 2019; Cleland et al., 2014;
Willis et al., 2008; Møller et al., 2008). In many ways, our results
parallel a recent study of Lepidoptera in the Britain (MacGregor et al.,
2019). In both studies, changes in some phenological metric were
associated with increases in abundance, especially in multivoltine
species. However, in contrast to Britain, trends in the mean timing of
activity of Massachusetts butterflies were not strongly associated with
trends in abundance. Instead, changes in the flight period were the best
predictor of trends in abundance. In Massachusetts, changes in flight
period were driven more by extended late-season activity than advanced
onset of activity, which also contrasts with past work in Britain (Roy
and Sparks, 2000). Our results highlight the importance of evaluating
the entire temporal distribution of a phenological phenomenon, as
opposed to simply using trends in mean observation dates.

One possible mechanistic explanation for the association between
extended flight period and increased abundance is an increase in the
number of generations per year. In other insect species, climate warming
may allow for faster larval development and the production of additional
direct-developing generations before adverse winter conditions
(MacGregor et al., 2019; Miton \& Ferrenberg, 2012; Altermatt, 2009;
Tobin et al., 2008). Visual inspection of representative species in our
data suggest that at least some populations are adding, at least, a
partial generation (Fig. 1b, c \& e). For example, \emph{Satyrdotes
appalachia} (Appalachian Brown, Fig. 1b) may be transitioning from
univoltine to bivoltine in our region (Stichter, 2015). If additional
generations are the causal mechanism between changes in flight period
and increased abundance in multivoltine species, it would suggest that
beneficial effects of higher population growth rates (Kerr et al., 2020)
and/or more rapid evolution (Chevin et al., 2010) outweigh the potential
costs of developmental traps (Van Dyck et al., 2015; Levy et al., 2015).
In general, multivoltine butterfly species seem to have higher
population growth rates than univoltine species in north-temperate
climates (MacGregor et al., 2019; Wepprich et al., 2019; Breed et al.,
2013). In one study from a Mediterranean climate, however, multivoltine
species showed significantly steeper declines than obligate uni- or
bivoltine species. In this ecosystem, seasonal adversity in summer
months via drought may have led to a developmental trap (Melero et al.,
2016).

Past studies of butterflies in Massachusetts (Breed et al., 2013) and of
other Lepidoptera (Pöyry, 2009; Parmesan et al., 1999) have revealed
strong patterns of poleward range shifting in response to climate
change. In Massachusetts, southern species are increasing in relative
abundance, and northern species tend to be declining (Breed et al.
2013). Our structural equation model suggests one causal pathway by
which southern species are increasing in relative abundance to some
extent because of their capacity to extend their overall flight period.
Although the SEM supports variability in voltinism as the strongest
driver of trends in abundance, there were also direct effects of
species' range type on trends in abundance. These effects presumably
reflect other limiting factors such as physiological tolerance or trends
in interacting species. For example, \emph{Nymphalis vaualbum} (Compton
Tortoiseshell, Fig. 1a) is a northern species that has recently declined
rapidly in the state of Massachusetts. The abrupt shift in abundance
suggests a physiological threshold could have been crossed in the past
decade.

For at least two species in our data, changes in flight period appeared
to be due in part to changes in the ability to overwinter in
Massachusetts, as opposed to voltinism~\emph{per se}.~\emph{Poanes
zabulon} Zabulon skipper (Fig. 1b) and~\emph{Battus philenor} Pipevine
swallowtail (Fig. 1c) were formerly considered vagrants. In other words,
in the past, individuals of these species moved annually to
Massachusetts from populations further south during the summer, but
never established resident populations (Stichter, 2015). It seems likely
that these species can now successfully overwinter in Massachusetts, as
shown by an apparent early season generation (see Fig. 1b \& c). These
trends contrast with two other migratory species,~\emph{Vanessa
cardui~}(Painted Lady, Fig. 1d) and~\emph{Danaus plexippus} (Monarch,
Fig. 1f), both which are non-resident, multi-generational, obligatory
migrants, far north of overwintering sites. Both species have delays
(not advances) in the onset of spring activity in Massachusetts, and
only~\emph{V. cardui} showed delayed end of flight activity, while there
was no substantial trend in the end of~\emph{D. plexippus} flight in
Massachusetts. Both species showed negligible changes in relative
abundance in Massachusetts.

Across the community of butterflies in Massachusetts, trends in
abundance estimated from counts, as well as trends from the list length
analysis, were centered on zero. In other words, approximately equal
numbers of Massachusetts species are increasing in abundance as are
declining in abundance. Our results contrast with negative average
trends reported across butterfly communities in four other regions:
Ohio, USA (Wepprich et al. 2019), California, USA (Forister et al.
2011), the UK (Brereton et al. 2018), and Spain (Melero et al. 2016).
Unlike Ohio and California, both of which have agricultural development
at industrial scales, Massachusetts may be subject to lower levels of
agricultural land use intensification. Systematic monitoring rather
than~\emph{ad hoc} citizen science observations would be needed to
rigorously test whether butterfly communities in Massachusetts are truly
more stable in abundance when compared to other regions. Nonetheless,
our data suggest a hopeful hypothesis that butterflies may be able to
adjust to climate change, if freed from the additional pressure of
anthropogenic changes in land use.

Our results highlight the potential for phenological changes to affect
population dynamics, and the general importance of phenological
flexibility for insect population dynamics in changing environments.
Past studies of phenology have rarely evaluated shifts in the end of
activity (Zipf al., 2017; Gallinat et al., 2015), despite the
demographic importance of late-season events such as entry into diapause
(Kerr et al. 2020). Unlike traditional metrics of phenological change,
quantile regression (Cade \& Noon, 2003) enabled us to estimate trends
in both the onset and end of flight across Massachusetts butterflies. In
Massachusetts, changes in the total flight period and end of activity
were associated more strongly with trends in abundance than onset or
mean trends. In Lepidoptera, a diversity of cues and mechanisms, both
genetic and plastic (Lindestad et al., 2019; Kozak et al., 2019), shape
patterns of phenology and voltinism. Understanding these mechanisms
would be a valuable direction for future research if our goal is to
predict changes in phenology and abundance in new systems. This case
study highlights how phenological changes relate to population dynamics,
and the importance of phenological changes for population viability in
changing environments.

\textbf{Acknowledgements}

The authors thank the Massachusetts Butterfly Club (especially Mark
Fairbrother, Sharon Stichter \& Steve Moore) for graciously providing us
with their extensive database records. This work was funded by a DOD
SERDP award (RC-2700) to E. Crone and NSF DEB award (1257251) to E.B.D.

\textbf{References}

Altermatt, F. (2010). Climatic warming increases voltinism in European
butterflies and moths.\emph{Proc. R Soc. Lond. {[}Biol{]},} 277,
1281--1287.

Bale, J. S. (2002). Insects and low temperatures: from molecular biology
to distributions and abundance.~\emph{Philos. Trans. R. Soc. B.} ,~357,
849-862

Bale, J. S., \& Hayward, S. A. L. (2010). Insect overwintering in a
changing climate.~\emph{J. Exp. Biol.} ,~213, 980-994.

Barton, G. G., \& Sandercock, B. K. (2018). Long-term changes in the
seasonal timing of landbird migration on the Pacific
Flyway.~\emph{Condor,} 120, 30-46.

Bates, D., Maechler, M., Bolker, B., \& Walker, S. (2014). lme4: Linear
mixed-effects models using Eigen and S4.~\emph{R package version} ,~1,
1-23.

Both, C., Bouwhuis, S., Lessells, C. M., \& Visser, M. E. (2006).
Climate change and population declines in a long-distance migratory
bird.~\emph{Nature} ,~441, 81.

Breed, G. A., Stichter, S. \& Crone, E. E. (2013). Climate-driven
changes in northeastern US butterfly communities. \emph{Nat Clim
Change,} 3, 142--145.

Brereton, T. M., Botham, M. S., Middlebrook, I., Randle, Z., Noble, D.,
\& Roy, D. B. (2014). United Kingdom butterfly monitoring scheme report
for 2014.~\emph{Centre for Ecology and Hydrology and Butterfly
Conservation} .

Brooks, S. J., Self, A., Toloni, F., \& Sparks, T. (2014). Natural
history museum collections provide information on phenological change in
British butterflies since the late-nineteenth century.~\emph{Int. J.
Biometeorol.} ,~58, 1749-1758.

Cade, B. S. \& Noon, B. R. (2003). A gentle introduction to quantile
regression for ecologists.\emph{Front. Ecol. Envrion.,} 1, 412-420.

Chevin, L. M., \& Lande, R. (2010). When do adaptive plasticity and
genetic evolution prevent extinction of a density\selectlanguage{english}-regulated
population?~\emph{Evolution,} ~64, 1143-1150.

Cleland, E.E., Allen J. M., Crimmins, T.M., Dunne, J.A., Pau, S.,
Traver, S.E., Zavaleta, E.S., Wolkovich, E.M. (2012). Phenological
tracking enables positive species responses to climate
change.\emph{Ecology,} 93, 1765-1771.

Fric, Z. F., Rindo\selectlanguage{polish}š, \selectlanguage{english}M., \& Konvi\selectlanguage{polish}č\selectlanguage{english}ka, M. (2020). Phenology responses of
temperate butterflies to latitude depend on ecological
traits.~\emph{Ecol. Lett.} , 23, 172-180.

Forister, M. L., Shapiro, A. M. (2003). Climatic trends and advancing
spring flight of butterflies in lowland California. \emph{Glob. Change
Biol.,} 9, 1130--1135.

Forister M. L., Jahner J. P., Casner K. L., Wilson J. S., Shapiro A. M.
(2011). The race is not to the swift: long-term data reveal pervasive
declines in California's low-elevation butterfly fauna. \emph{Ecology} ,
92, 2222--2235

Franklin, D. C. (1999). Evidence of disarray amongst granivorous bird
assemblages in the savannas of northern Australia, a region of sparse
human settlement. \emph{Biol. Cons.,} 90, 53--68.

Forrest, J., \& Miller-Rushing, A. J. (2010). Introduction: Toward a
synthetic understanding of the role of phenology in ecology and
evolution.~\emph{Philos. Trans. R. Soc. B.} , 3101-3112.

Gallinat, A. S., Primack, R. B., \& Wagner, D. L. (2015). Autumn, the
neglected season in climate change research.~\emph{Trends Ecol. Evolut.}
,~30, 169-176.

Gimesi, L., Homor\selectlanguage{ngerman}ódi, R., Hirka, A., Szaboki, C., \& Hufnagel, L.
(2012). The effect of climate change on the phenology of moth abundance
and diversity.~\emph{Appl. Ecol. Environ. Res.,} 10, 349-363.

Gordo, O., Tryjanowski, P., Kosicki, J. Z., \& Fulín, M. (2013). Complex
phenological changes and their consequences in the breeding success of a
migratory bird, the white stork Ciconia ciconia.~\emph{J. Anim. Ecol.,}
~82, 1072-1086.

Grace, J. B., Scheiner, S. M., \& Schoolmaster Jr, D. R. (2015).
Structural equation modeling: building and evaluating causal models:
Chapter 8. (\emph{Ecological Statistics: Contemporary theory and
application} ) Oxford Scholarship Online, pp. 168-199.

Grevstad, F. S., \& Coop, L. B. (2015). The consequences of
photoperiodism for organisms in new climates.~\emph{Ecol. Appl.} ,~25,
1506-1517.

Hulme, P. E. (2011). Contrasting impacts of climate\selectlanguage{english}-driven flowering
phenology on changes in alien and native plant species
distributions.~\emph{New Phytol.} , 189, 272-281.

Inouye, D. W., Barr, B., Armitage, K. B., \& Inouye, B. D. (2000).
Climate change is affecting altitudinal migrants and hibernating
species.~\emph{PNAS} ,~97, 1630-1633.

Inouye, B. D., Ehrlen, J., \& Underwood, N. (2019). Phenology as a
process rather than an event: from individual reaction norms to
community metrics.~\emph{Ecol. Monograph} ,~89, e01352.

Karlsson, B. (2014). Extended season for northern
butterflies.~\emph{Int. J. Biometeorol.} , 58, 691-701.

Kerr, N. Z., Wepprich, T., Grevstad, F., 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. \emph{Glob. Change Biol.,} 26, 2014-2027\emph{.}

Kerr, N. Z.,~Crone, E. E.,~Chew, F. S.~(2019). Life history trade\selectlanguage{english}-offs
are more pronounced for a noninvasive, native butterfly compared to its
invasive, exotic congener.~\emph{Pop. Ecol} ., 62: 119--~133.

Koenker, R., Portnoy, S., Ng, P. T., Zeileis, A., Grosjean, P., \&
Ripley, B. D. (2019). Package `quantreg'.

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. \emph{J
Evol. Biol.,} 28, 40--53.

Lindestad, O., Wheat, C.W., Nylin, S. and Gotthard, K., (2019). Local
adaptation of photoperiodic plasticity maintains life cycle variation
within latitudes in a butterfly. \emph{Ecology} , 100, p.e02550.

Macgregor, C. J., Thomas, C. D., Roy, D. B., Beaumont, M. A., Bell, J.
R., Brereton, et al. (2019). Climate-induced phenology shifts linked to
range expansions in species with multiple reproductive cycles per
year.~\emph{Nature,} ~10, 1-10.

Melero, Y., Stefanescu, C., \& Pino, J. (2016). General declines in
Mediterranean butterflies over the last two decades are modulated by
species traits.~\emph{Biol. Cons.,} ~201, 336-342.

Mitton, J. B., \& Ferrenberg, S. M. (2012). Mountain pine beetle
develops an unprecedented summer generation in response to climate
warming.~\emph{Amer. Nat.,} ~179, 163-171.

Moller, A. P., Rubolini, D. \& Lehikoinen, E. (2008). Populations of
migratory bird species that did not show a phenological response to
climate change are declining. \emph{PNAS,} 105, 16195-16200.

Ozgul, A., Childs, D. Z., Oli, M. K., Armitage, K. B., Blumstein, D. T.,
Olson, L., et al. (2010). Coupled dynamics of body mass and population
growth in response to environmental change.~\emph{Nature,} ~466,
482-485.

Parmesan, C., Ryrholm, N., Stefanescu, C., Hill, J. K., Thomas, C. D.,
Descimon, H., et al. (1999). Poleward shifts in geographical ranges of
butterfly species associated with regional warming.~\emph{Nature,} ~399,
579-583.

Parmesan, C., \& Yohe, G. (2003). A globally coherent fingerprint of
climate change impacts across natural systems.~\emph{Nature} ,~421,
37-42

Polgar, C. A., \& Primack, R. B. (2011). Leaf\selectlanguage{english}-out phenology of temperate
woody plants: from trees to ecosystems.~\emph{New phytol.} ,~191,
926-941.

Polgar, C. A., Primack, R. B., Williams, E. H., Stichter, S., \&
Hitchcock, C. (2013). Climate effects on the flight period of Lycaenid
butterflies in Massachusetts.~\emph{Biol. Cons.,} ~160, 25-31.

Poyry, J., Luoto, M., Heikkinen, R. K., Kuussaari, M., \& Saarinen, K.
(2009). Species traits explain recent range shifts of Finnish
butterflies.~\emph{Glob. Change Biol.,} ~15, 732-743.

R Core Team (2019). \emph{R: A language and environment for statistical
computing} . R Foundation for Statistical Computing, Vienna, Austria.
URL https://www.R-project.org/

Ramula, S. Johansson, J., Linden, A. \& Jonzen, N. (2015). Linking
phenological shifts to demographic change. \emph{Clim. Res.} , 63,
135-144.

Rathcke, B., \& Lacey, E. P. (1985). Phenological patterns of
terrestrial plants.~\emph{Annu. Rev. Ecol. Evol. Syst.} ,~16, 179-214.

Roy D. B. \& Sparks, T. H. (2000). Phenology of British Butterflies and
climate change. \emph{Glob. Change Biol.,} 6, 76-81.

Saino, N., Ambrosini, R., Rubolini, D., von Hardenberg, J., Provenzale,
A., Huppop, K., et al. (2011). Climate warming, ecological mismatch at
arrival and population decline in migratory birds.~\emph{Proc R Soc Lond
{[}Biol{]},} ~278, 835-842.

Stichter, S. (2015).\emph{Butterflies of Massachusetts} . Available at:
https://www.butterfliesofmassachusetts.net/index.htm. Last accessed 28
January 2020.

Szabo, J. K., Vesk, P. A., Baxter, P. W. J. \& Possingham, H. P. (2010).
Regional avian species declines estimated from volunteer-collected
long-term data using List Length Analysis. \emph{Ecol. Appl.,} 20,
2157--2169.

Tobin, P. C., Nagarkatti, S., Loeb, G., \& Saunders, M. C. (2008).
Historical and projected interactions between climate change and insect
voltinism in a multivoltine species.~\emph{Glob. Change Biol.} ,~14,
951-95.

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?~\emph{Oikos} ,~124, 54-61.

Wepprich, T., Adrion, J. R., Ries, L., Wiedmann, J., \& Haddad, N. M.
(2019). Butterfly abundance declines over 20 years of systematic
monitoring in Ohio, USA.~\emph{PloS ONE} , 14 e0216270.

Westwood, A. R., \& Blair, D. (2010). Effect of regional climate warming
on the phenology of butterflies in boreal forests in Manitoba,
Canada.~\emph{Environ. Entomol.} ,~39, 1122-1133.

Willis, C. G., Ruhfel, B., Primack, R. B., Miller-Rushing, A. J. \&
Davis, C. C. (2008). Phylogenetic patterns of species loss in Thoreau's
woods are driven by climate change. \emph{PNAS,} 105, 17029--17033.

Zipf, L., Williams, E. H., Primack, R. B. \& Stichter, S. (2017).
Climate effects on late-season flight times of Massachusetts
butterflies. \emph{Int. J. Biometeorol.,} 61, 1667--1673.

\par\null

\textbf{Figure legends}

\textbf{}

\textbf{Figure 1} Paired day of year observations and change in presence
for select species over the study period (1992-2018). For each species,
the top plot graphs observations by the year and numeric day of year of
every sighting throughout the study period, slightly jittered for
visibility; The upper line represents the end of flight (estimated as
the 0.9 quantile day of year observation) and the bottom line represents
the onset of flight (estimated as the 0.1 quantile day of year
observation). For each species, the bottom plot describes the change in
abundance over the study period; each point represents the total
proportion of trips on which the species was present for each year, the
slopes of these plots reflect trends in relative abundance incorporated
into the list length model. Axis ticks for all plots represent the years
2000 and 2010.

\par\null

\textbf{Figure 2.} Correlations between annual changes in relative
abundance and 0.1 quantile (a), 0.9 quantile (b), flight period (c).
Each point (a, b \& c) represents a single species, shaped by voltinism
and colored by range type. Y axis values represent the annual change in
relative abundance estimated from list length analysis. (a) change in
start of flight is measured as the slope of the 0.1 quantile (advances
or delays); (b) change in the end of flight is measured as the slope of
the 0.9 quantile (advances or delays); (c) change in total flight period
is measured as the difference between the 0.9 and 0.1 quantiles
(elongations or contractions).

\par\null

\textbf{Figure 3.} Structural equations models. Lines indicate the
direction of causality. Listed adjacent to lines are the scaled
regression coefficients, the relative magnitude of the effects with
respect to one changing standard deviation of the predictor. The \emph{a
priori} model is represented by solid lines, and \emph{post
hoc}modifications are represented by dashed lines.

\emph{}

\emph{Figure 1.}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig1-a-850x450/Fig1-a-820x520}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/Fig1-b-820x520}

\end{center}
\end{figure}

\emph{}

\emph{Figure 2.}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/Fig-2-Authorea}

\end{center}
\end{figure}

\emph{Figure 3}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.42\columnwidth]{figures/Fig-3-Authorea/Fig-3-Authorea}
\caption{{This is a caption
{\label{534704}}%
}}
\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

