Phenotypic plasticity is aligned with phenological adaptation on both micro- and macroevolutionary timescales

In seasonally-variable environments, phenotypic plasticity in phenology may be critical for adaptation to fluctuating environmental (temperature) conditions. Using an 18-generation longitudinal dataset from natural insect (damselfly) populations, we show that phenology has strongly advanced. Individual fitness data suggest this is likely an adaptive response towards a thermally-dependent fitness optimum. A laboratory experiment revealed that developmental plasticity qualitatively matches environmental-dependence of selection, partially explaining observed phenological advance. Expanding our analysis to the macroevolutionary level, we use a database of over 1-million occurrence records and spatiotemporally-matched temperature data from 49 Swedish Odonate species to infer macroevolutionary dynamics of phenology. Phenological plasticity was more closely aligned with adaptation for species that have recently colonized northern latitudes, with more phenological mismatch at lower latitudes. Our results show that phenological plasticity plays a key role in microevolutionary dynamics within a single species, and such plasticity may have facilitated post-Pleistocene range expansion in this insect clade.


Introduction 64
In many organism, life is characterized by dramatic life-history transitions between discrete 65 stages that correspond to the unique demands of resource acquisition versus reproduction. For Here we combine three independent, large-scale datasets from the insect order Odonata 127 (damselflies and dragonflies). Odonates are semi-aquatic insects with a phylogenetically 128 conserved complex life cycle characterized by metamorphosis from aquatic nymph (henceforth, 129 'larvae') to aerial reproductive adult (Stoks & Córdoba-Aguilar 2012). Our aim was to address 130 the above two related questions on phenological plasticity's role in adaptation at short and deep 131 timescales. 132

Long term study of Ischnura elegans 134
We analyzed data on individual phenology and fitness in a set of populations (approximately 16) 135 of the damselfly Ischnura elegans in province of Skåne, southern Sweden, for 18 years (2000-136 2017). Each year, these populations were surveyed in a standardized approach during the 137 reproductive season of I. elegans. All captured individuals were sexed and their copulation status 138 (mated or not) was recorded. In this study, we used data on one major male fitness component:  Maternal egg clutches (full sib clutches) were obtained from females caught in copula in the 146 field in the general population survey described above. Individual eggs and larvae from a total of 147 101 full-sib clutches were reared in a split-clutch GxE design across both environments. 148 149

