Loss of leaf-out and flowering synchrony under global warming

The temporal overlap of phenological stages, phenological synchrony, crucially influences ecosystem functioning. For flowering, among-individual synchrony influences gene flow. For leaf-out, it affects interactions with herbivores and competing plants. If individuals differ in their reaction to the ongoing change in global climate, this should affect population-level synchrony. Here, we use climate-manipulation experiments, Pan-European long-term (>15 years) observations, and common garden monitoring data on up to 72 woody and herbaceous species to study the effects of increasing temperatures on the extent of within-population leaf-out and flowering synchrony. Warmer temperatures reduce in situ leaf-out and flowering synchrony by up to 55%, and experiments on European beech provide a mechanism for how individual genetic differences may explain this finding. The rapid loss of reproductive and vegetative synchrony in European plants predicts changes in their gene flow and trophic interactions, but community-wide consequences remain largely unknown.

The structure and functioning of ecosystems crucially depends on the timing of annually repeated length sensitivity differences among individuals (as documented for F. sylvatica; Fig. 3b): Under 167 cold winter conditions, days are already long when spring warming occurs, reducing the effect of 168 a tree's day length sensitivity on its leaf-out time, whereas with early spring warming, days are 169 still short, preventing day-length sensitive trees from budburst. In natural populations, leaf-out 170 advancement in day length-sensitive individuals, but not in day length-insensitive individuals, 171 will thus increase the period of leaf-out under short day conditions. Both the experimental and the 172 PEP in situ data confirm this idea, showing that (i) phenological variation among individuals 173 strongly decreases under short day conditions (Figs. 2b and 3b) and (ii) genetic differences in 174 day-length requirements are the single most important factor explaining variation in budburst 175 times (Fig. 4c,d). 176 This insight explains why, especially in Fagus sylvatica, in which day length has the 177 most pronounced effect on spring phenology (10,11), LOS is strongly affected by preseason 178 temperatures (Fig. 1c). By contrast, in day-length insensitive species, such as silver birch Betula 179 pendula and Norway spruce Picea abies (11), preseason warming has a smaller (but still 180 significant) effect on LOS, suggesting that heritable differences in day-length sensitivity are a 181 major driver of within-population phenological variation. In our common garden data, the 182 standard deviation of inter-individual leaf-out times increased by 0.09 ± 0.02 (mean ± CI) days 183 per decrease in one chilling day, and the standard deviation of inter-individual forcing 184 requirements increased by 0.23 ± 0.06 degree-days per decrease in one chilling day (lower panel 185 What biological consequences can be expected from less synchronized leaf-out and 188 flowering of the individuals of a species? With regard to vegetative development, precocious leaf unfolding under warm springs increases the risk of late frost damage (16-18), but 190 also potential carbon gain due to earlier photosynthetic activity (19). This risk-return trade-off 191 will affect selection on suitable genotypes under future conditions, and the increasing spread of 192 leaf-out should increase the selective importance of spring phenology. Whether opportunistic 193 phenological strategies (relying on temperature as the main trigger) or conservative strategies 194 (relying on day length and/or winter chilling as a buffer against highly variable spring 195 temperatures) will be favored in the future will be region-specific, depending on the relative 196 advancement rates of spring warming and late frost events. In continental regions, where the 197 advent of spring is relatively invariable (low late frost risk), phenological strategies reliant on 198 temperature should be favored (20). 199 With regard to flowering, decreased synchrony among individuals, as already strongly 200 evident in Alnus glutinosa (Fig. 1f), should lead to reduced inter-individual pollen transfer. 201 Strong divergence in flowering times among individuals also might lead to assortative mating 202 (depending on incompatibility systems), possibly promoting local adaptation (21-23) and should 203 act as a buffer against climate change-induced phenological mismatch between plants and leaf-204 feeding or pollen-collecting insects (24). Rapid adaptive responses, for instance a filtering out of 205 extreme phenotypes through increased mortality or reduced reproduction, might counteract 206 warming-induced losses of inter-individual synchrony. Such selection of the standing variation 207 can occur very rapidly, at least in herbaceous plants (25). 208 While our results show that climate warming causes a loss of phenological synchrony 209 among the individuals of a population, a study of leaf-out along elevational gradients in four 210 European tree species, between 1960-2016, revealed that leaf-out times at higher and lower 211 elevations are today compressed into a shorter time window compared to 58 years ago (26). 212 These findings do not contradict those of the present study because populations growing at high elevations were able to advance their phenology more than those at lower elevations for which 214 chilling and/or day-length requirements are no longer fulfilled (Fig. S8). As a result, the leaf-out 215 times of high-and low-elevation populations are converging (26). At the same time, however, 216 differences in day-length sensitivity (as well as chilling and temperature sensitivity) among the 217 individuals at any one elevation under climate warming are resulting in diverging flowering and 218 leaf-out times. 219 The overall prediction from the present findings is that human-caused climate warming is 220 leading to plant phenologies that are more heterogeneous within populations and more uniform 221 among populations (over altitude or latitude). The rapid loss of reproductive and vegetative 222 synchrony in European plants also predicts changes in their gene flow and trophic interactions, 223 although community-wide consequences are presently unknown. 224 225

