Latitudinal clines in the phenology of floral display associated with adaptive evolution during a biological invasion

26 Premise: The timing and pattern of flowering -- flowering phenology -- strongly influences 27 reproductive success in plants. Days to first flower is easy to quantify and widely used to 28 characterize phenology, but reproductive fitness depends on the full schedule of flower 29 production over time. 30

We examined floral display evolution associated with rapid adaptive evolution and 31 range expansion among thirteen populations of Lythrum salicaria populations, sampled along a 32 10-degree latitudinal gradient in eastern North America. Growing these collections in a common 33 garden field experiment at a mid-latitude site, we used Principal Coordinates Analysis (PCoA) 34 and quantitative metrics analogous to central moments of probability distributions (i.e., mean, 35 variance, skew, and kurtosis) to quantify flowering schedules.. 36 Key Results: Consistent with adaptation to shorter growing seasons at higher latitudes, we found 37 that populations from higher latitudes had earlier start and mean flowering day for both 38 population average (i.e., among individuals) and population aggregate (i.e., pooled flowers). 39 Clines in population average and aggregate flowering skew increased with latitude while kurtosis 40 decreased, consistent with a bet-hedging strategy arising in biotic environments with more 41 herbivores and greater competition for pollinators. Although clines were generally similar at the 42 average vs aggregate level, we found that individual flowering schedules from more southern 43 populations tended to flower later than their population aggregate. 44

Heritable genetic changes in flowering schedule are consistent with biotic and 45 abiotic responses to natural selection, potentially contributing to rapid evolution and range 46 expansion of this invasive species. 47 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 3

48 Phenology affects how organisms interact with their environment, generating important 49 consequences for their survival and reproduction. As targets of natural selection and 50 determinants of population vital rates, the timing of growth, reproduction, and other phenological 51 traits are crucial to understanding eco-evolutionary dynamics. For example, in sexual outcrossing 52 plants, flowering phenology governs interaction with mutualists (Fenster et al., 2004), 53 antagonists (Strauss & Whittall, 2006) and abiotic factors (Jin et al., 2019). Consequently, plants 54 experience spatial and temporal variation in natural selection on flowering phenology (Kudo, 55 2007), sometimes resulting in the evolution of local adaptation (Anderson et al., 2012; Fitchett et 56 al., 2015; D. W. Inouye, 2022; Panchen, 2022; Panchen & Gorelick, 2017; Prather et al., 2023). 57 As such, flowering phenology is a useful composite trait for investigating evolutionary responses 58 to past and future environmental challenges. 59 Days to first flower (i.e., flowering time) is a simple yet informative metric with which to 60 study the ecology and evolution of phenology in plant populations. For example, latitudinal and 61 altitudinal clines in days to first flower are common to a variety of plant taxa (Alexander et al., 62 2009; Colautti & Lau, 2015; Halbritter et al., 2018), and reciprocal transplant experiments have 63 confirmed that these clines are an important component of local adaptation (Ågren & Schemske, 64 2012; Anderson & Gezon, 2015; Colautti & Barrett, 2013; Ensing & Eckert, 2019; Griffith & 65 Watson, 2005). Flowering time clines are consistent with local adaptation models characterized 66 by stabilizing selection within populations but a shift in optimal flowering time along 67 environmental gradients (de Villemereuil et al., 2020; Gauzere et al., 2020; Kirkpatrick & 68 Barton, 1997). In contrast to the prediction of stabilizing selection, empirical measurements of 69 selection suggest that flowering time is often under directional selection within populations 70 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 4 (Munguía-Rosas et al., 2011). Moreover, long-term phenology records show a variety of 71 responses to climate change -- from earlier to later to no change in flowering time, depending on 72 the species (Rafferty et al., 2020; Wolkovich et al., 2012). Thus, there appears to be 73 disagreement between theoretical models and empirical observations of geographic clines and 74 local adaptation on the one hand and observed measurements of natural selection and phenotypic 75 changes in natural populations on the other (Austen et al., 2017). Antagonistic selection on traits 76 correlated with the onset of flowering offers one potential resolution to this apparent 77 contradiction. 78 Day of first flowering may trade off with other aspects of flowering phenology such as 79 the shape of the flowering schedule. These other aspects of flowering phenology may not be well 80 represented by flowering time (Newstrom et al. 1994; CaraDonna et al., 2014; B. D. Inouye et 81 al., 2019) but have important effects on reproductive fitness (Fox, 2003). The shape of the 82 flowering schedule can be quantified by metrics analogous to central moments of probability 83 distributions, which include the mean, variance, skew, and kurtosis (Box 1). Not to be confused 84 with the central moments of a statistical population or sample, the analogous characteristics of a 85 flowering schedule quantify temporal variation in relative reproductive investment through time. 86 These characteristics use the same mathematical equations (Table 1) but are not shaped by 87 probabilistic process and therefore may not correlate with onset or duration of flowering (Fig. 88 1A). These metrics have rarely been used in studies of flowering phenology, yet they have the 89 potential to offer a novel perspective on selection and trade-offs affecting the evolution of 90 phenology. 91 Central moments of individual flowering schedules describe the timing and intensity of 92 reproductive investment, and thus are likely to experience strong selection from abiotic or biotic 93 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 5 features of the environment. We are not aware of studies quantifying evolutionary divergence of 94 entire flowering schedules among naturally occurring populations, but there are reasons to expect 95 populations to evolve differences. For example, individuals that experience competition for 96 pollinators later in the growing season may benefit from showy inflorescences that enhance 97 pollinator attraction (O’Neil, 1997; Thomson, 1980) and be characterized by a flowering 98 schedule with a more positive skew and longer duration. Herbivores and seed predators can be 99 unpredictable, which may select for a bet-hedging strategy and possibly delay in the onset of 100 flowering (Augspurger, 1981; Elzinga et al., 2007; Wright & Meagher, 2003). These conditions 101 might further result in a longer flowering duration with a higher variance and negative skew in 102 the flowering schedule if late-season ramp-up in flowering is limited by harsh growing 103 conditions at the end of a season. 104 In addition to biotic interactions, flowering schedules are likely to be shaped by abiotic 105 factors that determine the timing and duration of environmental conditions favorable to the 106 development of reproductive organs. Shorter growing seasons typically favour an earlier onset of 107 flowering (Austen & Weis, 2015; Colautti & Barrett, 2013; Griffith & Watson, 2005), and may 108 favour reproduction only within a shorter growing season; this would result in flowering 109 schedules with lower duration and variance. In marginal environments, plants may combine 110 strategies, concentrating flower production into a narrow window but also producing a few 111 flowers over a long timeframe as a bet-hedging strategy (Forrest & Thomson, 2010; Simons, 112 2011; Tufto, 2015), which would produce a long-tailed leptokurtic or skewed flowering 113 schedule. The above examples illustrate the complex ways that biotic interactions and the abiotic 114 environment could shape selection on individual flowering schedules, and how these differences 115 can be quantified using equations developed for probability distributions. 116 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 6 Natural selection acts on phenological schedules of individuals, but variation among 117 individuals within a population affect eco-evolutionary interactions with pollinators, herbivores, 118 and potential mates. Just as individuals have flowering schedules, the aggregated flower 119 production of all individuals within a population can be characterized (Box 1), which could yield 120 emergent properties that are relevant to understanding eco-evolutionary feedbacks. We define 121 emergent properties here as deviations in the aggregate population from the average of individual 122 schedules comprising it (Table 1). For example, a population can have a high aggregate variance 123 but low average individual variance in the flowering schedule if individuals vary in onset and 124 duration of flowering. In cases like this, the central moments of the aggregate flowering schedule 125 show characteristics that are emergent and distinct from the typical individual comprising it (Fig. 126 1). Although emergent properties can arise from adaptive processes acting on individuals, they 127 can also have important ecological consequences for species communities. For example, 128 population aggregate flowering schedules are sometimes positively skewed with an abrupt start 129 and a long tail of flowering (Blionis et al., 2001; Rabinowitz et al., 1981), which may be an 130 adaptation to enhance pollinator visitation to the population (Stemkovski et al., 2023; Thomson, 131 1980). Additionally, strong frequency-dependent selection may concentrate flowering within a 132 narrow but synchronous timeframe to attract more pollinators (Rathcke & Lacey, 1985), 133 increase mating opportunities (Ison & Weis, 2017), or prevent interspecific cross-pollination 134 (Fantinato et al., 2018). Therefore, in pollinator or mate-limited environments, population 135 aggregate flowering schedules may be very similar to the average flowering schedule. Thus, 136 emergent properties (or a lack thereof) have the potential to provide insight into eco-evolutionary 137 processes and ecological consequences of aggregate flowering schedules. 138 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 7 Environmental gradients associated with latitude can be useful for studying phenotypic 139 evolution because selection often varies predictably due to shorter growing seasons and weaker 140 biotic interactions towards polar latitudes (reviewed in Zvereva & Kozlov, 2021). Adaptation 141 along gradual environments is commonly manifested as clines in phenotypic traits that are 142 observable when genetically related individuals are grown in a common garden environment 143 (Olsson & Ågren, 2002; Panchen, 2022; Sobral et al., 2013). For example, there is good evidence 144 for the rapid evolution of adaptive flowering time clines in the perennial wetland herb Lythrum 145 salicaria (purple loosestrife) during its invasion of North America. Selection splines (Colautti & 146 Barrett, 2010), trait correlations (Colautti & Barrett, 2011), reciprocal transplants (Colautti & 147 Barrett, 2013), and herbarium records (Wu & Colautti, 2022) collectively support a model of 148 constrained selection in which shorter growing seasons toward the northern range limit favour 149 earlier flowering. In contrast, longer growing seasons in the southern part of the introduced range 150 relax selection on flowering time and favors plants that flower later but grow larger in stature. 151 However, by focusing on the first day of flowering, previous studies have ignored other aspects 152 of the flowering schedule that may be important for understanding adaptation and constraint 153 during evolution and invasion of L. salicaria. 154 Here, we investigate clines in flowering schedules among thirteen populations sampled 155 across ten degrees of latitude in eastern North America and grown in a common garden field 156 study at the Koffler Scientific Reserve north of Toronto, Ontario, Canada. We quantify both 157 individual and population aggregate flowering schedules to address the following questions: (1) 158 How well do common phenological metrics, like the onset and duration of flowering, predict 159 often overlooked features of the flowering schedule such as the variance, skew and kurtosis? (2) 160 Are there latitudinal clines in flowering schedules that are consistent with the adaptive 161 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 8 hypotheses outlined above? Specifically, we predicted that longer growing season lengths, 162 stronger competition for pollinators, and higher rates of herbivory in southern populations will 163 have favoured genotypes with: (i) a longer flowering duration and (ii) a more negative kurtosis 164 with a concentrated peak and rare flowering events far from the average flowering day, but also 165 (iii) an increase in the flowering schedule variance and (iv) a more negative skew with flowering 166 focused later in the schedule. (3) Are there emergent properties of population aggregate 167 flowering schedules that differ from the population average of individual flowering schedules? 168 Because natural selection acts on individual flowering schedules, emergent properties represent 169 (i) potential biases in evolutionary inference when aggregate characteristics are used in place of 170 individual or population average characteristics of flowering schedules, and (ii) ecological 171 phenomena that are shaped by selection and potentially affect herbivore and pollinator 172 communities. 173 174

