Urbanization extends flight phenology and leads to local adaptation of seasonal plasticity in Lepidoptera.

Urbanization is gaining force globally, which challenges biodiversity, and it has recently also emerged as an agent of evolutionary change. Seasonal phenology and life cycle regulation are essential processes that urbanization is likely to alter through both the urban heat island effect (UHI) and artificial light at night (ALAN). However, how UHI and ALAN affect the evolution of seasonal adaptations has received little attention. Here, we test for the urban evolution of seasonal life-history plasticity, specifically changes in the photoperiodic induction of diapause in two lepidopterans, Pieris napi (Pieridae) and Chiasmia clathrata (Geometridae). We used long-term data from standardized monitoring and citizen science observation schemes to compare yearly phenological flight curves in six cities in Finland and Sweden to those of adjacent rural populations. This analysis showed for both species that flight seasons are longer and end later in most cities, suggesting a difference in the timing of diapause induction. Then, we used common garden experiments to test whether the evolution of the photoperiodic reaction norm for diapause could explain these phenological changes for a subset of these cities. These experiments demonstrated a genetic shift for both species in urban areas toward a lower daylength threshold for direct development, consistent with predictions based on the UHI but not ALAN. The correspondence of this genetic change to the results of our larger-scale observational analysis of in situ flight phenology indicates that it may be widespread. These findings suggest that seasonal life cycle regulation evolves in urban ectotherms and may contribute to ecoevolutionary dynamics in cities.

Global urban cover is projected to expand by 180 to 590% during the 21st century (1). Correspondingly, urbanization has recently emerged as a strong source of selection and a potential agent of evolutionary change, with clear evidence of the adaptive evolution of organisms to cities in some plant and animal species (2–4). A major challenge, however, of studying urban evolution has been attributing observed, phenotypic changes to genetic evolution rather than just phenotypic plasticity, and few studies have tested for the evolution of plasticity itself in urban environments (4).

Seasonal plasticity provides an excellent opportunity to test for the evolution of phenotypic plasticity in urban environments. Seasonal plasticity is a classic example of a trait that can rapidly evolve in response to anthropogenic change and relies on the ability to use cues, such as photoperiod, to predict seasonal environmental changes. The relationship between cues and environmental changes can be easily disrupted by human activities, leading to developmental traps in which organisms produce inappropriate plastic phenotypes because of cue–environment mismatches (5, 6). Research on climate change (7, 8) and range-expanding invasive species (9–11) has shown that rapid evolution has the potential to correct for mismatches between photoperiodic cues and seasonality, at least in these contexts. Nevertheless, the potential for photoperiodic reaction norms to evolve in response to urban environments has previously not been tested, despite the potential for altered urban phenology across a variety of insects (e.g., refs. 12 and 13).

Two major consequences of urbanization that could lead to the evolution of photoperiod-cued seasonal plasticity are the urban heat island effect (UHI) and artificial light at night (ALAN). The UHI—urban areas being warmer than nearby rural areas—potentially shapes natural selection and evolution in urban environments, especially for ectotherms (14, 15). For example, the UHI selects for increased heat tolerance in urban water fleas and ants (16, 17). Moreover, understanding how populations evolve in response to rapid temperature rise in cities can also provide insight into the contemporary adaptation to climate change, given that the magnitude of urban warming (averaging 1 to 5 °C) approximates current and forecasted levels of climatic warming across the globe (18, 19).

Whereas the UHI is likely to directly select on phenology by altering seasonal weather conditions, ALAN may instead disrupt the ability to predict those conditions by altering daylength independently of the season itself. This could lead to a developmental trap in which organisms continue nondiapause development into winter conditions, with severe fitness consequences (5, 6). For example, simulated ALAN substantially reduced diapause incidence in mosquito eggs and flesh fly pupae (20, 21), although urban and rural mosquito demes responded similarly (20). In moths, ALAN can disrupt larval development and pupal diapause (e.g., ref. 22), indicating that ALAN could also negatively affect fitness of day-flying Lepidoptera through impacts on earlier life stages (23). For instance, in horse-chestnut leaf miners, the number of diapausing pupae halved in continuously illuminated host plants (24).