Conclusion 226
The synchrony of developmental stages among organisms is a critical aspect of ecosystem 227 functioning. Here, based on massive ground observations and climate-manipulation experiments, 228 we show that global warming is altering within-population synchrony of leaf-out and flowering 229 dates in temperate plants, with warmer temperatures reducing inter-individual synchrony by up to 230 55%. Experiments suggest that individual genetic differences in the sensitivity to day-length 231 and/or winter chilling underlie the loss of synchrony, and future climate warming is expected to 232 further strengthen this trend. These results predict consequences for gene flow and trophic 233 interactions, but also emphasize the importance of adaptation when forecasting future plant 234 growth and productivity. 235 236

Materials and Methods 238 239
Analysis of leaf-out and flowering synchrony (LOS and FLS)  Chemische Industrie). For the two conifers Larix decidua and Picea abies leaf-out was defined as 247 the date when the first needles started to separate ("mouse-ear stage"; BBCH 10). Flowering was 248 defined as the date of beginning of flowering (BBCH 60). We removed (i) individuals, for which 249 the standard deviation of phenological observations across years was higher than 25 and (ii) leaf-250 out and flowering dates that deviated from an individual's median more than 3 times the median 251 absolute deviation (moderately conservative threshold) (26). 252

Analysis.
To test for an effect of spring temperature on inter-individual leaf-out synchrony (LOS) 253 and flowering synchrony (FLS), we divided the study area into pixels of one degree resolution 254 (~110 x 85 km), an area that can reasonably be considered as reflecting populations, at least for 255 wind-pollinated woody species (see discussion on herbs in the main text). To allow for within-256 pixel comparisons of LOS and FLS between years, data from the same individuals had to be used 257 each year. To achieve this, we kept only pixels for which there were at least three individuals 258 with data for the same 15 years. For each pixel, we deleted all (i) individuals growing at altitudes 259 that deviated by >200 m from the average altitude of all individuals within the pixel, and (ii) 260 years that had less than 90% plant-coverage, i.e., data from at least 90% of the individuals within the pixel had to be available for the respective year, otherwise the year was excluded from the 262 analysis. This data cleaning left us with a total of 12,536 individuals, 317,672 phenological 263 observations (individuals x year), and a median time-series length of 25 years (minimally 15 264 years, maximally 48 years). The number of individuals within pixels (per species and 265 phenological stage) ranged between 3 and 53 (median = 12). See Figs. S1b and S2b for 266 information on the number of pixels used per species. 267 For each year and species, LOS and FLS within pixels were then calculated as the requirements. Individual forcing requirements until leaf-out were calculated as the sum of degree-273 days (DD) from 1 January until leaf-out or flowering using 5°C as base temperature (e.g., ref.

276
where DD sum is the accumulated degree days until leaf unfolding, t LO is the day of leaf unfolding, 277 T t is the mean daily temperature on day t, and t 0 is the start date for forcing accumulation, which 278 was fixed at 1 January. For each year and species, LOS-DD and FLS-DD within pixels were then 279 calculated as the standard deviation of forcing requirements until leaf-out or flowering dates. 280 The daily mean air temperature at each site was derived from a gridded climatic data set 281 of daily mean temperature at 0.5º spatial resolution (approximately 50 km, ERA-WATCH) (28). 282 For each year, preseason temperature within pixels was defined as the average temperature 283 during the 60 days prior to the average leaf unfolding or flowering date within the respective 284 pixel, which is the period for which the correlation coefficient between phenological event and 285 temperature is highest (29). 286 To test if shortened photoperiods and/or reduced winter chilling explain the decrease in 287 phenological synchrony under warmer preseasons, for each year, pixel, and species, we 288 calculated the average chilling hours until leaf-out or flowering and the average photoperiod (PP) 289 at the date when the average forcing requirements until leaf-out or flowering were fulfilled. 290 Chilling hours were calculated on basis of 6-hourly temperature data (CRU-NCEP, spatial 291 resolution of 0.5°; https://crudata.uea.ac.uk/cru/data/ncep/), as the sum of hours from 1 292 November until leaf-out/flowering with an average temperature between 0°C and 5°C (e.g., ref 293 where Ch sum is the sum of chilling hours until leaf unfolding, t LO is the day of leaf unfolding, T t is 296 the hourly mean temperature on hour t, and t 0 is the start date for chilling accumulation, which 297 was fixed at 1 November in the year before leaf unfolding. 298 PP was calculated as a function of latitude and DOY (30)

Statistical analyses. 308
Within each pixel we applied linear models to test for an effect of preseason temperature, 309 photoperiod, and winter chilling on phenological synchrony (LOS, LOS-DD, FLS and FLS-DD). 310 We then determined the frequency distributions for the correlation coefficients between 311 phenological synchrony and preseason temperature across all species and sites. For each species, 312 we applied t-tests to detect whether the average of all correlation coefficients obtained for each 313 pixel differs from zero. To model changes in the distribution of within-population leaf-out and 314 flowering dates (means and standard deviations) in response to temperature, we applied mixed-315 effects models using average leaf-out / flowering dates or LOS / FLS as response variables, 316 preseason temperature as explanatory variable, and site as a random effect to control for the use 317 of different sites in the model. 318 To test for the relative effects of preseason temperature on (i) inter-individual variation 319 in leaf-out/flowering date (LOS / FLS) and (ii) inter-individual variation in forcing requirements 320 until leaf-out/flowering (LOS-DD / FLS-DD) we applied hierarchical Bayesian models. To test 321 for the effects of winter chilling and day-length on phenological synchrony, we applied 322 hierarchical Bayesian models including both winter chilling until leaf-out and day length at the 323 date when the average forcing requirements until leaf-out or flowering were fulfilled as predictor 324 variables. The use of a Bayesian framework allowed us to fit slope parameters across traits 325 simultaneously without concerns of multiple testing or P-value correction. All models included 326 random effects for (i) species (to address within-species rather than between species phenological 327 synchrony) and (ii) pixels (to address within-population rather than between-population 328 phenological synchrony). To allow for direct effect size comparisons, all continuous variables were standardized by subtracting their mean and dividing by 2 SD before analysis (31). The 330 resulting posterior distributions are a direct statement of the probability of our hypothesized 331 relationships. Effective posterior means ± 95% confidence intervals are shown in Fig. 2. 332 To parameterize our models, we used the JAGS implementation (32)  Materials Table 1 for a list of species). An individual was scored as having leafed out when at 344 least three branches had unfolded leaves pushed out all the way to the petiole (37). To test 345 whether the trends observed in the PEP analysis are consistent with our common garden data, the 346 same parameters (LOS, LOS-DD, preseason temperature, winter chilling, and photoperiod) were 347 calculated as described above (Analysis of leaf-out and flowering synchrony (LOS and FLS) 348 using the PEP database). We then applied hierarchical Bayesian models including species 349 random effects (see paragraph above) to test for the effects of preseason temperature, winter 350 chilling, and day-length on LOS and LOS-DD. 351

Twig cutting experiments and phenological scoring 354
To study the extent of intraspecific variation in leaf-out strategy (within-species variation in 355 photoperiod, chilling, and forcing requirements) and its implications under climate warming, we 356 Temperatures inthe climate chambers were held at 12°C during the night and 20°C during 367 theday, with an average daily temperature of 16°C to simulate forcing temperatures. Illuminance 368 in the chambers was about 8 klux (∼100 µmol s −1 m −2 ). Relative air humidity was held between 369 40%and 60%. To account for within-individual variation, we used 10 replicate twigs per 370 individual treatment and monitored bud development every second day. For each individual and 371 treatment, we then calculated the mean leaf-out date out of the first eight twigs that leafed out. A 372 twig was scored as having leafed out when three buds had unfolded leaves pushed out all the way 373 to the petiole (37). Forcing requirements until leaf-out were calculated as the sum of degree-days 374 [outside of and in climate chambers] from 10 December (1 st collection date) until leaf-out using 5°C as base temperature (e.g., ref. 27). Chilling hours were calculated as the sum of hours from 1 376 November until leaf-out with an average temperature between 0°C and 5°C. 12°C during night to 20°C during day, with an average daily temperature of 16°C. Day length in 384 the chambers was set to 8h, 12h, or 16h. 385 Individual photoperiod sensitivity was defined as the slope of the function between day-386 length treatment and accumulated degree days (>5°C) until leaf-out (twigs were collected on 21 387 March; see Fig. 3b). The steeper the slope, the stronger the effect of photoperiod on the amount 388 of warming required for leaf-out. A flat slope indicates that photoperiod has no effect on the 389 timing of leaf-out. 390 Individual chilling sensitivity was defined as the slope of the function between chilling 391 treatment (collection date) and accumulated degree days (>5°C) until leaf-out when twigs were 392 kept under constant 16-h day length (see Fig. 3c). The steeper the slope, the stronger the effect of 393 chilling on the amount of warming required for leaf-out. 394 Individual forcing requirement was defined as the accumulated degree days (>5°C) until 395 leaf-out under long chilling (21 March collection) and constant 16-h day length (see Fig. 3a). 396 Under such conditions, chilling requirements and photoperiod requirements should be largely 397 met, and thus the remaining variation in leaf-out dates should be largely attributable to 398 differences in forcing (warming) requirements.
In winter 2015/2016, twigs from the same 11 individuals were harvested every two weeks (from 401 10 December until 21 March) and kept under the same temperature conditions applied in 402 experiment 1 (12°C during night to 20°C during day), with natural day length. This allowed us to 403 test if those individuals with no/little photoperiod sensitivity would advance their leaf-out more 404 under short winter conditions than photoperiod-sensitive individuals, and to determine the 405 relative effect of individual variation in photoperiod requirements, chilling requirements and 406 forcing requirements on leaf-out variation under different winter/spring conditions (Fig. 4). 407 Within-species leaf-out synchrony (LOS) was calculated as the standard deviation of individual 408 leaf-out dates. To analyze which leaf-out cues (photoperiod, chilling, and forcing requirements) 409 best explain leaf-out variation among individuals, we applied a multivariate linear model, 410 including individual forcing, photoperiod, and chilling requirements (as inferred from experiment 411 1) as explanatory variables. To express the total variation in leaf-out dates that can be attributed 412 to each trait, we used ANOVA sums of squares (see Fig. 4d). 413 To infer which percentage of the variation in leaf-out dates is due to treatment effects, 414 between-individual variation, or within-individual variation, we calculated variance components 415 by applying a random-effects-only model including treatments and individuals as random effects 416 (individuals nested within treatments). Results show that of the total leaf-out variation among 417 twigs, 52% can be explained by between-individual variation, 33% by treatments, and only 15% 418  sylvatica and (f)   chilling on inter-individual leaf-out synchrony. Hierarchical Bayesian linear models were applied 614 using information on 13 (upper), 9 (middle), and 59 species (lower panels). To account for 615 within-species rather than among-species synchrony, all models include species random effects. 616 The models using the PEP data (upper and middle panels) additionally include site random 617 effects (1° pixels) to address within-population phenological synchrony. All variables were 618 standardized to allow for direct effect size comparisons.

Leaf-out Flowering
Figure S7 | In species in which preseason temperature has little effect on the mean flowering date, preseason temperature also has little effect on FLS. a, Positive correlation between species' mean temperature sensitivity of flowering date (days advance in flowering per one degree warming) and the mean temperature sensitivity of FLS (increase in the standard deviation of inter-individual flowering times per one degree warming). b, No correlation between species' mean temperature sensitivity of leaf-out date (days advance in leaf-out per one degree warming) and the mean temperature sensitivity of LOS (increase in the standard deviation of inter-individual leaf-out times per one degree warming). The effects of preseason temperature on mean flowering date, mean leaf-out date, FLS, and LOS were inferred from mixed-effects models including site as a random effect.  Figure S8 | Schematic representation of within-and among-population phenological synchrony in response to climate warming. As demonstrated in this study, inter-individual synchrony within a population will decrease under warmer preseason temperatures because individuals differ in their sensitivity to temperature. Within-population variation under ambient or warmed preseason temperatures is illustrated by the solid blue and red arrows, respectively. By contrast, phenological synchrony among populations is expected to increase, given that populations in warm regions (Population 3) will advance their phenology less than populations in cold regions (Population 1). This is illustrated by the dashed blue and red arrows, showing that the difference in the average phenological date between Population 1 and 3 is smaller under warmer preseasons (red dashed arrow) than under ambient preseason temperatures (blue dashed arrow).  Figure S9 | Percent variation in leaf-out dates attributable to treatment effects and between-and within-individual variation within treatments. Data from experiment 2 (see Methods). Variance components were inferred from random-effects-only models, including leafout date as the dependent variable and treatment and individuals as nested random effects.