Combining range and phenology shifts offers a winning strategy for boreal Lepidoptera

Species can adapt to climate change by adjusting in situ or by dispersing to new areas, and these strategies may complement or enhance each other. Here, we investigate temporal shifts in phenology and spatial shifts in northern range boundaries for 289 Lepidoptera species by using long- term data sampled over two decades. While 40% of the species neither advanced phenology nor moved northward, nearly half (45%) used one of the two strategies. The strongest positive population trends were observed for the minority of species (15%) that both advanced flight phenology and shifted their northern range boundaries northward. We show that, for boreal Lepidoptera, a combination of phenology and range shifts is the most viable strategy under a changing climate. Effectively, this may divide species into winners and losers based on their propensity to capitalize on this combination, with potentially large consequences on future community composition.


INTRODUCTION
Biodiversity is threatened by a multitude of anthropogenic factors, including ongoing climate change (Díaz et al., 2019;IPCC, 2014). Under environmental change, species have two avenues for escaping decline and consequent extinction: adapting in situ through plastic or evolutionary responses or moving to areas where conditions are more favourable (Davis et al., 2005). Phenology and range shifts are the most conspicuous responses of species to rapid environmental change (Parmesan & Yohe, 2003;Root et al., 2003). The spatial distribution of a species and the temporal manifestation of life-history events reflect its realized niche (Socolar et al., 2017). Shifts in distribution and phenology represent the mechanisms through which a spatial or temporal change in the utilization of the niche can be observed (Amano et al., 2014). Shifts in phenology are often related to positive population trends and increased demographic stability (Cleland et al., 2007;Franks et al., 2018;Møller et al., 2008;Saino et al., 2011). Species that are able to shift their geographic distribution to track climate change and thus remain within their climatic niches are less likely to suffer overall population declines across their distribution (Cooper et al., 2011;Devictor et al., 2012;Urban, 2015).
Both range and phenology shifts offer species a means to track the thermal niche as climate changes and remain within the optimal temperature during important life-history events. Indications for complementarity in thermal niche tracking via space or time has been found for both plants and birds (Amano et al., 2014;Socolar et al., 2017). That species would use either range or phenology shift to respond to changing climatic conditions (Hypothesis 1; Figure 1a) is likely explained by constraints in their abilities to adjust phenological timing versus their abilities to disperse. Mobile species would not have as strong a need to advance their phenology as they can track suitable conditions in space. Less mobile species, however, will experience stronger selection pressure to adapt in situ and adjust their phenology in order to maintain suitable thermal conditions during critical life-history events (Amano et al., 2014;Socolar et al., 2017). If both of these responses increase population viability on their own, using either strategy should be reflected in overall stable or positive population trends. Potentially more optimal, however, would be to both advance phenology and move to cooler regions (Hypothesis 2; Figure 1a), as the strategies can provide enhancing benefits that together buffer species against climate warming and increased variation in extreme events. The ability to adapt in one dimension of niche utilization is likely to be correlated with high responsiveness to variability in others, not least through positive feedback loops that increase population persistence (Willis et al., 2010). Phenological timing is perhaps the most important aspect of life history that affects species distributions (Chuine, 2010), as it defines where and how successfully individuals of a population can proliferate. Hypothesis 2 is thus based on the assumption that adaptive in situ responses in phenology increase the fitness of individuals, leading to higher survival rates and more offspring (Cleland et al., 2007). Evidence suggests that stable or positive population trends, that is, no change or increase in abundance, are a prerequisite for species to expand their ranges (Mair et al., 2014) as emigration is higher from larger populations (Glorvigen et al., 2013;Pärn et al., 2012). Also, the probability of successful establishment increases with the summed contribution of individuals from neighbouring source populations (Hanski & Ovaskainen, 2003;Hanski et al., 1995). Successful colonization of new available habitat would further increase the probability of survival and reproductive success of individuals, which should have a positive effect on species abundance. If this is the case, combining both responses should offer more positive effects on population size than using only one of the strategies. Finally, a species that is not able to make use of either of the strategies (Hypothesis 0; Figure 1b) is expected to show declines in abundance.
In this study, we assess how 289 species of Lepidoptera in Finland have responded to ca. 20 years of climate change. Lepidoptera have been shown to be responsive to climate change, as exemplified by the observed range shifts towards higher latitudes (Chen et al., 2011;Parmesan et al., 1999;Pöyry et al., 2009) and altitudes (Konvicka et al., 2003;Wilson et al., 2005), as well as phenology shifts (e.g., advance in the date of first appearance [Roy & Sparks, 2000;Stefanescu et al., 2003;Diamond et al., 2011] and increased voltinism [Altermatt, 2010;Pöyry et al., 2011]). Previous studies have shown that the distribution of butterflies in Finland is mainly determined by climatic factors (Luoto et al., 2006) and that the phenology of moths tends to be controlled by temperature (Valtonen et al., 2014). During the last few decades, the mean temperature in Finland has risen by 0.2 to 0.4°C per decade (Mikkonen et al., 2014; Figure  S1), with springs starting earlier and the timing of phenological events being advanced for many species (Hällfors et al., 2020;Helama et al., 2020). Therefore, we expect Lepidoptera to respond to changing climatic conditions by either shifting their ranges or adjusting their phenology in situ, or both. However, it is unclear how these different strategies influence species' population trends. Here, we ask whether Lepidoptera have advanced their mean flight period (=adaptive phenology shift) or shifted their northern range boundary (NRB) northwards (=adaptive range shift). Using these response estimates, we test two main hypotheses: do the same species tend to (1) either advance their mean flight period or shift their NRB northwards, or (2) both advance their mean flight period and shift their NRBs northwards (Figure 1b). To gain understanding on the potential causes and effects of the response combinations, we further relate the responses to life-history traits and population trends. Although we cannot test for causality between traits, responses, and population trends, such an approach can reveal potential links between phenomena that could later further be tested experimentally.