Both the UHI and ALAN are predicted to lead to an extended flight season with adult activity later in the year, but they would have very different fitness consequences, leading to opposite evolutionary changes in the underlying reaction norms (Fig. 1). The UHI should extend the length of the growing season in urban environments, similar to the broader consequences of climate change for phenology (13, 25, 26). This would increase the fitness of genotypes that can effectively use this extra time, leading to selection for the direct development to continue later in the year at shorter photoperiods. In contrast, ALAN will lengthen the photoperiod, making species that use photoperiod as a cue mistakenly develop directly later in the year than they would under ordinary conditions (a developmental trap). This should lead to selection for the induction of diapause at longer photoperiods in urban environments, the opposite direction of the change predicted by the UHI. ALAN could also select for reduced sensitivity to photoperiod as a cue, leading to a shallower slope of the photoperiodic reaction norm. Here, using a butterfly and a moth as models, we set out to test for phenological shifts in urban areas and whether the underlying reaction norm for seasonal plasticity has evolved in the way predicted by the UHI or ALAN, as a mechanism for this phenological shift. Both the multivoltine green-veined white (Pieris napi) and the bivoltine latticed heath (Chiasmia clathrata) are widespread lepidopterans, entering pupal diapause primarily using photoperiod as a cue but also responding to temperature (27–30).

Fig. 1.

Predicted outcomes of adaptive evolution on the photoperiodic reaction norm for diapause, following selection by the UHI or ALAN. The UHI should speed up development and extend the growing season, increasing the fitness of genotypes that can continue direct development later in the season, while ALAN should extend the experienced photoperiod cue, leading to direct development too late in the season to complete development (i.e., developmental trap) and decreasing the fitness of those genotypes. Thus, if the UHI is the main factor selecting for change, CDL—daylength at which 50% of a population enters diapause—should be shorter under laboratory common garden conditions in urban than in rural populations and vice versa if ALAN is the main factor (orange versus yellow dashed vertical lines, respectively).

First, we test whether the flight season in urban environments has lengthened, as predicted by both the UHI and ALAN hypotheses, in the three largest cities of both Finland and Sweden (Fig. 2). For both species, we amassed long-term standardized data from monitoring schemes and citizen science data from observation schemes and compared the flight curves of urban and nearby rural populations for the six city regions. Then, we compared urban and rural populations from the two largest cities (Helsinki and Stockholm) to test for the genetic shifts in their photoperiodic reaction norms for diapause predicted by the UHI or ALAN. For each species, we used a common garden split-brood experiment, in which we reared offspring from wild-caught females along a controlled, photoperiodic gradient to establish among-population variation in the response to photoperiod. If the UHI was the primary factor selecting for change in these reaction norms, we expected shorter critical daylengths (CDLs) (i.e., daylength at which 50% of a population enters diapause) in urban populations, whereas if ALAN was the primary factor, we either expected longer CDLs and/or shallower photoperiodic reaction norms in urban populations (Fig. 1).

Fig. 2.

Map showing the locations of six cities (black dots) selected for phenological analyses based on citizen science data from monitoring and observational schemes. The colored dots indicate the urban (overlapping with black dots) and rural populations where C. clathrata (magenta) and P. napi (blue) were sampled for the experimental measurement of photoperiodic reaction norms of diapause induction.

Results

Monitoring and Citizen Science-Based Phenology.

To verify whether, and to what extent, the phenology of P. napi and C. clathrata has changed in urban environments, we gathered phenological data on adults of both species from various standardized monitoring and citizen science-based observation schemes. Using generalized additive models (GAMs), we fitted yearly phenological flight curves to these daily count data, separately for the urban versus rural context (0 to 10 km versus 20 to 50 km from the city center, respectively), of the three largest city regions in each of Finland and Sweden (Fig. 2). Linear mixed-effects model analysis of these flight curves showed that the flight periods of both study species were typically longer in urban than in nearby rural environments for the six city regions analyzed (Fig. 3 and ). However, some of these differences were not statistically significant at the 95% confidence level, contributing to environment by city interactions, with even an opposite difference for C. clathrata in Göteborg (Fig. 3 and ). Both species had significantly longer flight periods in urban Helsinki than in surrounding rural sites (estimated rural–urban contrast [95% CI]: C. clathrata: 21.3 [12.2, 30.4] days, z = 4.60, and P < 0.0001; P. napi: 37.1 [20.4, 53.7] days, z = 4.36, and P < 0.0001). For the Stockholm region, however, the effect was not statistically significant (C. clathrata: 8.68 [−11.9, 29.2] days, z = 0.828, and P = 0.41; P. napi: 8.02 [−6.07, 22.1] days, z = 1.12, and P = 0.26).

Fig. 3.