175 Experimental design -- 176 Details of the field experiment are provided in Colautti & Barrett (2010) and the seeds used were 177 collected along a 10-degree latitudinal gradient in eastern North America, as described in 178 Montague et al. (2008). Briefly, 13 populations were chosen to represent a latitudinal gradient 179 from Timmins, Ontario (48.48 oN) to Easton, Maryland (38.75 oN) and were grown at the 180 Koffler Scientific Reserve at Jokers Hill, located approximately at the geographic center of the 181 latitudinal gradient (44.03 oN). Prior to inflorescence development, plants were sprayed with 182 insecticide to reduce incidence of herbivorous insects. The primary purpose of this experiment 183 was to characterize flowering time clines and genetic trade-offs in a multi-year (2006-2008) 184 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 9 common garden study, as reported in Colautti & Barrett (2010, 2011). In 2008 we additionally 185 recorded all open flowers on every stem of a subset of 369 individuals of 216 unique seed 186 families representing the 13 populations (Appendix S1). Our daily census records of this subset 187 of plants occurred over 125 days on a rotating schedule such that each individual was surveyed 188 every five days from the onset of flowering in July until September. A 5-day cycle was chosen 189 because stigma receptivity is nearly zero after 24 h of anthesis (Waites & Ågren, 2006), after 190 which unpollinated flowers begin to senesce and abscise while pollinated flowers begin to close 191 and fuse at the calyx, limiting the risk of counting the same flowers on multiple census dates. All 192 calculations and statistical analyses were conducted in R version 4.2.2 (R Core Team, 2022) in R 193 Studio v4.2.2 (RStudio Team, 2022), and are available with the original data from the Dryad 194 database (Dryad DOI: TBD(ROB)). 195 196 Characterizing flowering schedules -- 197 We characterized the flowering schedules of each individual plant using the equations shown in 198 Table 1, standardized to the number of calendar days after the first flower was observed in the 199 experiment. The number of fruits produced over a growing season is variable in this experiment 200 because individuals that are locally adapted to mid-latitude balance flowering time with growth 201 to produce more flowers compared to small, early flowering phenotypes more characteristic of 202 northern populations or large, late-flowering phenotypes characteristic of more southern 203 populations. To account for variation in total flower production, we standardized the number of 204 open flowers every five days for each individual as a proportion of the total number of open 205 flowers censused over 125 days of observation. We then calculated a weighted census day by 206 multiplying each day by the proportion of flowers open on that day. The proportion of flowers 207 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 10 and the weighted census days were summed to calculate the weighted mean and other 208 characteristics of individual flowering schedules analogous to central moments (variance, skew 209 and kurtosis) (Box 1; Table 1). In addition to calculating the characteristics of every individual’s 210 flowering schedule, we calculated the average of each characteristic for each population. 211 In addition, the start of flowering, duration, and central moments, we used a Principle 212 Coordinates Analysis to characterize flowering schedules of individuals as described by Austen 213 et al., (2014). Briefly, we computed a dissimilarity matrix from raw flower counts using 214 Kolmogorov-Smirnov (KS) distances (Sokal & Rohlf, 1995) where each matrix cell was a KS 215 statistic measuring the maximum distance between two flowering schedules. We chose KS 216 distance because it is robust to differences in flower number and census schedules. However, the 217 KS statistic is sensitive to differences in scale, so flowering durations were standardized to run 218 from 0 to 1 by subtracting the shortest duration time from each individual and dividing by the 219 difference of the longest and shortest durations (Austen et al., 2014). We used the KS matrix as 220 the input for the cmdscale function to generate the PCoA axes. Since KS distances do not meet 221 the triangle inequality condition of PCoA, we consider only eigenvectors with positive 222 eigenvalues (Austen et al., 2014; Legendre & Legendre, 1983). It is important to note that the 223 PCoA axes characterize variation in the shape of flowering schedules among individuals but not 224 variation in onset and duration of flowering. 225 In addition to individual flowering schedules, and the average of all individuals within 226 each population, we calculated characteristics of the aggregate population flowering schedules. 227 These calculations use the same equations after first pooling the flower counts for each day 228 across all individuals within a population (Table 1). For example, the aggregate start of flowering 229 of a population would be the first day that any individual in the population flowered, and 230 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 11 aggregate duration of flowering would be the number of days until the last open flower was 231 observed in the population. 232 We used bootstrap models with 999 iterations to compute a bootstrapped mean and 95% 233 confidence intervals for both population average and population aggregate flowering schedules. 234 In each iteration, individuals were hierarchically sampled with replacement from within each 235 population, and then each metric in Table 1 was calculated either from the pooled set of 236 individuals in each population to calculate the population aggregate characteristics, or for each 237 individual within the population to calculate the population average characteristics. Original data 238 with reproducible R code for the bootstrap models are available from the Dryad database (Dryad 239 DOI: TBD(ROB)) 240 241 C haracteristic correlations, clines and emergent properties -- 242 Our first research question addresses whether conventional metrics (i.e., onset and duration of 243 flowering) correspond to more detailed characteristics of the flowering schedule (i.e., mean, 244 variance, skew, kurtosis, and PCoA axes). To investigate these relationships, we calculated 245 pairwise Pearson correlations to compare three scales of analysis, namely (A) average schedules 246 of each population (n = 13), (B) the deviation of each individual from its population mean (n = 247 369), and (C) aggregate schedules calculated on pools of individuals for each population (n = 248 13). In addition to comparing conventional metrics with other flowering schedule characteristics, 249 the correlational structure of these matrices can also be instructive for testing evolutionary 250 hypotheses because evolutionary divergence of flowering schedules would produce correlations 251 among average (A) or aggregate (C) populations that are not evident among individuals within 252 populations (B). 253 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 12 We used simple linear regression models of bootstrap means to address our second 254 question concerning whether there were latitudinal clines characteristic of both average and 255 aggregate flowering schedules. Each univariate linear regression included the bootstrap mean of 256 a flowering schedule characteristic as the dependent variable, with latitude as a linear and 257 quadratic predictor, and significance was tested using the anova function in R to find the 258 minimal adequate model. 259 To quantify emergent properties of population aggregate flowering schedules that differ 260 from the population average, we subtracted the average value of each population from the 261 aggregate value for the same population. Therefore, a positive value for an emergent property 262 means that a population had an aggregate value that was higher than the average individual 263 comprising it. For example, all populations will have a negative emergent property for onset of 264 flowering because the average onset of flowering will never be earlier than the day of first flower 265 in the aggregate population. As with average and aggregate characteristics, we tested for clines in 266 emergent properties using linear models as described in the previous paragraph. 267 Finally, we also considered the false discovery rate between 6 and 24 separate statistical 268 models. We consider a range of degrees of freedom because there are 24 separate statistical 269 models to test for clines (i.e., six characteristics across three types of comparisons) but the 270 aggregate, average, and emergent characteristics are each calculated from the same population 271 data, resulting in no more than 6 independent characteristics. 272 273 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 13