Measures for phenology and range shifts
The start of the growing season in Europe, including Finland, has advanced over the past decades (Helama et al., 2020;Menzel & Fabian, 1999), and suitable | 1621 F I G U R E 1 Chart describing processes and predictions of the hypotheses. Panel (a) describes the underlying processes that may give rise to the patterns predicted by the outlined hypotheses. Underpinning Hypothesis 1 (either phenology or range shift) is the assumption that species differ fundamentally in their abilities to adjust either in situ or via dispersal. Assuming that these strategies are adaptive, being able to use either strategy will lead to an increased probability of presence, which should be reflected in positive population trends. Positive feedback loops through larger population size further enhance the ability of both strategies to function. Underpinning Hypothesis 2 (both phenology and range shift), on the other hand, is the assumption that adaptive in situ responses in phenology increases the fitness of the individuals, leading to higher rates of survival and/or more offspring. This in turn increases the probability of presence (stronger population trends) and thus higher colonization rates which leads to the species being able to expand into habitats becoming suitable as climate changes (=shift in the northern range boundary [NRB]). A successful colonization of new available habitat further increases the probability of survival and reproductive success of individuals, which again has a positive effect on species abundance. In this study, the hypothesized underlying processes are investigated through proxies for range shift, phenology shift, and probability of presences as depicted by derived estimates in yellow versus blue font in the process charts: shift in NRB as a measure of species range shift; change in adult flight period as a proxy for phenology shift; and population trends as a proxy for probability of presence across the distribution. Panel (b) describes the expected patterns in the data, i.e., the combinations of responses, as regards NRB and phenology shift estimates, that would support Hypotheses 1 (either advanced phenology or northwards shifting NRB), 2 (both advanced phenology and northwards shifting NRB), and 0 (neither advanced phenology nor northwards shifting NRB). Although these proxies do not allow us to infer evidence for the underlying processes, they can inform us of the patterns across a wide sample of species. By combining them with information on population trends, we can infer how successful the strategies likely are on their own and in combination for species experiencing climate change, and what may be the consequences if species cannot utilize either of the strategies conditions for seasonally recurring life-history events now often occur earlier than before. To investigate the strategy of adjusting phenology in situ, we focus on the phenological timing of adult flight period, with the assumption of advanced flight period mirroring an adaptive response.
The area of Finland extends over approximately 1100km latitudinal gradient across boreal and subarctic climates and is topographically relatively homogenous. The climatic isoclines roughly follow latitudes (Ahti et al., 1968), but they are currently shifting northwards due to climate change (Jylhä et al., 2010). Species range shifts as a response to climate change should thus be observable in Finland through colonization of new suitable areas north of their previous distribution (cf. Pöyry et al., 2009). To study the strategy of colonizing new suitable areas, we focus on measuring shifts in the latitudes of the northern range boundary (NRB). Although elevation can affect this response, as the distance to the next temperature isocline is shorter across altitudes than over latitudes, we do not take elevation into account in this study. The terrain in Finland is rather flat, and most of the data that we use are distributed towards the southern, less elevated parts of Finland (Figure 2), and only a minority of data come from altitudes over 500 m ( Figure S1b,c).