Field observation–based estimates of the flight period length (Top) and end (Bottom) in C. clathrata (Left) and P. napi (Right) in urban and rural areas in and around the three largest cities in both Finland and Sweden (note the difference in y-axis scaling between the species). In each panel, points are fitted values from a linear mixed-effects model (fixed effects presented in ; model-averaged fixed effects—that did not indicate an interaction between city region and environment—are used in the bottom right panel [see Statistical Analyses and for details]), and whiskers indicate 95% CIs of the fitted values, with lines connecting urban and surrounding rural areas for each of six differently colored cities. The city- and environment-specific locations from where we have additional experimental data (results illustrated in Fig. 4) are indicated with squares.

Fig. 4.

Population-specific probability of direct development in C. clathrata (Left) and P. napi (Right) in relation to experimentally manipulated daylength. Each curve is for a different population (i.e., Helsinki rural, Helsinki urban, Stockholm rural, and Stockholm urban) and is drawn on the grounds of model-averaged fixed effects (Table 1) of generalized linear mixed-effects models, explaining the data. The shaded regions around the curves are 95% CIs. The fitted regression curves are for females. For males, the curves are similar, with no significant effect of sex in C. clathrata and only a small independent effect of sex in P. napi (see Table 1). Although centered daylength was used in the statistical analyses, we relocated the regression curves to the original daylengths for an easier interpretation of the figure. The horizontal dashed line indicates the probability of 0.5 so that the fitted regression curves intersect the dashed line at the CDL. Points indicate observed proportions of directly developing individuals in each population and daylength treatment, with circles being used for overlapping points.

Linear mixed-effects model analysis of the flight curves also showed that the end of the flight period was typically later in urban than in rural zones in both species (Fig. 3 and ). In C. clathrata, later last flight days in urban than in rural zones were statistically supported (95% confidence level) in three of the city regions, while Göteborg showed again the opposite difference (Fig. 3 and ). In P. napi, last flight days were consistently later in urban zones, the urban–rural difference being similar in all city regions (Fig. 3 and ). Both species had significantly later last flight days in the urban zone of the Helsinki region (C. clathrata: 16.8 [9.88, 23.7] days, z = 4.77, P < 0.0001; P. napi: 5.06 [0.395, 9.72] days, z = 2.13, and P = 0.034), but in the Stockholm region, this was the case only for P. napi (C. clathrata: 3.47 [−12.1, 19.0] days, z = 0.436, and P = 0.66; P. napi: 4.91 [0.676, 9.15] days, z = 2.27, and P = 0.023).

Photoperiod Experiment.

For the two biggest cities, Helsinki and Stockholm, we compared photoperiodic reaction norms of diapause induction, larval development time, and pupal mass between urban and rural populations. We used photoperiod gradients that were tailored to the known CDL of diapause induction in P. napi in Stockholm (31), suggestive experimental data on C. clathrata in southern Finland (29, 32), and known phenologies of the species in the two cities. For P. napi, we had access to families from urban and rural populations for both city regions, but we did not have access to C. clathrata from rural Stockholm. We used binomial, generalized linear mixed-effects models with a logistic link function to analyze these data. In both study species, the probability of direct development increased in the order Helsinki rural < Helsinki urban < Stockholm rural < Stockholm urban across the daylength gradient studied (Table 1 and Fig. 4 and ). Consequently, CDLs were shorter in Stockholm than in Helsinki populations, in accordance with the 0.9° latitudinal difference, but were also shorter in urban than in rural populations (CDL estimated from photoperiodic response curves of the female sex presented in Fig. 4 [hours; 95% CI]: P. napi: Helsinki: urban: 18.3 [18.1 to 18.6], rural 18.5 [18.3 to 18.7], Stockholm: urban: 17.8 [17.6 to 18.0], rural: 18.1 [17.9 to 18.3], C. clathrata: Helsinki: urban: 19.9 [19.5 to 20.4], rural: 21.0 [20.5 to 21.5], and Stockholm: urban: 19.0 [18.5 to 19.4]). In C. clathrata, the photoperiodic reaction norm for direct development was steeper in the two urban populations than in the rural one (Fig. 4), as indicated by the interaction between daylength and population (Table 1), with an urban–rural CDL difference of 62 min for the Helsinki populations, corresponding to a 9-d between-population difference in the timing of CDL based on civil twilight (Fig. 4 and ). While the urban–rural difference between the two Helsinki P. napi populations was not statistically significant at the 95% confidence level (Table 1 and Fig. 4; CDL difference: 11 min, corresponding to a 4- and 2-d between-population difference in the timing of CDL based on daylength and civil twilight, respectively; ), the two Stockholm P. napi populations clearly differed (Fig. 4; model-averaged difference [95% CI]: 1.02 [0.118, 1.93], z = 2.21, P = 0.027, with the Stockholm rural population as the reference; CDL difference: 19 min, corresponding to a 6-d and 3-d between-population difference in the timing of CDL based on daylength and civil twilight, respectively; ).