274 Flowering schedule characteristics -- 275 Pearson product-moment correlations revealed significant relationships among commonly 276 measured traits of the flowering schedule (i.e., start and duration) and more detailed 277 characteristics (i.e., mean, variance, skew, and kurtosis), across three types of comparison (Fig. 278 2). At the aggregate level, populations that started flowering later also had a later mean flowering 279 date. Additionally, populations with longer flowering duration had a higher coefficient of 280 kurtosis indicating more flowering far from the mean flowering date (Fig. 2C). In contrast, there 281 was no significant correlation between start date and duration among population average 282 flowering schedules, but there was a significant negative correlation with skew indicating that 283 populations with a later start date, on average, tend to have individuals with imbalanced 284 flowering clustered more at later dates (Fig. 2A). 285 Flowering schedule variance was significantly positively correlated with flowering 286 duration on average but not among aggregate flowering schedules. Residual, within-population 287 correlations were generally weaker when compared to correlations among population average 288 and aggregate flowering schedules (Fig. 2B). In addition, the correlation between onset and 289 duration of flowing among individuals within populations was in the opposite direction to 290 population average and aggregate. In other words, plants that started flowering earlier relative to 291 other individuals within the same population also flowered for a shorter duration (Fig. 2B), 292 whereas populations with earlier flowering individuals, on average, also tended to have 293 individuals with longer flowering durations (Fig. 2A). 294 295 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 14 Principal Coordinates Analysis -- 296 In addition to the above characteristics, we use Principal Coordinates Analysis (PCoA) to fully 297 characterize similarities and differences among flowering schedules after controlling for 298 variation in start date and duration. The PCoA of all 369 individuals resulted in 369 independent 299 eigenvectors, of which 192 had positive eigenvalues. However, of these 192, only the first two 300 PCoA axes had eigenvectors that explained more than 10% of the variation in flowing schedule 301 shape, and together they explained 35% of the variation (Table 2). By comparison, the next ten 302 eigenvectors combined explained another 35% of variation, with the remaining 180 accounting 303 for 30% of the variation. Thus, we next focus on the first two eigenvectors. 304 For both the population aggregate and average, the first axis (PCoA1) was positively 305 correlated with skew and negatively correlated with the mean day of flowering (Fig. 2). PCoA1 306 was also negatively correlated with duration of flowering when comparing population average 307 but not population aggregate flowering schedules. In contrast, within-population residual 308 correlations involving PCoA1 were generally weaker than among population average or 309 aggregate, with significant correlations for mean, skew and kurtosis (Fig. 2B). Correlations 310 involving PCoA2 were quite different than PCoA1 within and among all three comparisons. For 311 example, PCoA2 was not significantly correlated with any other metrics in aggregate 312 distributions (Fig. 2C) but was negatively correlated with start and kurtosis and positively 313 correlated with skew and duration among individuals within the same populations (Fig. 2B) and 314 was positively correlated with variance and duration of average distributions (Fig. 2A). By 315 mathematical definition of the PCoA, there was no correlation between PCoA1 and PCoA2 316 among all individuals, nor was there a significant correlation among residuals of individuals 317 within populations (Fig. 2B). In other words, PCoA1 and PCoA2 change in a correlated way. 318 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 15 This can be seen clearly in Fig. 3, where northern populations cluster at intermediate values of 319 PCoA1 and PCoA2 compared to mid-latitude (higher PCoA1 but lower PCoA2) and southern 320 populations (lower PCoA1 and lower PCoA2). 321 322 Latitudinal clines and emergent properties -- 323 The correlated changes in population aggregate and average flowering schedules are reflected in 324 significant changes with latitude. For example, the start, mean, and kurtosis of flowering 325 schedules were each significantly negatively related to latitude for both aggregate and average 326 schedules (Table 3, Fig. 4). This means that southern populations start flowering later, with a 327 later average flowering date, and more ‘outlier’ flowers opening far from the average flowering 328 date, both in aggregate and on average, when compared to northern populations. However, 329 aggregate and average schedules vary in two key characteristics – aggregate but not average 330 schedule duration declined with latitude, whereas average but not aggregate skew was negative 331 at lower latitudes (Fig. 4). In other words, individuals have similar flowering durations on 332 average but are more spread out in southern populations, resulting in a longer aggregate duration 333 towards the south. Similarly, individual schedules are more heavily weighted toward later 334 flowering dates towards the south, but variation in start dates average out such that aggregate 335 schedules are more evenly distributed and not significantly related to latitude. 336 To quantify emergent properties, we compared aggregate to average schedules and this 337 procedure revealed considerable variation among populations and two latitudinal clines (Table 4; 338 Fig. 4). Although the aggregate start of flowering must always be earlier than the average, 339 resulting in a negative emergent property score, we found a significant positive correlation with 340 latitude. Likewise, the duration of flowering must always be longer in aggregate than on average, 341 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 16 but we found a marginal negative relationship with latitude (Fig. 4, right column). Together, 342 these clines demonstrate more synchronous flowering schedules in northern populations 343 compared to those occurring in the south. We also found a nonlinear cline in variance 344 representing a lower aggregated flowering schedule variance than the average individual 345 variance at intermediate latitudes. 346 Overall, we expect about one false positive cline at α < 0.05 with 24 independent 347 statistical tests. But of the 10 significant clines, only five remained significant after a Bonferroni 348 adjusted α = 0.0083 for six independent tests or α = 0.0021 for 24 tests. Despite uncertainty in 349 the detection of clines, bootstrap estimates show a relatively high level of precision for most 350 population means, as shown by the small error bars in Fig. 4. 351 352