Systematic data for phenological analyses and population trends
To analyse phenological change in adult flight periods, we use ca. 20 and 25 years of systematically collected data from butterfly transects and moth traps, respectively. The raw data and how they were refined for the purposes of phenological analyses and population trend estimation in this study are described more in depth in Text S1. The final data for phenological analyses and population trend estimations consisted of 725,106 and 355,970 abundance observations of the 289 species, respectively.
Observational data for shift in the northern range boundary (NRB) To analyse shift in the NRB, we used observational data on the 289 species that were selected for the study through defining data for the phenological analyses (see above and Text S1). 2,474,498 observations on the F I G U R E 2 Spatial distribution of data across latitudes and elevation. Panel (a) shows the location of butterfly transects (n = 158) and moth traps (n = 93) where phenological data on the 289 species used in this study had been collected during 1999-2017 and 1993-2017 for the butterflies and moths, respectively. Panel (b) and (c) show the locations with observational data for T 1 and T 2 , respectively, with samples size per coordinate depicted by different shades of blue, on the 289 species used for calculating shifts in the northern range boundary. Also see Figure S1 for the temporal span of the data sets used in this study. The elevation data at 25-m resolution with a precision of 7 m is based on the European Digital Elevation Model (EU-DEM), version 1.1. Published on April 20, 2016 https://land.coper nicus.eu/image ry-in-situ/eu-dem/eudem-v1.1?tab=mapview | 1623 289 study species available through the Insect Database and National Butterfly Monitoring Scheme (Saarinen et al., 2003) were sourced from the Finnish Biodiversity Information Facility (FinBIF) in December 2019 (Text S2).
To measure shift in NRB over time, we compared NRB between two 5-year periods: 1992-1996 (hereafter T 1 ) and 2013-2017 (hereafter called T 2 ). The total number of presence cells for T 2 was, however, substantially higher than for T 1 due to overall increased sampling effort over time. To avoid effects caused by differences in overall observation intensity, we randomly subsampled the observations in T 2 across species to equal the number of observations in T 1 (Text S2). The temporal span of the distributional data in relation to the systematic data used to estimate phenological shift is depicted in Figure S1.

Data on life-history traits
We used data on four key life-history traits to potentially explain shifts in phenology and NRB and in population trends. The trait variables used included body size (continuous variable), overwintering stage (factor with four levels: Adult, Larvae, Pupa, Egg), specificity of larval host plant use (factor with two levels: Generalist, Specialist), and voltinism (factor with two levels: Univoltine, Multivoltine). These traits were chosen as they were linked to variation in responses to climate and environmental change in previous empirical studies on Lepidoptera (Betzholtz et al., 2013;Pöyry et al., 2009Pöyry et al., , 2017WallisDeVries, 2014) and because information on these traits is available for both moths and butterflies. See Figure S3 for further details and data sources.

Shift in phenology
We analysed phenological shift with the following linear mixed effects model: where Y is the day of year that species i was observed, a is the intercept (average day of year for observing species i), X is the year when the observation occurred (continuous, centered variable), b is the effect of Year on day of year(Y), Z is a design matrix for the site-level random effects u, and ɛ are the residuals weighted by the abundance of species i on the observation day. The model was fitted separately for each species and resulted in a slope coefficient describing the mean annual shift in adult flight timing over the study period.