Table 1.

Statistically significant (risk level 0.1) model-averaged (full average) fixed effects of generalized linear mixed-effects models (with binomial error distribution and logistic link function), explaining the probability of direct development in C. clathrata and P. napi in relation to experimentally manipulated daylength (centered daylength: centered DL) in urban and rural populations within Helsinki and Stockholm, with “Helsinki rural” and “female” coded as reference categories

Species	Model parameter	Averaged estimate	95% CI*	z value	P value
C. clathrata	Intercept†	0.0360	−0.293, 0.365	0.214	0.83
	Centered DL	0.707	0.516, 0.898	7.27	<0.0001
	Population (Helsinki urban)	1.11	0.546, 1.67	3.85	0.00012
	Population (Stockholm urban)	2.02	1.41, 2.63	6.52	<0.0001
	Centered DL × Population (Helsinki urban)	0.388	0.0470, 0.729	2.23	0.026
	Centered DL × Population (Stockholm urban)	0.319	0.00932, 0.629	2.02	0.043
P. napi	Intercept†	−1.63	−2.33, −0.930	4.57	<0.0001
	Centered DL	3.16	2.35, 3.97	7.62	<0.0001
	Population (Helsinki urban)	0.601	−0.326, 1.53	1.27	0.20
	Population (Stockholm rural)	1.26	0.300, 2.22	2.58	0.010
	Population (Stockholm urban)	2.29	1.40, 3.18	5.02	<0.0001
	Sex (male)	−0.813	−1.41, −0.22	2.68	0.0073

95% CIs derived by using the adjusted SE.

This is the predicted probability (in the scale of the linear predictor) when centered DL is 0, population is “Helsinki rural,” and sex is “female.”

Larval development time and pupal mass did not differ between urban and rural populations in either species (, including Figs. S1–S6 and Tables S1–S6).

Discussion

In our common garden experiments, we found a shift toward a lower photoperiod threshold for direct development in urban populations, indicated by consistently shorter (11 to 62 min) CDLs. This result is evidence for genetic shifts in photoperiodic reaction norms following urbanization. A shift in the same direction occurred in all three urban–rural comparisons tested, involving two different cities and two different species, demonstrating the potential to replicate our results spatially and phylogenetically. These laboratory results can also help explain the results of our monitoring and citizen science-based analysis of in situ flight phenology, which showed an extended urban flight season for both species in most cities considered. For instance, in P. napi, the laboratory-estimated CDL difference of 11 to 19 min corresponds to a 4- to 6-d later occurrence of CDL calculated from daylength in the urban sampling locations, which matches well with the on average 5 d later end of the flight period in urban environments in the in situ flight curves. Specifically, the shift toward direct development at shorter photoperiods indicates that the UHI may be driving the observed, phenological changes rather than a maladaptive response to ALAN. Although both photoperiod and temperature are known to affect the phenology and life cycle regulation of multivoltine Lepidoptera and other ectotherms (33–35), the genetic changes we found support a stronger role of urban warming than of ALAN in the evolution of seasonal plasticity in urban environments. This is also emphasized by the steeper slopes of the reaction norms in urban than in rural C. clathrata, opposite to the prediction of the hypothesized ALAN effect. Nonetheless, we cannot rule out that ALAN is also contributing to the observed phenological shifts, as implied by the diapause-averting effects of ALAN in other insects (20–22). Also, it is worth noting that ALAN has increasingly recognized biological effects on organismal fitness and behavior that could lead to the evolutionary differentiation of urban and rural populations (36, 37).