Occurrence records of Swedish Odonata 150
We obtained public, spatiotemporally referenced observational records of adult Swedish 151 odonates identified to species level from the Swedish Species Observation System (Artportalen: 152 https://www.artportalen.se/). We used data from a ten year period from 2006-2015, the range of 153 years for which we had accurate temperature data and for which substantial observational 154 records were available in this public database. We obtained spatiotemporally matched estimates 155 of spring mean surface temperature from the publicly accessible NOAA Global Forecast System 156 Analysis (GFS ANL) database (https://www.ncdc.noaa.gov/data-access/model-data/model-157 datasets/global-forcast-system-gfs). We obtained the mean surface temperature from 158 corresponding lat-long grid cell in the GFS ANL database for two time windows: the first 120 159 days of the year and the month of April. The first time window corresponds to a substantial 160 period of pre-metamorphosis larval development and covers the last instar for all the species in 161 our database, and so represents a long-term average temperature environment experienced during 162 a key stage of development. The second time window represents a shorter time slice from late 163 spring, just prior to the flight period of many early season species. Finally, we note that the scale 164 and purpose of our analysis precludes use of a "sliding window" approach; Our interest is to 165 compare developmental-temperature reaction norm variation across species (and so any such 166 comparison requires use of similar or the same environmental cues), rather than to identify a 167 specific time slice of spring temperature variation that is most important for a single species. 168 169

Methods -Statistical Analysis 170
Long term field study of Ischnura elegans 171 To infer selection on phenology, we modify the Lande-Arnold ( where and 2 are estimated across time points within binned latitude ( G>HI(EFG) ). 209 This bivariate model was fit separately for each species in our dataset. We then used the 210 corresponding estimates of B and b in phylogenetic mixed models (Hadfield & Nakagawa 2010) 211 to explore the relationship between species mean latitude, species mean phenology, and the 212 strength and form of inferred phenological plasticity. Our models accommodated uncertainty in 213 our estimates of B and b. 214 215

Results 216
We found evidence of significant phenological advance during the 18 year study period of I. 217 elegans in the field (F1, 47526 = 1231.38, P <0.0001, Figure 2B). Although there was significant 218 among-year heterogeneity, on average the emergence date has advanced by .58 days (SE = 219 0.017) per year, corresponding to a predicted total advance of 10.44 days across the entire 18-220 year period. 221 Using data on male mating success (N = 21,384 individuals for which we also had 222 annual temperature data, consisting of 5,487 mated males and 15,897 non-mated males), we find 223 evidence of statistically significant but annually variable (Figure 3B Celsius of spring warming ( Figure 3C). We also found evidence of negative B using an 238 independent dataset of 1695 observations from the Artportalen database when we inferred B 239 indirectly from interspecific data on phenology, although this indirect estimate was significantly 240 shallower than direct estimate we obtained from this long-term population study (Supplemental 241 material, Figure S2). 242 We obtained individual and family emergence time records from 521 individuals from 243 101 full-sib genetic families across two temperature treatments in our laboratory experiment. 244 We found significant acceleration of development in response to our experimental manipulation 245 of temperature (F1, 367 = 65.61, P <0.0001, Figure 3A (although the magnitude of these estimates differs substantially; Figure S2). 255 We analyzed paired temperature and phenology data from 1,027,952 records from 49 256 species of Swedish Odonates (see Table S1). We explored two temperature windows as 257 predictors of phenological variation: the mean temperature of the first 120 days of the year, and 258 the mean April temperature. We found substantial support for the 120 day window, as these 259 models consistently fit better than the April mean models for nearly all species ( Figure S3). We 260 then kept the 120 day average models for downstream analysis. 261 Reaction norm estimates varied across species, although most estimates where within the 262 range of ± five days per °C of spring temperature variation ( Figure 4A, Table S1). Notably, there 263 were a mix of both positive and negative reaction norm estimates and many species that showed 264 positive reaction norm slopes (Fig. 4). These analyses show that temperature during the last 265 instar larval stage can either advance or delay emergence, depending on species, and some 266 species were relatively insensitive to spring temperature (Fig. 4A). On average, the estimated 267 thermal reaction norms explained 1.4% of the all variation in phenology (Table S1). In 268 comparison to our estimates of b, estimates of B exhibited substantial sampling variance (Table  269   S1). 270 We found evidence of a significant phylogenetic signal in the phenological reaction norm 271   Figure 5C). When |B-b|/|B| is less 285 than unity plasticity can be considered adaptive; when above unity plasticity can be described as Here we combine three independent datasets, together representing well over 1 million 297 individual insect records, to test two related questions on the role of phenotypic plasticity in the 298 evolution of phenology. First, we demonstrate pronounced recent phenological advance in a 299 single damselfly species, Ischnura elegans, using a long term study of individuals in the wild. 300 We find that thermal plasticity in developmental time in the laboratory is substantial, 301 corresponding in sign and to some degree in magnitude with the relationship between 302 temperature and the optimum phenology of the same species estimated in the field. This suggests 303 an important role for thermal plasticity in adaptation to seasonal environmental variability, and 304 also suggests plasticity could contribute to advancing phenology that corresponds to recent 305 anthropogenically-caused climate change in Sweden (SMHI 2015). At the macroevolutionary 306 scale, we used a large interspecific database of individual emergence records of 49 species of 307 odonates spanning the extensive latitudinal gradient of Sweden. Combining these occurrence 308 records with spatiotemporally-matched estimates of spring mean temperature, we found that 309 species-specific mean phenological reaction norms are more positive and more closely aligned 310 with indirect inferences of the environmental dependence of the optimum for northern than for 311 southern species. We also find that plasticity contributes more to variation in phenology at 312 Northern latitudes. Our combined intra-and interspecific results indicate that temperature-313 dependent plasticity in phenology underlies both microevolutionary adaptation to contemporary 314 climate change as well as the macroevolutionary expansion into harsh high-latitude climatic 315 conditions that took place when species colonized these areas after the last Ice Age in 316 Fennoscandia, about 10,000 years before present. 317 We found a phenological advance in a single species (Ischnura elegans) of 0.58 318 days/year, almost six days per decade or 10.44 days over the entire 18-year study period ( Fig. 1). 319 This is a dramatic phenological advance during the first two decades of the century in Sweden. 320 In fact, this is almost twice as large advancement of phenology as was documented previously Britain, the new results presented here suggest temperature-dependent phenological advance in 324 odonates may be more extreme than previously appreciated. 325 We find some correspondence between point estimates of the environmental dependence 326 of selection in Ischnura elegans and the population mean reaction norm, using three independent 327 datasets and approaches. Our direct estimate of B using data on individual fitness from the long-328 term study on I. elegans and our indirect estimates using a public database of individual 329 occurrence records both suggest that the optimum emergence time advances substantially for 330 every centigrade degree of spring warming in I. elegans. Our direct estimate of the phenological 331 reaction norm matches with our direct estimate of the environmental dependence of selection, 332 and the correspondence is also found using a bivariate mixed model approach to indirectly 333 estimate these parameters (see Supplemental Material). Although absolute values of B and b 334 differ substantially across these two approaches (direct vs indirect), the qualitative conclusion -335 that phenotypic plasticity can explain the close match between local optima set by environmental 336 variation -is similar across these different approaches and datasets. One the one hand, this 337 correspondence in qualitative interpretation (alignment of B and b) is encouraging and suggests 338 that our space-for-time substitution approach to estimating environmental dependence of optimal 339 phenologies may be robust to inevitable violation of the assumptions (e.g., local adaptation and 340  2001). Although our results and analyses do indeed suggest that selection can be strong, we also 376 find little evidence of persistent directional selection of the magnitude required to produce such a 377 dramatic response. In fact, across all years, net selection on phenology is actually slightly 378 positive rather than negative as would be expected if earlier phenology was consistently favored 379 (Fig. S1). Therefore, it is unlikely that genetic evolution alone can entirely explain the strong 380 phenological advance of more than 10 days across the entire 18-year period, given the substantial 381 among-year variation (Figure 2) of I. elegans in southern Sweden. Instead, our analyses suggest 382 that phenotypic plasticity may have an important role in driving both inter-annual variation in 383 phenology and the progressively earlier flight seasons that we have documented here. These data 384 suggest that I. elegans can track changing phenological optima through plasticity, which will in 385 effect reduce or counteract directional selection for earlier phenology, a conclusion that was also 386 reached in a recent study of 21 bird and mammal species with long-term data (Villemereuil et al. Across all our interspecific dataset, all but one species in our comparative analysis over-winter as 410 either aquatic larvae or eggs, and so are exposed to spring temperatures at a similar 411 developmental stage. The only exception to this rule is the winter damselfly Sympecma fusca, 412 which overwinters as an adult in vegetation refugia. Consistent with this, S. fusca exhibited one 413 of the most extreme negative estimate of B (Table S1) plasticity in our study is not necessarily adaptive. Interestingly, there was some degree of mis-449 alignment between b and B in the southern species I. elegans (Fig. 3C), which is broadly 450 consistent with the greater mis-alignment in southern, compared to northern species at the 451 macroevolutionary level (Fig. 5C). More generally, plasticity can be adaptive and selected   in this system. We also note that female multiple mating is rare due to mate guarding in copula.  and these data are also subsequently exported to the Global Biodiversity Information Facility 957 (GBIF: https://www.gbif.org/). We used data from a ten year period from 2006-2015, the range 958 of years for which we had accurate temperature data (see below) and for which substantial 959 observational records were available in this public database. We excluded observations of 960 larvae, as well as for rare taxa for which records were unavailable in some years in this range, as 961 well as taxa that were not included in our phylogeny (see below). Our final dataset consists of 962 1,027,952 records from 49 species, which is the majority of odonate species (75 %) that have so 963 far been recorded in Sweden (currently a total of 65; 43 dragonflies (Anisoptera) and 22 964 damselflies (Zygoptera)). 965 We obtained spatiotemporally matched estimates of spring mean surface temperature 966 from the publicly accessible NOAA Global Forecast System Analysis (GFS ANL) database 967 (https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs). 968 The NOAA GFS ANL provides global weather data estimated four times per day gridded at .5° 969 resolution (corresponding to approximately 25km in Sweden), using a system of satellites and 970 globally-distributed weather stations. Comparison of these surface temperature estimates to 971 direct temperature records from a weather station in Malmö, Sweden, obtained from the Swedish 972 Meteorological and Hydrological Institute (https://www.smhi.se/en), indicates the GFS ANL 973 data are an accurate representation of surface temperature (r > 0.99 between temperature 974 records). For each record in our dataset of Swedish odonates, we obtained the mean surface 975 temperature (GFS ANL database variable: tmp2m) from corresponding lat-long grid cell in the 976 GFS ANL database for two time windows: the first 120 days of the year and the month of April. 977 The first time window corresponds to a substantial period of pre-metamorphosis larval 978 development and covers the last instar for all the species in our database, and so represents a 979 long-term average temperature environment experienced during a key stage of development. 980 The second time window represents a shorter time slice from late spring, just prior to the flight 981 period of many early season species. We focused on mean temperature, because daily 982 fluctuations in air temperature are unlikely to correlate strongly with water temperature, although 983 average temperatures are a strong predictor of average water temperatures (Livingstone & Lotter 984 1998). We note that surface air temperatures are strong a strong predictor of water temperatures 985 in freshwater bodies, particularly in the upper 1m of the water column (Livingstone & Lotter 986 1998) where odonate larvae are most abundant. Finally, we note that the scale and purpose of our 987 analysis precludes use of a "sliding window" approach; Our interest is to compare 988 developmental-temperature reaction norm variation across species (and so any such comparison 989 requires use of similar or the same environmental cues), rather than to identify a specific time 990 slice of spring temperature variation that is most important for a single species. 991 We obtained phylogenic information for all 49 species in our dataset, using a recently 992 published large dated molecular phylogeny of the odonates (Waller & Svensson 2017 ). This 993 θ 7 = − + − 2 Table S1. Summary statistics for 49 species of Swedish Odonata. B is the indirect estimate of the environmental dependence of the optimum, b is the indirect estimate of the temperature-phenology reaction norm. R2 is the variance in phenology explained by plasticity. Parameter values are posterior modes; standard errors are the square root of the posterior variance See text for further details.