Shift in NRB
To estimate shifts in NRB, we calculated the average latitude of the ten northernmost grid cells and the prevalence (proportion of occupied 10 × 10 km 2 cells for each species out of all cells with observations of any of the species) in both in T 1 and T 2 (following Thomas & Lennon, 1999). We then subtracted the NRB in T 2 from that in T 1 and similarly the change in prevalence by subtracting the prevalence in T 2 from that in T 1 . The statistical significance of an average shift in NRB across all species can be estimated by modelling the change in kilometers between the periods as a function of change in prevalence. This approach was first presented by Thomas and Lennon (1999) and has been used in similar analyses (e.g., Brommer, 2004;Mason et al., 2015;Pöyry et al., 2009). If the intercept is positive and significantly different from zero, the inference is that the species group has, overall, shifted their NRBs more towards the north than expected purely from their change in prevalence. We used this approach to obtain a linear effect estimate of shift in km as a function of change in prevalence ( Figure S4a). To obtain a corrected measure of NRB shifts per species, where the effect of prevalence change is reduced, we extracted the residuals of the model. This estimate thus describes the residual shift in NRB that is not explained by the linear effect of prevalence across the studied species. These residuals correlate strongly with the raw shift in kilometers ( Figure S4b), but are a more conservative metric as the linear effect of prevalence change across species has been removed. To ensure that a change in degrees of freedom caused by using the residuals instead of the raw estimates do not affect our conclusions, we also report the model statistics of models on the raw estimates including change in prevalence as an extra covariate for all models where we used the residual shift in NRB as a response (reported in Table 2, Tables S2 and S3). Further, to ensure that the results and conclusions were not biased by underlying differences in sampling patterns between the two periods, we conducted the same analysis with systematically collected data to confirm that the observational data does not overestimate the direction of shift in NRB (Text S3).

