Materials and methods

Study system

Lathyrus vernus (L. Bernh) is a perennial herb that is distributed over Europe and North-West Asia (POWO 2019). In the study area, L. vernus is mainly found in the understorey of deciduous or mixed-deciduous forests. Each plant individual produces one to several shoots that emerge early in spring and that can reach up to 40 cm in height (Ehrlén 2002). The leaves are pinnate and consist of 2-4 pairs of leaflets that unfold starting with the basalmost leaf on the shoot and the pair of leaflets closest to the stem. Each flowering individual produces 1 to 5 racemes, each having 1 to 9 flowers with the basalmost flowers opening first (Ehrlén 2002). The racemes emerge from leaf axes and start development when the neighbouring leaf unfolds. Racemes and their adjacent leaf develop simultaneously, and although flowers normally open after that leaflets have unfolded, the first flower might sometimes open before the first leaflet has unfolded. A previous study with this species showed that plants with an earlier leaf development initiate flowering earlier than plants with a later leaf development (Sola and Ehrlén 2007). In South-East Sweden, where this study was carried out, flowering usually starts in late April to early May. Fruits contain up to 18 ovules, and the mature seeds are dispersed ballistically about 2 months after flowering (Ehrlén 1992). L. vernus reproduces only sexually.
The reproductive performance of L. vernus is influenced by several abiotic and biotic factors, some of which have been linked to selection on flowering time. Cold April temperatures are associated with weaker selection for early flowering, possibly because frost damage reduces the benefit of flowering early (Ehrlén and Valdés 2020). Grazing by roe deer (Capreolus capreolus ) favours later flowering because early-flowering individuals experience the highest levels of damage (Fogelström and Ehrlén 2019). L. vernus individuals rely on pollinators (Bombus sp.) for their reproduction and seed production is sometimes pollen limited (Ehrlén 1992), suggesting that there are fitness benefits of synchronising flowering with a high availability of bumblebees. Bruchus atomarius (Chrysomelidae) larvae can damage a large proportion of the seeds before dispersal in the study area. B. atomarius seed predation is sometimes linked to the timing of flowering in L. vernus, but the relationship is generally weak and its direction varies among years (Ehrlén 1996, Ehrlén and Münzbergová 2009).

Data collection

This study was conducted in a mixed deciduous forest at Kålsö, Sweden (58°56’N, 17°39’E), in the years 2013 to 2015. We recorded all individuals in the population from shoots emergence in spring each year. We only included individuals that had not been damaged (e.g. by trampling, falling tree branches or by molluscs) before the first recording. Just before the first leaves in the population started to unfold, we increased the frequency of recordings to once or twice per week and continued recordings until both leaf-unfolding and flowering were terminated. During the study, individuals started to leaf-out 23 April – 20 May 2013, 2 April – 14 May 2014, and 26 March – 5 May 2015. Individuals initiated flowering 10 – 28 May in 2013, 18 April – 25 May in 2014, and 20 April – 23 May in 2015, respectively. The average time between leaf-out and first flowering day was 10 days in 2013, 18 days in 2014 and 20 days in 2015.
We used the estimated day of year when the first leaflet unfolded as a measure of vegetative spring phenology. Leaflets open sequentially inL. vernus , and we recorded the number of unfolded leaflets at each visit. Each individual was assigned a leaf-out day within the interval between the recording when the first leaflet was observed to be unfolded and the recording prior. We then used the proportion of open leaflets at the first recording to estimate the most probable first day with any open leaflets within that interval. We assumed that individuals with a larger proportion of open leaflets had started to leaf out earlier than individuals with a smaller proportion, and used a linear model including the total number of leaflets and its squared term to predict the number of unfolded leaflets at the first recording (Appendix S1, Table S1). Based on the deviation of the observed proportion of unfolded leaflets from this predicted proportion, we assigned individuals to a most likely leaf-out day within each recording interval, i.e. plants with larger deviations from the predicted proportion of unfolded leaflets were assigned earlier leaf-out dates than individuals with lower values (Appendix S1, Figure S1). We assigned a leaf-out day to all individuals that were undamaged at leaf-out and that had been recorded on at least one occasion before leaf-out.
We used the day of year when the first flower opened (first flowering day) as a measure of flowering phenology. We recorded the size of the largest bud for each flowering individual at each visit up to flowering. The size of the largest bud at the recording before the recording of the first open flower was then used to assign a first flowering day value when first flowering day occurred between two recordings (Appendix S2). First flowering day for individuals grazed before the first flower was observed was estimated using information about the relationship between first flowering day and bud size, day of year for the bud size observation, and aboveground volume of intact individuals (Appendix S2). We used the number of days from of leaf-out to first flowering day as a measure of relative timing. We refer to this timespan as “development time” henceforward.
The total number of shoots and the height and diameter of the focal shoot was measured in early July in each year. We then calculated the size of each individual in terms of the aboveground volume, using the formula for the volume of a cylinder, applying this to a focal shoot, and multiplying the resulting volume with the number of shoots (aboveground volume = (0.5 × shoot diameter)2 × shoot height × π × number of shoots). During each visit in spring, we recorded incidences of roe-deer grazing. The final height of individuals grazed before the last recording in spring was estimated from the shoot height-diameter relationship of non-grazed plants (Appendix S3). For individuals that had been grazed between the last recording in spring and the final shoot recording in early July, we used the maximum recorded height in spring as an estimate of final shoot height.
We used the number of developed seeds that had not been damaged byB. atomarius larvae, as a measure of plant fitness. We collected all fruits in the population in early July each year. The number of seeds in fruits that had not yet opened was counted in the field, and each seed was checked for B. atomarius entry holes. The number of seeds in fruits that had opened prior to collection in the field was estimated at the lab at Stockholm University. We counted the total number of seeds from clearly visible indentations made by the seeds in the pod walls. We also registered the number of B. atomariusentry holes in the pod wall, and estimated the proportion of seeds that had been preyed upon from the relationship between the proportions of seeds preyed upon and the number of B. atomarius entrance holes in the pod wall of intact fruits. We used the relationship between these two variables established in a previous study with L. vernus to estimate seed predation (proportion of seeds preyed upon = 1-e-1.218 B. atomarius entry holes per seed; Supporting Information Appendix S3 in Fogelström and Ehrlén 2019). We were able to get estimates of size, vegetative phenology, reproductive phenology and fitness from 198, 207 and 207 individuals in 2013, 2014 and 2015, respectively.