Even a moderate UHI (0.2 to 1.5 °C), as we show here for two medium-sized cities (1.3 to 1.6 M inhabitants) (), probably results in different optima for the induction of diapause. Warmer urban sites allow for faster development rates and longer periods of growth and reproduction than their cooler rural surroundings (13, 38, 39). Consequently, one expects these urban sites to allow more individuals to successfully complete an extra generation near the end of the flight season (25, 26). Nonetheless, fully demonstrating that the genetic changes we found are adaptive will require reciprocal transplant experiments that would contrast the fitness of urban versus rural populations in both settings (4, 40), similar to those recently performed for acorn ants (41). To provide further evidence of the parallel evolution of seasonal plasticity to urban environments, we can test how widespread our observation is in other species and other cities, since each city and species will likely involve independent and rapid selection events. Indeed, urban-to-urban dispersal of urban-adapted variants is unlikely in most nonmigratory Lepidoptera because of limited dispersal ranges (42), and selective pressures likely differ among cities because of variation in urbanization intensity and strong, natural latitudinal gradients for selection on photoperiodic reaction norms (11, 43). Both reciprocal transplant experiments and observations of parallel evolution by independent populations across urban environments would be considered strong evidence for adaptation by natural selection (41, 44) and would suggest that the genetic signal we observed is present in other species and other city regions (45).

Although it does not provide concrete evidence for a genetic basis for these changes, our complementary monitoring and citizen science-based analysis of in situ flight phenology suggests that similar changes have occurred in several additional cities, indicating a strong signal despite the sometimes relatively low amount of urban compared to rural observations. For instance, the overall number of urban citizen science records and corresponding abundance amounted to only 52 and 13% of the rural records and to 32 and 7% of the rural abundance for P. napi and C. clathrata, respectively. While these lower urban amounts fit the generally adverse impact of urbanization on butterfly and moth abundance (46, 47), it also implies that, with a more equal number of observations in both urban and rural settings, we could expect an even stronger effect of extended flight periods in urban settings. This analysis shows typically longer flight periods for both study species that end later in the season in urban areas versus their rural surroundings, suggestive of a partial extra generation. Recent global warming has led to similar responses of prolonged flight duration and/or higher investment in late or extra generations for a significant proportion of European Lepidoptera species (25, 26). Consequently, one can expect that the UHI may further increase investment into late generations and increased voltinism in Lepidoptera. From an evolutionary perspective, it is this increased voltinism that is the most important response, since an extra generation may accelerate population growth and facilitate further adaptation (25). Nonetheless, producing an extra generation at high latitudes also entails risks; entire cohorts may perish if winter arrives before individuals reach a stage that allows successful overwintering (26). In general, urbanization can interact with regional temperature and precipitation, leading to different phenology patterns in cool versus warm and in wet versus dry regions (48, 49).

The six city regions we studied included Helsinki and Stockholm (1.3 to 1.6 M inhabitants), for which we showed experimental evidence for a genetic change following urbanization, but also four smaller cities (0.3 to 0.6 M inhabitants), for which observed responses may or may not have a genetic basis. As we made sure to only include sites so that urban–rural distances were short (10 to 40 km), our evidence for differences appears remarkable, especially for P. napi given its considerably higher mobility, and hence higher gene flow, than C. clathrata (42, 50). However, heritable divergence in the face of ongoing gene flow is achievable when a strong selection pressure leads to divergent selection and local adaptation to distinct selective optima, and indeed, divergence between populations in the face of gene flow is considered strong evidence for local adaptation (51, 52). Since the diapause trait is essential for winter survival in temperate and boreal insects, it is under strong selection, which could counteract such gene flow (53, 54). For instance, the swift expansion of the globally invasive Asian tiger mosquito across the United States was accompanied by the rapid evolution of diapause and CDL along a latitudinal cline, matching that of its native range after only 23 y (11).

Seasonal plasticity is a classic example of a trait that can rapidly evolve because of its importance for ensuring proper seasonal timing (7), while urbanization is a new source of a mixture of anthropogenic changes whose evolutionary implications we are just beginning to appreciate (3). Here, we provide evidence for the evolution of the reaction norm for diapause induction in urban environments, using two species of Lepidoptera. Lepidoptera are of fundamental importance to ecosystem functioning given their role as pollinators, herbivores, prey, and hosts for parasitoids (23, 55), so changes in their ecoevolutionary dynamics can potentially lead to disruption across trophic levels (56, 57). Understanding how Lepidoptera can adapt to the altered seasonality found in urban environments is thus critical to understanding the resilience of their roles in urban ecosystems. With increasing attention being paid to the impacts of climate change on phenology and whether evolution can correct for it (5, 58), we show here that these same questions need to be asked for urbanization as well. Our results in that regard may also apply to the wide range of other species globally that need to maintain appropriate phenologies in the face of the UHI and ALAN, and whose contemporary evolution may hence be shaped by urban seasonal change (2).

Materials and Methods