353 There is growing evidence for rapid evolution in invasive species (van Kleunen et al., 2018), but 354 it remains uncertain how often adaptive evolution contributes to their spread (Colautti & Lau, 355 2015). Rapid evolution of a latitudinal cline in flowering time in response to local growing 356 conditions likely increased the rate and extent of invasion by Lythrum salicaria across North 357 America (Colautti & Barrett, 2013; Wu & Colautti, 2022). However, previous studies focusing 358 on the onset and duration of flowering provided only a partial view of how plants allocate 359 flowering resources (Box 1), raising uncertainty about the extent to which phenological 360 evolution has been adaptive. We used a common garden study with detailed metrics to 361 characterize entire flowering schedules and tested hypotheses about adaptative latitudinal clines. 362 We found that the onset and duration of flowering were not good predictors of other 363 characteristics of individual flowering schedules (Fig. 2B). In contrast, the average flowering 364 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 17 schedule characteristics were highly correlated among population averages (Fig. 2A) and also 365 among population aggregate flowering schedules (Fig. 2C). Furthermore, that these among-366 population correlations were largely due to correlated evolution along latitudinal clines (Fig. 4). 367 Emergent properties were highly variable among populations and formed clines indicative of 368 more synchronized flowering schedules in northern populations. Below we discuss the 369 implications of these results for understanding the evolutionary causes and ecological 370 consequences of flowering schedule diversity. 371 372 Evidence for correlated evolution -- 373 Date of first flower and duration of flowering are relatively easy to measure and are commonly 374 reported as focal traits in phenological studies. However, reproductive fitness depends on the full 375 schedule of flower production, and it has not been clear how well the onset or duration of 376 flowering correlates with other characteristics of the flowering schedule. In L. salicaria, plants 377 that flowered later also flowered for a longer duration with a higher flowering variance and 378 slightly lower skew, relative to other individuals sampled from the same population (Fig. 2B). 379 This result is consistent with a life history trade-off between flowering earlier at a small size or 380 delaying reproduction to accumulate more resources for later floral investment – a mechanism 381 that was proposed previously for this species (Colautti & Barrett, 2013; Wu & Colautti, 2022) 382 but was not validated with floral display data until now. 383 In contrast to individual flowering schedules, measurements of the flowering schedules of 384 whole populations provide more limited insights into the evolutionary mechanisms shaping 385 flowering schedules. The different correlational structure for average (Fig. 2B) versus aggregate 386 schedules (Fig. 2C) produced significant variation in emergent property characteristics among 387 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 18 populations (Fig. 4C), suggesting that aggregate characteristics are not generally indicative of 388 individual flowering schedules within a population. Studies that have not considered the full 389 schedule shape of individual plants may therefore neglect key aspects of phenology that are 390 crucial for investigating ecological and evolutionary hypotheses. For example, there is 391 considerable among-population and among-species variability in phenological shifts due to 392 climate change (Anderson et al., 2012; Fitchett et al., 2015; D. W. Inouye, 2022; Panchen, 2022; 393 Panchen & Gorelick, 2017; Prather et al., 2023). This variation may be explained by selection on 394 characteristics of flowering schedules that are not correlated with the onset of flowering. 395 There were differences in both the strength and direction of correlations among 396 individuals within L. salicaria populations when compared with population average and 397 aggregate characteristics (Fig. 2). Among-population correlations were caused by either 398 statistical artifacts or genetic constraints, then we would expect to see similar correlations among 399 individuals within populations. Instead, the differing correlation matrices are more consistent 400 with correlated evolution among populations. We cannot rule out the possibility that 401 environmental effects mask genetic correlations within populations; however, this would require 402 a much larger experiment with sufficient replication to estimate correlations among relatives 403 within populations (Price et al., 1988; Rausher, 1992). Nor can we rule out the potential role of 404 stochastic processes without replicated sampling from independent geographical gradients 405 (Colautti & Lau, 2015). Nevertheless, our finding that flowering schedules change continuously 406 with latitude lends additional evidence for the correlated evolution hypothesis. Moreover, the 407 direction of the observed clines are consistent with adaptive responses to biotic and abiotic 408 selection, as discussed below. 409 410 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 19 Latitudinal clines and emergent properties -- 411 Latitudinal clines in the onset of flowering have evolved independently in several species across 412 North America (McGoey et al., 2020; Montague et al., 2008; Samis et al., 2012; Wu & Colautti, 413 2022). In our study of L. salicaria, we additionally found latitudinal clines in average mean, 414 skew and kurtosis of flowering schedules (Fig. 4, left column). These types of clines indicate that 415 individuals from southern environments have floral displays that peak later in the season but also 416 a few open flowers far from the mean flowering date. These findings correspond to higher 417 intensity of biotic interactions in southern environments reported in a meta-analysis of terrestrial 418 ecosystems (Zvereva & Kozlov, 2021). 419 Lythrum salicaria is a self-incompatible outcrosser (Colautti et al., 2010; Darwin, 1877; 420 O’Neil, 1997) and therefore individuals in southern populations likely experience strong inter-421 specific competition for pollinators. Additionally, several effective biological control insects 422 have been released into eastern North America, including Gallerucella beetles that feed on 423 meristem tissue as larvae and can drastically reduce flower production (Blossey & Notzold, 424 1995; Grevstad, 2006). The later average flowering and more negative skew in plants from 425 southern populations is therefore consistent with selection from higher herbivory and more 426 competition for pollinators under a longer growing season. We cannot easily parse the specific 427 effects of biotic and abiotic factors based on clines alone, but we can make the prediction that 428 clines in herbivore and pollinator activity would be more patchy and less predictable by latitude 429 compared to season length. It is therefore interesting that we detected strong clines in average 430 start (β = -0.856) and mean (β = -0.881) flowering dates, which we predict to be more sensitive 431 to season length, in contrast to kurtosis (β = -0.590), which we predicted as a bet-hedging 432 strategy for biotic interactions. Skew shows an intermediate pattern (β = 0.745) and could be an 433 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 20 adaptation to biotic selection but could also be partly an artifact of our mid-latitude common 434 garden site, as discussed below. 435 Characteristics of aggregate flowering schedules were generally very different from their 436 population averages, resulting in high variability in emergent properties among populations (Fig. 437 4, right column), but these differences did not generally correlate with latitude. Two important 438 exceptions are the significant and marginally significant clines observed in the onset and 439 duration of flowering respectively, which characterized more synchronized flowering schedules 440 in the north. However, the mid-latitude location of our common garden may also be responsible 441 for these trends. Cooler temperatures that are expected with increasing latitude not only delay the 442 start of favorable growing conditions, but they also restrict insect pollinator activity. Therefore, 443 lack of pollinator activity towards the end of the growing season at our mid-latitude field site 444 may represent an artificial truncation of flowering by southern populations leading to a perceived 445 difference in flowering schedule synchrony. Overall, variation in emergent properties among L. 446 salicaria populations has important implications for field research on components of floral 447 display. It has previously been established that individuals with flowering schedules 448 uncharacteristic of their populations can potentially bias ecological and evolutionary inferences 449 drawn from aggregate population schedules (Newstrom et al. 1994; B. D. Inouye et al., 2019). 450 The emergent properties in our study suggest that such cases are likely more common than is 451 generally appreciated and that they may also be geographically structured. As such, aggregate 452 schedules should be interpreted cautiously. 453 454 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 21

-- 455 Our study demonstrates that populations of L. salicaria in eastern North America have evolved 456 differences in flowering schedules within less than a century of invasion. This phenological 457 differentiation among populations includes predictable changes in the timing of resource 458 allocation to flower development and floral display and adds to the growing list of phenotypic 459 traits that are targets of selection as invasive species spread along environmental gradients 460 (MacDougall et al., 2006; McGoey et al., 2020; Molina-Montenegro et al., 2018; Urbanski et al., 461 2012). Based on the observed clines in our study, we predict that longer growing seasons 462 accompanying anthropogenic climate change will alter flowering schedules. These changes are 463 likely to involve a later start date to flowering, a later mean day of flowering, and a larger 464 kurtosis, which in turn should increase seed production and promote population growth and 465 range expansion of L. salicaria. The ecological and evolutionary insights drawn from 466 observations of flowering start dates and population aggregate characteristics should, however, 467 be interpreted cautiously. In future, more detailed features of individual flowering schedules can 468 be characterised by multivariate methods and metrics drawn from central moments theory. As we 469 have shown for L. salicaria from eastern North America, these more refined characteristics of 470 individual flowering schedules might help in developing a deeper understanding of the diverse 471 ways that species respond to changing environments. 472 473