Statistical analyses

All statistical analyses were carried out in R version 4.0.5 (R Core Team 2021). Analyses were carried out separately for each year because a large proportion of individuals (46.1%) were recorded in only one year, and including individual as a random effect in mixed effects models could render the models unstable (Harrison et al. 2018). Summary statistics for the variables used in the analyses are presented in Appendix S4.
To estimate the strength of the correlation between leaf-out day and first flowering day, we calculated Pearson’s r for the untransformed variables We estimated phenotypic selection on leaf-out and first flowering day using linear models (Lande and Arnold 1983). Before analysis, we standardized leaf-out day, first flowering day, and the development time between leaf-out and first flowering day (\((x-\overset{\overline{}}{x})/\sigma\)) to a mean of 0 and unit standard deviation. Aboveground volume was transformed to its natural logarithm before standardization. Fitness was relativized\((x/\overset{\overline{}}{x})\) to a mean of 1. For soft selection, we expect the selection surface to be related to the flowering time of a focal individual relative to other individuals within a given year, rather than to day number. For hard selection, it is possible that the selective surface is similar across years, e.g. in terms of the probability of cold weather conditions, but that individuals experience different parts of this surface in different years. To examine selection corresponding to these two scenarios, we ran two sets of selection models that were based on relativizing fitness and standardizing traits within and across years, respectively (De Lisle and Svensson 2017). These two sets of models yielded similar results, and we present the results based on local within-year relativizations and standardizations in the main text, and the results for global across years relativizations and standardizations in Appendix S5. A relatively large proportion of the flowering individuals produced no seeds (67% in 2013, 31% in 2014 and 48% in 2015), and thus the residuals from the linear models did not meet the assumption of normality. We therefore evaluated the significance for the estimates from all linear models by estimating 95% bias-corrected and accelerated (BCa) bootstrap intervals using the function ‘Boot’ from the package ‘car’ version 3.0-10 with 10000 replications (Fox and Weisberg 2019) and ‘boot.ci’ in the package ‘boot’ version 1.3-27 (Davison and Hinkley 1997, Canty and Ripley 2021) in R. Estimates for which the 95% BCa interval does not overlap 0 are considered significant at α = 0.05.
To estimate total selection on leaf-out and first flowering day, we estimated linear selection differentials from simple linear regression models for the two traits, with fitness as a response variable and the trait as the predictor variable. We estimated nonlinear selection by adding a squared term to each model. To reduce bias of the selection estimates caused by differences plant condition, we included plant size as a covariate in these regression models (cf. Rausher 1992). Below, we report the selection estimates from these models including size as “total selection”. The results of corresponding models without size are presented in Appendix S6.
To investigate to what extent selection on first flowering day was direct vs. indirect via leaf-out day, we estimated direct selection in terms of linear selection gradients for leaf-out and first flowering day from a multiple linear regression with fitness as a response variable, and aboveground volume as a covariate. We then examined whether these estimates of direct selection differed from the estimates of total selection on leaf-out and first flowering day. We considered estimates to differ if the estimate of direct selection was outside the BCa intervals for the estimate of total selection on each trait. Differences between the estimate of direct selection and the estimate of total selection on leaf-out day or first flowering day were interpreted as that selection was exerted indirectly, via the other phenological trait.
We investigated selection on the relationship between leaf-out and first flowering day in two ways: First, we estimated selection on the development time from leaf-out to first flowering day, using linear models with relative fitness as a response variable, development time as a predictor variable, and plant size as a covariate. Second, we estimated correlational selection by including the interaction leaf-out day × first flowering day in the multiple linear regression model with relative fitness as a response variable, leaf-out and first flowering day as predictor variables and with plant size as a covariate (cf. Brodie 1992). A statistically supported effect of the leaf-out day × first flowering day interaction indicates that there is correlational selection on the two traits, i.e. that the strength or direction of selection on flowering time and timing of leaf-out is dependent on the level of the other trait. Selection on development time, on the other hand, could occur regardless of whether leaf-out and flowering time are correlated or not, and regardless of whether selection on flowering time is dependent on the timing of leaf-out or not.
Lastly, we ran models with leaf-out day2 and first flowering day2 added to the model estimating correlational selection, to estimate the nonlinear selection gradients. All nonlinear selection estimates were doubled to obtain more accurate measures of the magnitude of nonlinear selection (Stinchcombe et al. 2008).