To determine how urbanization has affected the phenology of P. napi and C. clathrata, we amassed standardized monitoring data and citizen science-based observational data on the abundance and phenology of both species. For the observational data, we included data from 2005 onwards, while for the monitoring data we included all years available. We gathered the available phenological data on adults of both species from various Lepidoptera-monitoring schemes [Finnish Moth Monitoring Scheme (“Nocturna”): 1993 to 2020 (C. clathrata); Swedish Butterfly Monitoring Scheme: 2010 to 2020 (P. napi); Finnish Butterfly Monitoring Scheme: 1999 to 2020 (P. napi and C. clathrata), additionally including data (P. napi and C. clathrata) collected with the same sampling protocol for a separate ecological study in Helsinki (59)] and from two national observation portals (http://www.artportalen.se and http://www.laji.fi: 2005 to 2020 [P. napi and C. clathrata]). Cleaning of the portal data involved excluding observations referring to nonadult stages, treating the few observations that lacked count data as single individuals, and excluding a small minority of observations with counts >35, as higher abundances appeared to be estimates rather than real counts. Observations with spatial accuracy >1 km were not retained, while observations linked to grid coordinates were attributed to the respective grid centers. While the observation portals collate nonsystematic observational records, the monitoring schemes collect data systematically based on diurnal “Pollard walks” along fixed transects or based on nocturnal light sampling at fixed sites (e.g., refs. 60 and 61).

We used data associated with the three biggest city regions in each country (Finland: Helsinki, Tampere, and Turku; Sweden: Stockholm, Göteborg, and Malmö) (Fig. 2). For each of these six regions, we attributed sites as urban or rural if they were located within a 0- to 10-km versus 20- to 50-km radius from the historical city center, respectively (62). The 10- to 20-km radius buffer between urban and rural sites allows for clearer contrast between the urban and rural settings.

For both species, and for each day with data, abundance was summed per site. The observational, opportunistic citizen science data were summed within the abovementioned 0- to 10-km and 20- to 50-km radius areas, which were each treated as a single, separate site for analyses, while standardized, monitored transects and light trap locations were treated as additional separate sites. This means that the resulting flight curves are in most cases to a larger degree informed by phenology data from standardized monitoring (i.e., often several sites) than by phenology data from opportunistic citizen science observations (i.e., one site per city per environmental [urban/rural] type). Note that while all abundance data refer to exact dates, the Finnish moth–monitoring scheme data are abundances nonetheless lumped over typically 1 wk. Here, dates refer to the last day of such sampling weeks. Using the “rbms” package (63) in R 4.0.2 (64), we then fitted yearly flight curves on these (daily) count data, separately for the urban versus rural context of each of the six city regions. These flight curve computations are GAMs spline fitted on counts collected across one or multiple sites and standardized so that the area under each flight curve equals one. We structured site-specific time series with the first of April and 31st of October as the first possible start date and last possible end date of the flight season, respectively, and defined the time series daily for increased temporal resolution of the flight curves. R codes are provided in .

We adopted a two-step data retention criterion. First, in order to retain ideally a relatively large number of sites or transects jointly informing the GAMs, we imposed a minimum number of visits to a specific site or transect of only three, with minimally two occurrences of the species. Second, from the resulting flight curves (n = 613), we excluded curves (n = 68) with a flight length <30 d and/or with a rectangular shape, indicative of poor coverage of the total flight period. Thus, a total of 545 standardized flight curves were retained, for which we then calculated 1) the length of the flight period as the number of days with fitted occurrence >0.1% of the total flight curve and 2) the end of the flight period as the last day with fitted occurrence >0.1% of the total flight curve. To check the sensitivity of the last flight day estimates to the used threshold, we repeated 2) with the condition of >0.05%. Finally, to make sure that our estimates of the lengths and ends of flight periods are rather conservative than overconfident, we plotted the last day of the flight period with fitted occurrences >0.1% on those of >0.05%. This indicated only six and one flight curves of C. clathrata and P. napi, respectively, strongly disagreed in last flight day estimates between the two threshold values. Excluding these seven uncertain points, the estimates derived with the two alternative thresholds were strongly and positively correlated (Pearson correlation [95% CI]: C. clathrata: r = 0.953 [0.940, 0.996]; P. napi: r = 0.944 [0.929, 0.955]). Because of the uncertainty in the last flight day in the cases where the two estimates strongly disagreed, we removed each of these seven data points from further analyses, leaving 538 flight curves for final analyses (N = 246; N = 292).

Common Garden Photoperiod Manipulation Experiment.

