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).