474 The authors thank Z. Burivalova and R. Mackenzie for field assistance, and to E.J. Austen, A.E. 475 Weis and J. Friedman for feedback on the analysis. Funding was provided from NSERC 476 Discovery grants to S.C.H.B. and R.I.C. and to M.N.A. from NSERC CGSM. 477 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 22 478 AUTHOR CONTRIBUTIONS 479 M.N.A. led the data analysis with assistance from D.M., S.C.H.B., and R.I.C. The 480 conceptualization and experimental design were developed by S.C.H.B. and R.I.C., and all 481 authors contributed to the manuscript, based on an original draft written by M.N.A. 482 483 DATA AVAILABILITY STATEMENT 484 Raw data and fully reproducible R code for all figures, tables, and statistical analyses are 485 available in the Dryad database (DOI:) 486 487 LITERATURE CITED 488 Ågren, J., & Schemske, D. W. (2012). Reciprocal transplants demonstrate strong adaptive 489 differentiation of the model organism Arabidopsis thaliana in its native range. New 490 Phytologist, 194(4), 1112–1122. https://doi.org/10.1111/j.1469-8137.2012.04112.x 491 Alexander, J. M., Edwards, P. J., Poll, M., Parks, C. G., & Dietz, H. (2009). Establishment of 492 parallel altitudinal clines in traits of native and introduced forbs. Ecology, 90(3), 493 612–622. https://doi.org/10.1890/08-0453.1 494 Anderson, J. T., & Gezon, Z. J. (2015). Plasticity in functional traits in the context of climate 495 change: A case study of the subalpine forb Boechera stricta (Brassicaceae). Global 496 Change Biology, 21(4), 1689–1703. https://doi.org/10.1111/gcb.12770 497 Anderson, J. T., Inouye, D. W., McKinney, A. M., Colautti, R. I., & Mitchell-Olds, T. (2012). 498 Phenotypic plasticity and adaptive evolution contribute to advancing flowering 499 phenology in response to climate change. Proceedings of the Royal Society B: 500 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 23 Biological Sciences, 279(1743), 3843–3852. 501 https://doi.org/10.1098/rspb.2012.1051 502 Augspurger, C. K. (1981). Reproductive synchrony of a tropical shrub: Experimental studies 503 on effects of pollinators and seed predators in Hybanthus prunifolius (violaceae). 504 Ecology, 62(3), 775–788. https://doi.org/10.2307/1937745 505 Austen, E. J., Jackson, D. A., & Weis, A. E. (2014). Describing flowering schedule shape 506 through multivariate ordination. International Journal of Plant Sciences. 507 https://doi.org/10.1086/673934 508 Austen, E. J., Rowe, L., Stinchcombe, J. R., & Forrest, J. R. K. (2017). Explaining the 509 apparent paradox of persistent selection for early flowering. New Phytologist, 510 215(3), 929–934. https://doi.org/10.1111/nph.14580 511 Austen, E. J., & Weis, A. E. (2015). What drives selection on flowering time? An 512 experimental manipulation of the inherent correlation between genotype and 513 environment: causes of selection on flowering time. Evolution, 69(8), 2018–2033. 514 https://doi.org/10.1111/evo.12709 515 Blionis, G. J., Halley, J. M., & Vokou, D. (2001). Flowering phenology of Campanula on Mt 516 Olympos, Greece. Ecography, 24(6), 696–706. https://doi.org/10.1111/j.1600-517 0587.2001.tb00531.x 518 Blossey, B., & Notzold, R. (1995). Evolution of Increased Competitive Ability in Invasive 519 Nonindigenous Plants: A Hypothesis. Journal of Ecology, 83(5), 887–889. 520 https://doi.org/10.2307/2261425 521 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 24 CaraDonna, P. J., Iler, A. M., & Inouye, D. W. (2014). Shifts in flowering phenology reshape a 522 subalpine plant community. Proceedings of the National Academy of Sciences, 523 111(13), 4916–4921. https://doi.org/10.1073/pnas.1323073111 524 Colautti, R. I., & Barrett, S. C. H. (2010). Natural selection and genetic constraints on 525 flowering phenology in an invasive plant. International Journal of Plant Sciences, 526 171(9), 960–971. https://doi.org/10.1086/656444 527 Colautti, R. I., & Barrett, S. C. H. (2011). Population divergence along lines of genetic 528 variance and covariance in the invasive plant Lythrum salicaria in eastern North 529 America. Evolution, 65(9), 2514–2529. https://doi.org/10.1111/j.1558-530 5646.2011.01313.x 531 Colautti, R. I., & Barrett, S. C. H. (2013). Rapid adaptation to climate facilitates range 532 expansion of an invasive plant. Science, 342(6156), 364–366. 533 https://doi.org/10.1126/science.1242121 534 Colautti, R. I., & Lau, J. A. (2015). Contemporary evolution during invasion: Evidence for 535 differentiation, natural selection, and local adaptation. Molecular Ecology, 24(9), 536 1999–2017. https://doi.org/10.1111/mec.13162 537 Colautti, R. I., White, N. A., & Barrett, S. C. H. (2010). Variation of Self‐Incompatibility 538 within Invasive Populations of Purple Loosestrife ( Lythrum salicaria L.) from Eastern 539 North America. International Journal of Plant Sciences, 171(2), 158–166. 540 https://doi.org/10.1086/649023 541 Darwin, C. (1877). The Different Forms of Flowers on Plants of the Same Species. John 542 Murray, Albermarle Street. 543 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 25 de Villemereuil, P., Charmantier, A., Arlt, D., Bize, P., Brekke, P., Brouwer, L., Cockburn, A., 544 Côté, S. D., Dobson, F. S., Evans, S. R., Festa-Bianchet, M., Gamelon, M., Hamel, S., 545 Hegelbach, J., Jerstad, K., Kempenaers, B., Kruuk, L. E. B., Kumpula, J., Kvalnes, T., 546 … Chevin, L.-M. (2020). Fluctuating optimum and temporally variable selection on 547 breeding date in birds and mammals. Proceedings of the National Academy of 548 Sciences, 117(50), 31969–31978. https://doi.org/10.1073/pnas.2009003117 549 Elzinga, J. A., Atlan, A., Biere, A., Gigord, L., Weis, A. E., & Bernasconi, G. (2007). Time after 550 time: Flowering phenology and biotic interactions. Trends in Ecology & Evolution, 551 22(8), 432–439. https://doi.org/10.1016/j.tree.2007.05.006 552 Ensing, D. J., & Eckert, C. G. (2019). Interannual variation in season length is linked to 553 strong co-gradient plasticity of phenology in a montane annual plant. New 554 Phytologist, 224(3), 1184–1200. https://doi.org/10.1111/nph.16009 555 Fantinato, E., Del Vecchio, S., Giovanetti, M., Acosta, A. T. R., & Buffa, G. (2018). New 556 insights into plants co-existence in species-rich communities: The pollination 557 interaction perspective. Journal of Vegetation Science, 29(1), 6–14. 558 https://doi.org/10.1111/jvs.12592 559 Fenster, C. B., Armbruster, W. S., Wilson, P., Dudash, M. R., & Thomson, J. D. (2004). 560 Pollination syndromes and floral specialization. Annual Review of Ecology, 561 Evolution, and Systematics, 35, 375–403. 562 Fitchett, J. M., Grab, S. W., & Thompson, D. I. (2015). Plant phenology and climate change: 563 Progress in methodological approaches and application. Progress in Physical 564 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 26 Geography: Earth and Environment, 39(4), 460–482. 565 https://doi.org/10.1177/0309133315578940 566 Forrest, J., & Thomson, J. D. (2010). Consequences of variation in flowering time within and 567 among individuals of Mertensia fusiformis (Boraginaceae), an early spring 568 wildflower. American Journal of Botany, 97(1), 38–48. 569 https://doi.org/10.3732/ajb.0900083 570 Fox, G. A. (2003). Assortative mating and plant phenology: Evolutionary and practical 571 consequences. Evolutionary Ecology Research, 5(1), 18. 572 Gauzere, J., Teuf, B., Davi, H., Chevin, L.-M., Caignard, T., Leys, B., Delzon, S., Ronce, O., & 573 Chuine, I. (2020). Where is the optimum? Predicting the variation of selection along 574 climatic gradients and the adaptive value of plasticity. A case study on tree 575 phenology. Evolution Letters, 4(2), 109–123. https://doi.org/10.1002/evl3.160 576 Grevstad, F. S. (2006). Ten-year impacts of the biological control agents Galerucella pusilla 577 and G. calmariensis (Coleoptera: Chrysomelidae) on purple loosestrife (Lythrum 578 salicaria) in Central New York State. Biological Control, 39(1), 1–8. 579 https://doi.org/10.1016/j.biocontrol.2006.03.007 580 Griffith, T. M., & Watson, M. A. (2005). Stress avoidance in a common annual: Reproductive 581 timing is important for local adaptation and geographic distribution. Journal of 582 Evolutionary Biology, 18(6), 1601–1612. https://doi.org/10.1111/j.1420-583 9101.2005.01021.x 584 Halbritter, A. H., Fior, S., Keller, I., Billeter, R., Edwards, P. J., Holderegger, R., Karrenberg, 585 S., Pluess, A. R., Widmer, A., & Alexander, J. M. (2018). Trait differentiation and 586 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 27 adaptation of plants along elevation gradients. Journal of Evolutionary Biology, 587 31(6), 784–800. https://doi.org/10.1111/jeb.13262 588 Inouye, B. D., Ehrlén, J., & Underwood, N. (2019). Phenology as a process rather than an 589 event: From individual reaction norms to community metrics. Ecological 590 Monographs, 89(2), e01352. https://doi.org/10.1002/ecm.1352 591 Inouye, D. W. (2022). Climate change and phenology. WIREs Climate Change, 13(3), e764. 592 https://doi.org/10.1002/wcc.764 593 Ison, J. L., & Weis, A. E. (2017). Temporal population genetic structure in the pollen pool for 594 flowering time: A field experiment with Brassica rapa (Brassicaceae). American 595 Journal of Botany, 104(10), 1569–1580. https://doi.org/10.3732/ajb.1700210 596 Kirkpatrick, M., & Barton, N. H. (1997). Evolution of a species’ range. The American 597 Naturalist, 150(1), 1–23. https://doi.org/10.1086/286054 598 Kudo, G. (2007). Flowering phenologies of animal-pollinated plants: Reproductive 599 strategies and agents of selection. In Ecology and Evolution of Flowers. Oxford 600 University Press, Incorporated. 601 Legendre, L., & Legendre, P. (1983). Numerical Ecology. Elsevier Scientific Publishing 602 Company. 603 MacDougall, A. s., Boucher, J., Turkington, R., & Bradfield, G. e. (2006). Patterns of plant 604 invasion along an environmental stress gradient. Journal of Vegetation Science, 605 17(1), 47–56. https://doi.org/10.1111/j.1654-1103.2006.tb02422.x 606 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 28 McGoey, B. V., Hodgins, K. A., & Stinchcombe, J. R. (2020). Parallel flowering time clines in 607 native and introduced ragweed populations are likely due to adaptation. Ecology 608 and Evolution, 10(11), 4595–4608. https://doi.org/10.1002/ece3.6163 609 Molina-Montenegro, M. A., Acuña-Rodríguez, I. S., Flores, T. S. M., Hereme, R., Lafon, A., 610 Atala, C., & Torres-Díaz, C. (2018). Is the Success of Plant Invasions the Result of 611 Rapid Adaptive Evolution in Seed Traits? Evidence from a Latitudinal Rainfall 612 Gradient. Frontiers in Plant Science, 9, 208. 613 https://doi.org/10.3389/fpls.2018.00208 614 Montague, J. L., Barrett, S. C. H., & Eckert, C. G. (2008). Re-establishment of clinal 615 variation in flowering time among introduced populations of purple loosestrife 616 (Lythrum salicaria, Lythraceae). Journal of Evolutionary Biology, 21(1), 234–245. 617 https://doi.org/10.1111/j.1420-9101.2007.01456.x 618 Munguía-Rosas, M. A., Ollerton, J., Parra-Tabla, V., & De-Nova, J. A. (2011). Meta-analysis 619 of phenotypic selection on flowering phenology suggests that early flowering plants 620 are favoured. Ecology Letters, 14(5), 511–521. https://doi.org/10.1111/j.1461-621 0248.2011.01601.x 622 Newstrom, L. E., Frankie, G. W., & Baker, H. G. (1994). A new classification for plant 623 phenology based on flowering patterns in lowland tropical rain forest trees at Le 624 Selva, Costa Rica. Biotropica, 26(2) 151-169. 625 Olsson, K., & Ågren, J. (2002). Latitudinal population differentiation in phenology, life 626 history and flower morphology in the perennial herb Lythrum salicaria. Journal of 627 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 29 Evolutionary Biology, 15(6), 983–996. https://doi.org/10.1046/j.1420-628 9101.2002.00457.x 629 O’Neil, P. (1997). Natural selection on genetically correlated phenological characters in 630 Lythrum salicaria l. (Lythraceae). Evolution, 51(1), 267–274. 631 https://doi.org/10.2307/2410980 632 Panchen, Z. A. (2022). Plant reproductive phenology along an elevation gradient in the 633 extreme environment of the Canadian High Arctic. Plant Ecology & Diversity, 15(5–634 6), 213–226. https://doi.org/10.1080/17550874.2022.2147804 635 Panchen, Z. A., & Gorelick, R. (2017). Prediction of Arctic plant phenological sensitivity to 636 climate change from historical records. Ecology and Evolution, 7(5), 1325–1338. 637 https://doi.org/10.1002/ece3.2702 638 Prather, R. M., Dalton, R. M., barr, billy, Blumstein, D. T., Boggs, C. L., Brody, A. K., Inouye, 639 D. W., Irwin, R. E., Martin, J. G. A., Smith, R. J., Van Vuren, D. H., Wells, C. P., 640 Whiteman, H. H., Inouye, B. D., & Underwood, N. (2023). Current and lagged climate 641 affects phenology across diverse taxonomic groups. Proceedings of the Royal 642 Society B: Biological Sciences, 290(1990), 20222181. 643 https://doi.org/10.1098/rspb.2022.2181 644 Price, T., Kirkpatrick, M., & Arnold, S. J. (1988). Directional selection and the evolution of 645 breeding date in birds. Science, 240(4853), 798–799. 646 https://doi.org/10.1126/science.3363360 647 R Core Team. (2022). R: A language and environment for statistical computing. [Computer 648 software]. R Foundation for Statistical Computing. https://www.R-project.org/ 649 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 30 Rabinowitz, D., Rapp, J. K., Sork, V. L., Rathcke, B. J., Reese, G. A., & Weaver, J. C. (1981). 650 Phenological properties of wind- and insect-pollinated prairie plants. Ecology, 651 62(1), 49–56. https://doi.org/10.2307/1936667 652 Rafferty, N. E., Diez, J. M., & Bertelsen, C. D. (2020). Changing climate drives divergent and 653 nonlinear shifts in flowering phenology across elevations. Current Biology, 30(3), 654 432-441.e3. https://doi.org/10.1016/j.cub.2019.11.071 655 Rathcke, B., & Lacey, E. P. (1985). Phenological patterns of terrestrial plants. Annual 656 Review of Ecology and Systematics, 16(1), 179–214. 657 https://doi.org/10.1146/annurev.es.16.110185.001143 658 Rausher, M. D. (1992). The measurement of selection on quantitative traits: Biases due to 659 environmental covariances between traits and fitness. Evolution, 46(3), 616–626. 660 https://doi.org/10.1111/j.1558-5646.1992.tb02070.x 661 RStudio Team. (2022). RStudio: Integrated Development for R [Computer software]. 662 RStudio, PBC. http://www.rstudio.com/. 663 Samis, K. E., Murren, C. J., Bossdorf, O., Donohue, K., Fenster, C. B., Malmberg, R. L., 664 Purugganan, M. D., & Stinchcombe, J. R. (2012). Longitudinal trends in climate drive 665 flowering time clines in North American Arabidopsis thaliana. Ecology and 666 Evolution, 2(6), 1162–1180. https://doi.org/10.1002/ece3.262 667 Simons, A. M. (2011). Modes of response to environmental change and the elusive 668 empirical evidence for bet hedging. Proceedings of the Royal Society B: Biological 669 Sciences, 278(1712), 1601–1609. https://doi.org/10.1098/rspb.2011.0176 670 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 31 Sobral, M., Guitián, J., Guitián, P., & Larrinaga, A. R. (2013). Selective pressure along a 671 latitudinal gradient affects subindividual variation in plants. PLOS ONE, 8(9), 672 e74356. https://doi.org/10.1371/journal.pone.0074356 673 Sokal, R. R., & Rohlf, F. J. (1995). Biometry (3rd ed.). WH Freeman. 674 Stemkovski, M., Dickson, R. G., Griffin, S. R., Inouye, B. D., Inouye, D. W., Pardee, G. L., 675 Underwood, N., & Irwin, R. E. (2023). Skewness in bee and flower phenological 676 distributions. Ecology, 104(1), e3890. https://doi.org/10.1002/ecy.3890 677 Strauss, S. Y., & Whittall, J. B. (2006). Non-pollinator agents of selection on floral traits. In 678 Ecology and Evolution of Flowers (pp. 120–135). Oxford Biology. 679 Thomson, J. D. (1980). Skewed flowering distributions and pollinator attraction. Ecology, 680 61(3), 572–579. https://doi.org/10.2307/1937423 681 Tufto, J. (2015). Genetic evolution, plasticity, and bet-hedging as adaptive responses to 682 temporally autocorrelated fluctuating selection: A quantitative genetic model. 683 Evolution, 69(8), 2034–2049. https://doi.org/10.1111/evo.12716 684 Urbanski, J., Mogi, M., O’Donnell, D., DeCotiis, M., Toma, T., & Armbruster, P. (2012). Rapid 685 Adaptive Evolution of Photoperiodic Response during Invasion and Range 686 Expansion across a Climatic Gradient. The American Naturalist, 179(4), 490–500. 687 https://doi.org/10.1086/664709 688 van Kleunen, M., Bossdorf, O., & Dawson, W. (2018). The ecology and evolution of alien 689 plants. Annual Review of Ecology, Evolution, and Systematics, 49(1), 25–47. 690 https://doi.org/10.1146/annurev-ecolsys-110617-062654 691 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 32 Waites, A. R., & Ågren, J. (2006). Stigma receptivity and effects of prior self-pollination on 692 seed set in tristylous Lythrum salicaria (Lythraceae). American Journal of Botany, 693 93(1), 142–147. https://doi.org/10.3732/ajb.93.1.142 694 Wolkovich, E. M., Cook, B. I., Allen, J. M., Crimmins, T. M., Betancourt, J. L., Travers, S. E., 695 Pau, S., Regetz, J., Davies, T. J., Kraft, N. J. B., Ault, T. R., Bolmgren, K., Mazer, S. J., 696 McCabe, G. J., McGill, B. J., Parmesan, C., Salamin, N., Schwartz, M. D., & Cleland, 697 E. E. (2012). Warming experiments underpredict plant phenological responses to 698 climate change. Nature, 485(7399), Article 7399. 699 https://doi.org/10.1038/nature11014 700 Wright, J. W., & Meagher, T. R. (2003). Pollination and seed predation drive flowering 701 phenology in Silene latifolia (Caryophyllaceae). Ecology, 84(8), 2062–2073. 702 https://doi.org/10.1890/02-0676 703 Wu, Y., & Colautti, R. I. (2022). Evidence for continent-wide convergent evolution and 704 stasis throughout 150 y of a biological invasion. Proceedings of the National 705 Academy of Sciences, 119(18), e2107584119. 706 https://doi.org/10.1073/pnas.2107584119 707 Zvereva, E. L., & Kozlov, M. V. (2021). Latitudinal gradient in the intensity of biotic 708 interactions in terrestrial ecosystems: Sources of variation and differences from the 709 diversity gradient revealed by meta-analysis. Ecology Letters, 24(11), 2506–2520. 710 https://doi.org/10.1111/ele.13851 711 712 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 33 Box 1: Summary characteristics of flowering schedules Central moments (i.e., mean, variance, skew, and kurtosis) arise from probability theory and characterize distributions that can be visualized as histograms with observed values on the x-axis and probability density, frequency, or number of observations on the y-axis. Here, we use the same equations to describe analogous characteristics of flowering schedules that arise from developmental rather than probabilistic process. In contrast to probability histograms, flowering schedules can be visualized as a time series with number or proportion of flowers on the y-axis. Although the underlying processes are distinct, the equations are the same (Table 1). To understand the biological significance of mean, variance, skew and kurtosis, we contrast flowering schedules of different shape (Fig. 1). Standardizing to proportion of open flowers (P t) over time (t), rather than total flower number (Nt) accounts for variation in total flower number. Moreover, we can use proportions to calculate a weighted mean day of flowering. By analogy to the mean of a probability distribution, the weighted mean of a flowering schedule represents the ‘balance point’, which better captures the anthesis day of a typical flower. As demonstrated in Figure 1A, central moments capture biologically meaningful variation in flowering schedules that may not be correlated with onset or duration of flowering. For example, the variance parameter (𝜎𝜎 𝑖𝑖 2) describes how much flowering is spread out, even if the start and end dates are the same. An individual’s coefficient of skew (CSi) accounts for expected differences in variance and thus describes whether flowering is concentrated earlier (i.e., right skew, CSi > 0 as shown by the blue dashed line) or later in the schedule (i.e., left skew, CSi < 0 as shown by the purple dotted line). Finally, the coefficient of kurtosis (CKi) describes deviations from the Gaussian expectation, contrasting individuals with a concentrated cluster of open flowers and a few flowers over a long-time frame (i.e., leptokurtic distribution, CKi > 0 as shown in the gold dashed line) or individuals that spread out flowering more evenly over the entire schedule (i.e., platykurtic, CKi < 0 as shown by the dashed green line). In addition to calculating these characteristics for individual plants, all flowers (or proportion of flowers) on each day can be pooled to yield a population aggregate flowering schedule (Fig. 1B). The aggregate schedule of a population may not be representative of the individual schedules comprising it, resulting in emergent properties that we define as deviations from the average individual flowering schedule. .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 34 Table 1. Equations of flowering schedule characteristics for an average individual within each population and for the aggregate population of pooled individuals. Summary Statistic Individual (i) Population aggregate (j) Proportion (𝑃𝑃) of flowers (𝐹𝐹) produced at time 𝑡𝑡 𝑃𝑃𝑖𝑖,𝑡𝑡= 𝐹𝐹𝑖𝑖,𝑡𝑡 ∑ 𝐹𝐹𝑖𝑖,𝑡𝑡 𝑡𝑡 𝑃𝑃𝑡𝑡,𝑗𝑗 ∑ = 𝐹𝐹𝑡𝑡,𝑗𝑗 ∑ ∑ 𝐹𝐹𝑡𝑡,𝑗𝑗 ∑ 𝑡𝑡 = ∑ 𝐹𝐹𝑖𝑖,𝑡𝑡,𝑗𝑗 𝑖𝑖 ∑ ∑ 𝐹𝐹𝑖𝑖,𝑡𝑡,𝑗𝑗 𝑖𝑖 𝑡𝑡 Day of first flower 𝛼𝛼𝑖𝑖= min{𝑡𝑡 ∣ 𝐹𝐹𝑖𝑖,𝑡𝑡> 0} 𝛼𝛼𝑗𝑗 ∑ = min{𝑡𝑡 ∣ � 𝐹𝐹𝑖𝑖,𝑡𝑡,𝑗𝑗 𝑖𝑖 > 0} Day of last flower 𝛺𝛺𝑖𝑖= max{𝑡𝑡 ∣ 𝐹𝐹𝑖𝑖,𝑡𝑡> 0} 𝛺𝛺𝑗𝑗 ∑ = max{𝑡𝑡 ∣ � 𝐹𝐹𝑖𝑖,𝑡𝑡,𝑗𝑗 𝑖𝑖 > 0} Duration of flowering 𝛺𝛺𝑖𝑖− 𝛼𝛼𝑖𝑖 𝛺𝛺𝑗𝑗 ∑ − 𝛼𝛼𝑗𝑗 ∑ Day of mean flowering 𝜇𝜇𝑖𝑖= � 𝑡𝑡 𝑡𝑡 𝑃𝑃𝑖𝑖,𝑡𝑡 𝜇𝜇𝑗𝑗 ∑ = � 𝑡𝑡 𝑡𝑡 𝑃𝑃𝑡𝑡,𝑗𝑗 ∑ Flowering schedule variance 𝜎𝜎𝑖𝑖 2 = � (𝑡𝑡 − 𝜇𝜇)2 𝑡𝑡 𝑃𝑃𝑡𝑡 𝜎𝜎𝑗𝑗 2∑ = � �𝑡𝑡 − 𝜇𝜇𝑗𝑗 ∑� 2 𝑡𝑡 𝑃𝑃𝑡𝑡,𝑗𝑗 ∑ Flowering schedule coefficient of skew 𝐶𝐶𝑆𝑆𝑖𝑖= ∑ (𝑡𝑡 − 𝜇𝜇𝑖𝑖)3 𝑡𝑡 𝑃𝑃𝑖𝑖,𝑡𝑡 𝜎𝜎𝑖𝑖 3 𝐶𝐶𝑆𝑆𝑗𝑗= ∑ �𝑡𝑡 − 𝜇𝜇𝑗𝑗 ∑� 3 𝑡𝑡 𝑃𝑃𝑡𝑡,𝑗𝑗 ∑ 𝜎𝜎𝑗𝑗 3∑ Flowering schedule coefficient of kurtosis 𝐶𝐶𝐾𝐾𝑖𝑖= ∑ (𝑡𝑡 − 𝜇𝜇𝑖𝑖)4 𝑡𝑡 𝑃𝑃𝑖𝑖,𝑡𝑡 𝜎𝜎𝑖𝑖 4 − 3 𝐶𝐶𝐾𝐾𝑗𝑗= ∑ �𝑡𝑡 − 𝜇𝜇𝑗𝑗 ∑� 4 𝑡𝑡 𝑃𝑃𝑡𝑡,𝑗𝑗 ∑ 𝜎𝜎𝑗𝑗 4∑ .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 35 Table 2. Summary of the first five eigenvector axes of a Principal Coordinates Analysis (PCoA) using Kolmogorov-Smirnov distances between flowering schedules of 369 plants from 13 populations of Lythrum salicaria sampled along a latitudinal gradient and grown in a common garden at the Koffler Scientific Reserve in Newmarket, ON. Percent variation is calculated as the eigenvalue over the sum of all positive eigenvalues. Axis Eigenvalue Variation explained (%) PCoA1 13.4 23.5 PCoA2 6.6 11.5 PCoA3 3.9 6.9 PCoA4 3.2 5.6 PCoA5 2.6 4.7 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 36 Table 3. Standardized slope and test statistics of univariate linear regressions testing latitude as a predictor of the bootstrapped means of average and aggregate population flowering schedule traits, and the emergent property (i.e. population aggregate minus average) in populations of Lythrum salicaria. Bold text indicates significant values. Linear term Quadratic term Traits β F (1,11) P X2 P Start Average -0.856 30.0 < 0.001 1.73 0.189 Aggregate -0.591 5.92 < 0.05 0.01 0.907 Emergent 0.857 30.52 < 0.001 2.50 0.114 Duration Average -0.369 1.73 0.215 0.11 0.738 Aggregate -0.849 28.37 < 0.001 0.01 0.904 Emergent -0.470 3.13 0.105 2.96 0.085 Mean Average -0.881 38.0 < 0.001 0.47 0.494 Aggregate -0.845 27.5 < 0.001 0.12 0.724 Emergent 0.368 1.66 0.224 0 0.948 Variance Average 0.031 0.01 0.920 0.07 0.791 Aggregate -0.214 0.53 0.483 0 0.977 Emergent 0.156 0.276 0.610 6.83 < 0.01 Skew Average 0.745 13.7 < 0.01 1.03 0.310 Aggregate 0.382 1.87 0.198 0.07 0.797 Emergent 0.227 0.60 0.456 0.12 0.741 Kurtosis Average -0.590 5.88 < 0.05 0.20 0.659 Aggregate -0.592 5.94 < 0.05 0.88 0.347 Emergent -0.133 0.19 0.664 1.18 0.279 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 37 Figure 1. Hypothetical flowering schedules of (A) five individuals with different characteristics and (B) their pooled aggregate flowering schedule. Flowering schedules are scaled to the proportion of total flower production (P t) over time (t). .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 38 Figure 2. Matrices showing pairwise Pearson correlation coefficients among the start, duration, central moments (i.e., variance, skew and kurtosis), and two the first two eigenvectors of a principal coordinates analysis (PCoA) characterizing flowering schedules in populations of Lythrum salicaria grown in a common garden. The metrics were calculated using either (A) flowering schedule characteristics averaged across individuals within each population (n = 13), (B) the residual correlations calculated on deviations of individuals from their population means (n = 369), or (C) aggregated schedules for individuals pooled within populations (n = 13). Circles are shown for statistically significant correlations (p < 0.05), colour-coded by correlation coefficient (i.e., -1 to 1 as shown on the x-axis). .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 39 Figure 3. Principal Coordinates Analysis (PCoA) of Kolmogorov-Smirnov distances between flowering schedules of 369 plants from 13 populations of Lythrum salicaria sampled along a latitudinal gradient and grown in a common garden at the Koffler Scientific Reserve in Newmarket, ON. Point colors correspond to latitude of origin, ranging from dark (48° N) to light (38° N). Bar plots show representative flowering schedules of four individual plants. PCoA1 is inversely related to the central moments of variance and mean such that points on the right represent individuals with lower flowering schedule variance and an earlier peak. .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 40 Figure 4. Bivariate plots showing relations between latitude and flowering schedule characteristics of 13 populations of Lythrum salicaria grown in a common garden field experiment. The bootstrapped estimate and 95% confidence intervals of average individual flowering schedule characteristics are shown for each population along the left column. The bootstrapped estimate and 95% confidence intervals of the aggregate flowering schedules for all flowers within each population are shown in the middle column. The right column indicates the ‘emergent properties’ of each population, defined as the aggregate minus the average. .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 41 Appendix S1. The distribution of 369 individuals and 216 seed families representing 13 populations of Lythrum salicaria sampled along a latitudinal gradient and grown in a common garden at the Koffler Scientific Reserve in Newmarket, ON. Point colors correspond to latitude of origin ranging from dark (48° N) to light (38° N). Latitude (°) Number of seed families Number of individuals 38.75261 15 23 40.34025 15 33 41.31283 18 24 42.23822 20 34 42.52306 12 34 43.69719 19 21 44.49339 19 38 44.59489 19 22 45.26278 16 26 45.43194 17 36 45.49000 16 23 47.69081 14 32 48.47892 16 23 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint 42 FIGURE LEGENDS Figure 1. Hypothetical flowering schedules of (A) five individuals with different characteristics and (B) their pooled aggregate flowering schedule. Flowering schedules are scaled to the proportion of total flower production (Pt) over time (t). Figure 2. Matrices showing pairwise Pearson correlation coefficients among the start, duration, central moments (i.e., variance, skew and kurtosis), and two the first two eigenvectors of a principal coordinates analysis (PCoA) characterizing flowering schedules in populations of Lythrum salicaria grown in a common garden. The metrics were calculated using either (A) flowering schedule characteristics averaged across individuals within each population (n = 13), (B) the residual correlations calculated on deviations of individuals from their population means (n = 369), or (C) aggregated schedules for individuals pooled within populations (n = 13). Circles are shown for statistically significant correlations (p < 0.05), colour-coded by correlation coefficient (i.e., -1 to 1 as shown on the x-axis). Figure 3. Principal Coordinates Analysis (PCoA) of Kolmogorov-Smirnov distances between flowering schedules of 369 plants from 13 populations of Lythrum salicaria sampled along a latitudinal gradient and grown in a common garden at the Koffler Scientific Reserve in Newmarket, ON. Point colors correspond to latitude of origin, ranging from dark (48° N) to light (38° N). Bar plots show representative flowering schedules of four individual plants. PCoA1 is inversely related to the central moments of variance and mean such that points on the right represent individuals with lower flowering schedule variance and an earlier peak. Figure 4. Bivariate plots showing relations between latitude and flowering schedule characteristics of 13 populations of Lythrum salicaria grown in a common garden field experiment. The bootstrapped estimate and 95% confidence intervals of average individual flowering schedule characteristics are shown for each population along the left column. The bootstrapped estimate and 95% confidence intervals of the aggregate flowering schedules for all flowers within each population are shown in the middle column. The right column indicates the ‘emergent properties’ of each population, defined as the aggregate minus the average. .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted June 3, 2024. ; https://doi.org/10.1101/2024.05.30.596697doi: bioRxiv preprint