P. napi females were caught in urban areas of Stockholm and Helsinki [9 at Stockholm University campus (59.362° N 18.058° E), 31 May to 2 June 2020; 15 at three sites within Helsinki (Töölö [60.178° N, 24.931° E], Vermo [60.212° N, 24.846° E] and Malmi [60.246° N, 25.045° E], May 26, 27, 29, and 30, 2020) and at two rural counterparts at similar latitudes (seven at two sites at Rindö [59.397° N, 18.401° E], May 30 to 31; 16 at four sites within Raasepori: Öby Östervik [59.930° N, 23.186° E], Öby [59.925° N, 23.165° E], Höglandet [59.975° N, 23.267° E], and Finby gränd [59.987° N, 23.240° E], May 28, 2020) (Fig. 2). All females were taken or shipped to Stockholm University for egg laying. There, females were kept at conditions favorable for oviposition in 1-L translucent plastic cups with Alliaria petiolata as a host plant and Taraxacum officinale as a nectar source. Between 24 and 48 h after egg laying, eggs were moved to a room with controlled conditions of 17 °C and 12 h light:12 h dark (12L:12D) photoperiod. Within 48 h of hatching, 25 caterpillars from each of up to 10 families per population were evenly divided among five experimental photoperiod treatments: 17L:7D, 17.5L:6.5D, 18L:6D, 18.5L:5.5D, and 19L:5D. Only larvae from females that laid at least 25 eggs were used, resulting in initial per treatment sample sizes of 40 (8 families) for Urban Stockholm, 30 (6 families) for Urban Helsinki, 35 (7 families) for Rural Stockholm, and 50 (10 families) for Rural Helsinki. Caterpillars were all raised in growth chambers (Termaks KBP6395-L/KB8400-L) at 23 °C, separately according to photoperiod treatments. Caterpillars from the same family shared a 1-L cup, in which they were fed ad libitum with wild-collected A. petiolata. Caterpillars were checked daily for pupation; the date of which was recorded. Around 2 d after pupation, pupae were weighed (Precisa Instruments XB 120A; precision 0.1 mg), sexed based on sex-specific genital scars, placed in individual cups, and moved to a separate room kept at 23 °C and 22L:2D. Pupae were monitored daily for eclosion, and all adults emerging within 21 d were labeled as undergoing direct development.

C. clathrata females were similarly caught within Helsinki (Malmi [60.246° N, 25.045° E]), Stockholm (Bromma [59.362° N, 17.927° E]), and at a Finnish rural site (Somero, Häntälä [60.576° N, 23.351° E]) in 2019 (Fig. 2). As we only managed to collect a few males and no females from rural Stockholm, our comparisons in that species focus on the Finnish populations while still retaining the urban Stockholm population for context. All females were taken to the University of Oulu, where they were kept at 21 °C under 12L:12D conditions for egg laying in individual pots with access to sugar water. Newly hatched larvae were moved individually to cups (0.25-L translucent plastic). Each larva had ad libitum access to fresh shoots of the natural host plant Lathyrus pratensis and was reared until pupation at 21 °C under 12L:12D. A layer of garden peat was added to the cups for pupation when the larvae reached the final instar. Around 5 d after pupation, the pupae were excavated from the peat and placed in a cup with Sphagnum spp. moss. After ca. 10 mo of overwintering in a refrigerator room at 5 °C, the pupae were taken to 20 °C and 18L:6D. Eclosed females were placed with a nonsibling male from the same population in cups under 21 °C and 18L:6D with access to sugar water and allowed to mate and lay eggs. Up to 35 freshly eclosed F2 larvae from each of up to 14 females per population were singly placed in cups. Initial sample sizes per population were 420 (14 families) for Rural Helsinki, 305 (10 families) for Urban Helsinki, and 395 (13 families) for Urban Stockholm. The offspring of each female was evenly divided among five climate rooms (Arctest Oy) at 21 °C with the following photoperiod regimes: 18L:6D, 19.5L:4.5D, 21L:3D, 22.5L:1.5D, and 24L:0D. Date and nearest hour were noted when these larvae were placed in cups, in which they were reared similarly to the previous generation. Final instar larvae were checked at least twice a day, and the hour and date when the larvae burrowed into the peat for pupation were recorded. Around 5 d after burrowing, pupae were excavated from the peat, weighed (Mettler Toledo MT5; precision: 0.01mg), and sexed based on sex-specific genital scars. After pupation, each individual was kept under the same conditions as those experienced during larval development, at least for the first 5 d, and monitored daily for eclosion for a minimum of 2 wk. All individuals that did not eclose during that time were considered to be in diapause. At the point when the monitoring of pupae was ended, no adults had eclosed for 3 d.

Our common garden split-brood experimental approach, hence, allowed us to test to what degree photoperiodic reaction norms of 1) diapause induction, 2) body size (i.e., pupal mass), and 3) larval development time differ between urban and rural populations associated with two replicate cities.