Population trends
We calculated the population trends from their collated annual abundance indices. This was done separately for butterflies and moths.
For the butterflies, at each site, annual indices were computed from weekly counts, following the method described in Dennis et al. (2013Dennis et al. ( , 2016 and implemented in the rbms R package (Schmucki et al., 2020). Missing week counts were derived from a Poisson generalized linear model (GLM) that included the regional flight curve as an offset (Schmucki et al., 2016). Collated annual abundance indices were then estimated with a weighted Poisson regression, accounting for site, transect length, and using the proportion of the flight curve monitored as weight. Thereby, sites with many missing counts during the flight period had lower weight than well monitored sites. For each species, we calculated the long-term trends with a linear model that we fitted on the log 10 transformed collated annual indices, starting with the year the species was first recorded until the last within the 1999-2017 period.
For the moth species, population trends were estimated using the TRIM software (Pannekoek & Van Strien, 2005), as implemented in the rtrim R package (Boogart et al., 2020). TRIM uses Poisson regression to estimate annual abundance indices while accounting for missing observations, site differences, overdispersion, and temporal autocorrelation. As a long-term trend estimate, TRIM calculates a regression through the annual indices, and this linear trend slope (on the log scale; the "additive" slope in TRIM) was used as a measure of population trend for the moth species over 1993-2016. Four species appeared in the dataset after 1993, with the first year of occurrence marking the start of the timeframe for trend calculation.

Conversion of species responses into categories
The slope coefficients from the phenology models and the residuals of NRB shift formed our main results and are hereafter jointly referred to as the responses (as in climate change response). To enable a comparison of directionality in responses, we converted the continuous results into categories (Table 1). For the sake of visualization in Figure 3, and in order to calculate the percentage of species with population trends in different directions, population trends were also assigned into similar categories: a significantly positive trend was assigned into the "Positive" group, an insignificant trend with any sign into the "Stable" group, and a significantly negative trend into the "Negative" group.

Phylogenetic generalized linear models
To test the effect of traits on the responses and population trend, the responses on each other, and the effect of responses on population trends, we applied phylogenetic generalized least squares (PGLS; Freckleton et al., 2002; also called phylogenetic generalized linear models; Symonds & Blomberg, 2014) through the pgls function as implemented in the caper R-package (Orme et al., 2018). The PGLS method is an extension to generalized least squares where the phylogenetic relationships of species are incorporated into the modelling framework via estimation of covariance in multispecies data. Related species cannot be considered independent from each other in neither their life-history traits (Freckleton et al., 2002;Ives & Zhu, 2006) nor in their responses to environmental change (Davies et al., 2013;Fei et al., 2017). Thus, any model residuals of closely related species would often be more similar than by chance, which requires modification to the estimated slopes and intercepts of the models (Revell, 2010;Symonds & Blomberg, 2014 Notes: (a) The slope coefficients from the phenological analyses and NRB shift estimates were used to convert continuous responses into categorical response groups. (b) To enable tests of our hypotheses, we further combined these phenology and NRB response groups into (1) a four-level category describing combined response groups (RG) and (2) a three-level category describing hypothesis-wise responses (Hypotheses 0, 1, or 2), both corresponding to the potential combined responses shown in Figure 1b.

| 1625
To allow controlling for phylogenetic dependence, we constructed a phylogenetic tree for the 289 species based on the hypothesis derived by Pöyry et al. (2017) for Nordic Macrolepidoptera (available in the associated data: https://doi.org/10.5061/dryad.hhmgq nkgh). We measured the phylogenetic signal in our data and found the signals to be weak but significant for shift in NRB and population trend, which confirmed that controlling for phylogenetic relatedness in subsequent analyses was appropriate to account for phylogenetic nonindependence of the species. For more details on this and description of multicollinearity checks and scaling of variables, see Text S4.

Models for hypothesis testing
To test the degree to which Lepidoptera use range shifts and/or phenology shifts as a response to climate change, the effect of the responses and their combinations on population trends, and the role that life-history traits may play for identifying species able to capitalize on the responses, we use the PGLS models described above. For the purpose of statistical analyses, these categorical variables were converted into dummy variables, and the trait body size was square-root transformed. First, to measure average response across species, we conducted intercept-only PGLS models on both F I G U R E 3 Direction of shift in NRB, phenology shift, population trend, and hypothesis-related response per species, organized by their phylogeny. There was no systematic community-wide direction in any of the three responses (phenology shift: mean = 0.02 ± 0.03 days/year; t = 0.66; p = 0.51; NRB shift (residual km): mean = −40.8 ± 40.8 km; t = −1.00; p = 0.32; population trend: mean = −0.005 ± 0.01; t = −0.32; p = 0.75). We also found no difference in the estimates between the two main taxonomic groups, butterflies, and moths (Table S2); 45.3% (33.2 + 12.1) of the studied species either shifted their NRB northwards or phenology earlier (Hypothesis 2), 39.8% neither shifted their NRB northwards nor phenology earlier (Hypothesis 0), and 14.9% of the studied species shifted both their NRB northwards and phenology earlier (Hypothesis 1). We found weak, however significant, signals of phylogenetic effects on shift in NRB and population trends (Text S4). The northwards shift was prominent in Erebidae whereas a clear shift towards the south can be seen among parts of the butterfly and Geometridae. Negative population trends were most prominent for distinct parts of Geometridae and Noctuidae. Moth icon by Carpe Diem and butterfly icon by tulpahn, both from the Noun Project https://theno unpro ject.com/ continuous responses and the population trend. For this purpose, we used the unscaled versions of the variables to allow inference directly on the measured scale (days, kilometers). The model is defined as where Y stands for the continuous dependent variable, a is the intercept (average response), and ɛ are the residuals with covariance matrix C, which is optimized based on the phylogeny.
Second, to test our main hypotheses (Figure 1), we fitted five separate models. Hypothesis 1 suggests that either phenology advance or a northwards shift of NRB would be mirrored in positive population trends. Thus, we tested the effect of direction of shifts in (a) phenology and (b) NRB on population trend. Hypothesis 2 postulates that an advance in phenology would increase the probability that species can shift their NRB northwards. Thus, we tested the effect of (c) the direction of shift in phenology on shift in NRB. To test the effect of different combinations of the responses on population trends, and which hypothesis provides a more viable strategy for species (as opposed to which strategy is the more common), we tested the effect of the (d) combined responses and (e) hypothesis-wise responses on population trend.
Finally, to test the effect of life-history traits on the responses and population trend, we fitted three PGLS models, one on each continuous and scaled response and population trend as response variables, with all four lifehistory traits as explanatory variables. We also applied PGLS models to test for a potential difference in the responses and population trends of the major taxonomic groups (moths and butterflies).
Models (a)-(e), and models on the effect of traits and taxonomic groups are structured as follows: where Y stands for the continuous dependent variable, a is the intercept (average response), X is a 289 × k − 1 dimensional design matrix indicating the independent factorial variables, b are the effects of the independent variables, and ɛ are the residuals with covariance matrix C, which is optimized based on phylogenetic signal (Symonds & Blomberg, 2014).
All data management and analyses were conducted in R studio (Data management in R version 3.5.3; estimate derivation and analyses in R version 4.0.5; R Core Team, 2021).

R E SU LT S
We found, on average, no systematic shift in any direction across the species responses ( Figure 3; the estimates of the intercept-only PGLS-models were not significantly different from zero). Among the 289 species studied, 48.1% of species expanded their NRB towards the north and 27% of species advanced their phenology (Figure 3). By contrast, 41.5% of species contracted their NRB towards the south and 35.6% delayed their phenology. We also found no difference in the estimates between the two main taxonomic groups, butterflies and moths (Table S2).
Almost half of the species (45.3%) responded according to Hypothesis 1; that is, they either shifted their NRB northwards or their phenology earlier but not both. A minority of the studied species (14.9%) responded according to Hypothesis 2; that is, they shifted both their NRB northwards and phenology earlier. Finally, 49.8% of the species showed no assumedly adaptive response; that is, they neither shifted their NRB northwards nor phenology earlier (Hypothesis 0).
More than half of the studied species (61.5%) showed positive or stable population trends, but on average, there was no systematic trend for neither positive nor negative population trends across the studied species (Figure 3). Nevertheless, population trends differed between species that responded differently in NRB and phenology shifts. Species that advanced their phenology showed more positive (although insignificant) population trends over the study period than those that delayed or did not change their phenology (Table 2, model a). Species that shifted their NRBs further north (>20 km) showed significantly stronger positive population trends compared with other species (Table 2, model b; Figure 4a). In addition, species that advanced their phenology tended to move their NRB more towards the north, but this effect was not significant (Table 2, model c). The positive effect of a northwards shift in NRB on population trends in model (b) was also mirrored in model (d), which indicates that both combined responses including a northwards shift, no matter how the species reacted phenology-wise, showed stronger positive population trends (Table 2, model d). This effect was stronger for species that also advanced their phenology. Thus, species able to utilize a combined response as postulated by Hypothesis 2 (both northwards shift of NRB and advance in phenology) showed significantly stronger population trends (Table 2, model e; Figure 4b). An ability to utilize either of the responses (as postulated by Hypothesis 1) showed, on average, lower but also positive population trend (Table 2, model e; Figure 4b). The species that were not able to utilize either of the presumed adaptive responses (Hypothesis 0) showed the lowest and on average negative population trends (Table 2, model e; Figure 4b).
None of the four life-history traits tested showed a significant connection with population trends. Overwintering stage was the only trait that had an effect on shift in phenology and on NRB ( Figure S6; Table  S3). Species overwintering as adults were more likely to advance their phenology, whereas species that overwinter as pupae tended to retreat their NRB towards the south. Due to imbalance in the number of species

| 1627
representing different host plant-use categories, we had combined species feeding on lichen and fungi into the specialist group (eight species; 2.7% of studied species: see Section "Methods"). In an additional PGLS analysis treating lichen and fungi feeders as a separate group, the eight species that feed on lichen and fungi showed slight but insignificant shift in their NRB further north (t = 1.77; p = 0.07) and had more positive population trends (t = 3.91; p < 0.001; Figure S7).

DI SC US SION
Our analysis reveals that Lepidoptera in Finland most often use only one of the commonly assumed adaptive responses to climate change, as 45% of the studied species responded by either shifting their NRB northwards or by advancing their phenology. Nearly as large a proportion (40%) did not utilize either of the two strategies. Importantly, unresponsiveness coincided with more negative population trends. In contrast, the 15% of species that responded by both shifting their NRB northwards and advancing their phenology showed, on average, the strongest positive population trends. This minority of species, able to capitalize on both responses, advanced their flight period by 3.7 days/decade and shifted their NRBs 113.1 km further north between T 1 and T 2 , on average. Although this study cannot provide evidence for de facto underlying processes giving rise to the observed patterns, the results point to potential positive effects on species abundance when combining in situ adjustments with range shifts.
We found only few trait effects on the responses and population trends (cf. Angert et al., 2011;Coulthard et al., 2019;Pöyry et al., 2009), which raises concerns related to the potential to identify vulnerable species based on their life history. However, we found that adult overwintering species tended to advance their phenology and species overwintering as pupae were less likely to shift their ranges further north, with some even contracting F I G U R E 4 Relationship between responses and population trends. Panel (a) shows the mean population trend (±SEM) for species with different combinations of responses in phenology and NRB. The pgsl model (b) indicated that a shift in NRB had a significant positive effect on the population trend (Est. = 0.65; t = 3.36; p < 0.001). Panel (b) shows the mean population trend (±SEM) for species with different hypothesis-wise responses (dark grey = Hypothesis 0-the species neither shift NRB nor phenology; dark blue = Hypothesis 1-the species shifts either NRB or phenology; light blue = Hypothesis 2-the species shifts both NRB and phenology). Species that responded according to Hypothesis 1 showed stronger population trends (Est. 0.31; t = 2.51; p < 0.05), whereas species that responded according to Hypothesis 2 showed the strongest population trends (Est. = 0.52; t = 3.08; p < 0.01) their range southwards. Earlier studies have hypothesized that species overwintering as adults are among the species mostly benefiting from increased spring temperatures, whereas species overwintering as pupae are likely to increase the number of generations produced per year (Teder, 2020;Virtanen & Neuvonen, 1999). The lack of effect of voltinism was surprising in light of the general trend towards an overall increase in voltinism in European Lepidoptera (Altermatt, 2010;Pöyry et al., 2011), which would assumedly affect both species phenology and range shifts. However, the effect of added generations may not offer actual benefits to all species, rather "tricking" some species into another generation too late in the season, which could cause a so-called lostgeneration effect (cf. Pöyry et al., 2011).
Responding by advancing phenology was, overall, relatively rarely observed among the studied species, as a striking 73% of the species did not advance their phenology, or even delayed it. This result is surprising and in contrast to several reports that are based on similar data and show strong advances in Lepidoptera phenology (e.g., Diamond et al., 2011;Roy & Sparks, 2000;Stefanescu et al., 2003). However, most of the previous studies have been conducted in temperate regions, and regional differences in abiotic conditions and rates of climate change may induce different responses (Renner & Zohner, 2018). Several other species groups in Finland have, however, been reported to advance their phenology (e.g., Helama et al., 2020, for plants;Lehikoinen et al., 2019, for birds). It is possible that other environmental factors, such as light conditions in northern latitudes (Arietta et al., 2020;Hodgson et al., 2011), the phenology or palatability of host plants (Rytteri et al., 2021), or variation in weather conditions in the early season, limit the possibility for advanced phenology of Lepidoptera. These findings are, however, in line with those by Fric T A B L E 2 Model statistics of PGLS models a-e. Variables that showed a statistically significant effect in bold Notes: To test our main hypotheses (see Figure 1) we fitted five separate models. We tested the effect of direction of shifts in (a) phenology and (b) NRB (Table  1a) on population trend (continuous). We tested the effect of (c) categorical shift in phenology (Table 1a) on shift in NRB (continuous). To ensure that a change in degrees of freedom caused by using the residuals instead of the raw estimates did not affect our conclusions, we also report the model statistics of a model on the raw estimates including change in prevalence as an extra covariate (model c-x). To test the effect of different combinations of the responses on population trends, and which hypothesis provides a more viable strategy for species (as opposed to which strategy is the more common), we tested the effect of the (d) combined response groups (Figure 1b; RGs in Table 1b) and (e) hypothesis-wise responses (Figure 1b; Hypotheses in Table 1b) on population trend (continuous).

| 1629
et al. (2020) who observed less advancement and even delays in early flight periods of butterflies towards higher latitudes.
Although there was a positive connection between phenological response and population trends, the effect was not statistically significant. Radchuk et al. (2019) found that even though phenological advance is often stated as an adaptive strategy under climate change, this may not be the rule for all species, nor is phenological advance always enough to provide fitness benefits under ongoing rapid environmental change. Climate change also introduces more variable weather conditions (Rummukainen, 2012;Vasseur et al., 2014), whereby environmental cues may become less reliable and advanced phenology may not offer the expected fitness benefits but even cause declines in readily responding species. Additionally, declining populations may not only be less able to disperse and colonize new areas (fewer individuals that emigrate) but also have a lower potential for adjusting in situ, because of loss of genetic variability (Anderson, 2016). Large declines in insect populations have recently been reported, and although these trends vary greatly between regions and taxa (Crossley et al., 2020;van Klink et al., 2020;Pilotto et al., 2020), our results also point to comprehensive population declines among Finnish Lepidoptera as 38.5% of the studied species showed negative population trends.
In contrast to advanced phenology alone, northwards shifts in NRB was associated with significantly stronger population trends. Our results also show that the studied species are more often capitalizing on range shifts than phenology shifts. This is in line with previous studies that documented strong range shifts among Lepidoptera (Kharouba et al., 2009;Mason et al., 2015;Parmesan et al., 1999;Pöyry et al., 2009) and points to range shifts perhaps being a more readily available response for many species of Lepidoptera. Simultaneously, however, only less than half (48%) of the species studied here had shifted their NRBs northwards. Habitat availability plays a crucial role when species are moving as a response to climate change (Platts et al., 2019), and other abiotic factors than rising temperatures are likely to affect the ability of species to shift their ranges (Spence & Tingley, 2020). In Finland, decrease in the area and quality of suitable habitats is known to have substantial negative effects on butterflies (Ekroos et al., 2010;Kuussaari et al., 2007;Pöyry et al., 2018). This highlights the importance of considering species dispersal in land-use planning as it is one of the main pathways through which species can adjust to ongoing changes. Halting habitat decline and fostering the persistence and even reconstruction of large and connected habitat areas can help sustain large enough populations that can both colonize new area and harbour sufficient genetic and phenotypic variation to respond, in situ, to global changes. Policies like the European Union's Biodiversity Strategy for 2030 that aims at protecting at least 30% of terrestrial and aquatic areas (European Commission, 2020) could enable more species to combine the two strategies for staying within their thermal niche.
Our study highlights that combining advanced phenology and a northwards range shift is connected with the best potential for population viability. Among boreal Lepidoptera, however, only a small proportion of species are currently able to use both responses to form a winning strategy. Similarly to earlier studies (Amano et al., 2014;Socolar et al., 2017), we find that species most commonly use either one of the strategies. Our results, however, indicate that a complementary strategy is a good but not necessarily the most optimal response as measured through population trends, but instead may be the only available option for many species. Together with the large proportion of species that were not able to utilize either of the adaptive responses, this indicates that moths and butterf lies in Finland are presently on a track towards becoming either winners or losers and that this division is likely strongly affected by habitat availability and species' abilities to make use of newly available habitat and adjust appropriately within their ranges.

AC K NOW L E DGE M E N T S
We are grateful to all the volunteers who carried out the systematic sampling of butterflies and moths. We thank Tanja Lindholm and Bess Hardwick for assistance in data management and compilation of trait information. We acknowledge Manuel Frias for advice on visual elements and Mirkka Jones, Jari Oksanen, Jukka Sirén, and Torsti Schultz for statistical advice. We thank Dr. Oliver Schweiger and one anonymous reviewer for constructive comments on the manuscript, which helped to improve the paper. We thank the Finnish Ministry of the Environment for financial support for the moth monitoring scheme.

CON F L IC T OF I N T E R E ST
The authors declare no conflict of interest.

AU T HOR S H I P
MH conceived the study, and JH, MK, JP, and MS took part in designing the study. JP, JH, MK, and RL coordinated the data collection. JP and PS provided trait data, and JH, IK, and RS calculated population trends. MH carried out data management and statistical analyses. MH drafted the manuscript while all authors contributed substantially to the manuscript and gave final approval for submission.

PE E R R EV I EW
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ele.13774.

OPE N R E SE A RC H BA DGE S
This article has earned an Open Materials badge for making publicly available the components of the research methodology needed to reproduce the reported procedure and analysis. All materials are available at: https://doi.org/10.5061/dryad.hhmgq nkgh.