Statistical Analyses.

For all analyses, we used the R 4.0.2 software (64). Throughout we used the multimodel inference approach (e.g., ref. 65) to avoid problems in model selection and to reduce the risk of overconfident inferences. This was done by averaging over a set of models, including the “global model” (i.e., the most complex model considered) and all simpler models that follow the principle of hierarchy. R codes are provided in .

We analyzed the length of the flight period in both C. clathrata and P. napi by using linear mixed-effects models fitted with the function “glmmTMB” (66) and by setting a Gaussian error distribution and an identity link function. For both species, the fixed effects of the global model included the environment (urban/rural), city identity, and the interaction between them. Random effects included random intercepts for city–year combinations. We initially used numbers of observation sites as weights in the global models for both species, but this resulted in poor model goodness of fit, as indicated by diagnostic tools from package “DHARMa” (67) in P. napi. Removing these weights considerably improved model goodness of fit, so we used unweighted models for inferences in this species. In C. clathrata, the weighted model fitted the data well and was used for inferences. We derived the set of all meaningful, simpler models from the global models for both species by using the function “dredge” and derived model comparison metrics with the function “model.avg” [package “MuMIn”; (68)]. The global model was superior to any other model in both species (C. clathrata: difference in sample size-corrected Akaike information criterion [ΔAICc] = 21.1 between the two top-ranked models and global model’s Akaike weight = 1; P. napi: ΔAICc = 6.53 between the two top-ranked models and global model’s Akaike weight = 0.96). Thus, we based our inferences solely on the global models. Finally, we derived 95% CIs for the model predictions for each city–environment combination by using the function “predict.glmmTMB” (66).

We analyzed variation in the last flight day with a similar procedure, as explained in the previous paragraph. In these analyses too, weighting the observations with the number of observation sites resulted in a good model fit in C. clathrata, whereas in P. napi we used unweighted models for inferences, as they had a better fit. In C. clathrata, the global model was superior to any other model (ΔAICc = 19.0 between the two top-ranked models and global model’s Akaike weight = 1), so we based inferences solely on this global model. In P. napi, there were three models with a nonnegligible explanatory power (), so we used model averaging for inferences in this species. We used the function model.avg (68) for averaging over the set of models (including the global model) to obtain model-averaged coefficients (full average; i.e., a coefficient was assumed to be 0 in a model in which it was not included to avoid exaggerating statistical significance). The set of averaged models are shown in .

We analyzed the probability of direct development in both C. clathrata and P. napi by using generalized linear mixed-effects models [function glmmTMB; (66)] with a binomial error distribution and a logistic link function fitted with the maximum likelihood method. We included centered daylength (i.e., deviation from the middaylength), population, sex, and all interactions among them as fixed effects and family as a random effect (random family-specific intercepts) in the global model for both species. In P. napi’s global model, we also included rearing cup identity as a random effect (random intercepts) to consider potential phenotypic similarity due to a shared environment among individuals reared in the same cup. For both species, we also tested whether random family-specific slopes in relation to daylength would be needed in the global model. While the model with both random family-specific intercepts and slopes did not converge, a model with only random family-specific slopes converged and indicated negligible, family-level variation in random slopes for C. clathrata. Moreover, as AIC comparison suggested that the model with only random intercepts was better (ΔAIC = 0.75), we did not include random family-specific slopes in the global model for C. clathrata. For P. napi, AIC comparison indicated that the model with only family-specific random intercepts was better (ΔAIC = 2.1) than a model with both family-specific random intercepts and slopes in relation to daylength. Hence, for P. napi too, we included only random intercepts in the global model. We assessed the goodness of fit of the global models with the diagnostic tools available in package DHARMa (67). These tools indicated a good model fit.

Model averaging was performed as explained in Monitoring and Citizen Science-Based Phenology. The sets of averaged models are shown in ( for C. clathrata and P. napi, respectively). Based on the model-averaged coefficients, we derived 95% CIs for the population-specific regressions in relation to daylength. We used the function “predict.averaging” (68) and normal distribution approximation to derive CIs in the scale of the linear predictor, in which normal distribution approximation is valid. Next, we back transformed the fitted regressions and their CIs to the probability scale by using the inverse of the link function.