Data sets and R scripts used for the publication "Species' competitive ability, stress tolerance and phenological niche affect flowering curves in herbaceous species"

These data sets and R scripts were used for the analysis of the publication: Species’ competitive ability and phenology niche affects the flowering intensity in herbaceous species. Phenological data for this was collected between 2019 and 2024 in 15 botanical gardens. Data collection and analysis were carried out as part of the PhenObs project (https://www.idiv.de/research/projects/phenobs/). Abstract: Climate change may increase phenological mismatches in biotic interactions due to shifts in flowering times of plants. Most studies focussed on first flowering, the timing of entire flowering periods have hardly been studied though this information is important when evaluating phenological mismatches. Here, we explore variations in flowering curves across 258 perennial herbaceous species spanning the entire flowering period, and determine whether species-specific patterns are linked to species’ functional properties. The results clearly suggest that competitive species, early-flowering and late-flowering species tended to invest resources in single but intensive flowering events, with shorter flowering durations and more left-skewed curves (‘all-in-one’ strategy). In contrast, stress-tolerant species distributed resources over several flowering peaks (‘bet-hedging’ strategy). We conclude that information on species properties can be used to extract information on different flowering strategies, that can be used to evaluated impacts of climate change not only on flowering times but also on biotic interactions. Readme: To work with the data and scripts, please download all files and save them in one folder. Data preparation:To characterise the shape of the annual flowering curves per accession, we extracted four characteristics: Number of flowering peaks ('peaks'), maximum flowering intensity (‘intensity’), flowering duration (‘duration’) and skewness (‘skewness’). We used the findpeaks() function to determine the number of flowering peaks of the curves ('peaks'). We only included peaks that have reached at least one fifth of the maximum flowering intensity. Due to the observation interval of one week, a peak had to be at least seven days away from another peak to be counted as an additional peak. We applied both filtering rules to compensate for fluctuations owed to the smoothing of the curves and to determine the number of peaks more reliably. The highest amplitude of the curve represents the maximum flowering intensity in % (‘intensity’), and the difference between the day of the first and the last open flower is the flowering duration in days (‘duration’). Finally, for all the curves with a single peak, we determined the skewness (‘skewness’) using the skewness() function. To extract the properties ('peaks', ‘intensity’ and ‘skewness’) of the flowering curves, use the scripts:1_Flowering_Curves_NumberPeak_FloweringIntensity.R 1_Flowering_Curves_Skewness.R.Flowering duration is extracted from Phenology_stages_2019_2024.csv. Data analyses: We first calculated the mean number of flowering peaks of each accession across the years within one garden and afterwards the species’ mean across all botanical gardens. For all species we analysed the associations of the species' functional properties with the 'intensity' and 'duration' of the flowering curves, and additionally for the species with a single peak with the 'skewness' of the curves. To get an idea of the variation in the flowering curve parameters ‘intensity’, ‘duration’, and ‘skewness’, we determined the standard deviation for each species across all curves. Since the species were not distributed evenly across the botanical gardens and years, which can influence the species' mean values, we estimated the means using linear mixed-effects models accounting for variation among botanic gardens and years. We created models with phenological niche (expressed by flowering onset), ‘intensity’ and ‘duration’ for all species, and ‘skewness’ only for species with single peaks as response variables, species as fixed factor, and botanic garden and year as random factors. In order to exclude outliers from the analyses, we first removed observations with residuals deviating by more than two standard deviations from the species-specific residual mean. We then used the models mentioned above extracted the estimated marginal means (EEMs), correcting for variations among botanical gardens and years. To analyse associations between the number of peaks per species and the species’ functional properties, we applied logistic regressions (Generalised Linear Model, glm() function) with a binomial distribution to model the probability of a particular group of the categorical (one peak vs. more than one peak) dependent variable in relation to the continuous predictor variables (competitiveness, stress tolerance or phenological niche). To test for possible nonlinear relationships of the phenological niche (expressed by flowering onset), we included this variable in the model both as a linear and as a quadratic predictor (i.e. Phenological niche 2). To assess the model fit, we calculated McFadden's Pseudo-R2, providing a measure of the relative improvement in the model compared to the null model and reflecting the model's explanatory power regarding the observed data. To analyse whether ‘intensity’ and ‘duration’ for all species, and ‘skewness’ only for species with single peaks were related to species’ competitiveness, stress tolerance or phenological niche, we applied three independent polynomial regression models. We considered the calculated species’ means of the ‘intensity’, ‘duration’, and ‘skewness’ as response variables and each with the species’ means across all gardens of competitiveness, stress tolerance and phenological niche as explanatory variables (see Table S1). As with the analysis of the number of peaks, we included the quadratic predictor in the models when analysing the relationships between the phenological niche and the curve parameters. For running the models, we used the glmmTMB() function (Brooks et al. 2017) with a beta distribution for the analysis of ‘intensity’, as these values range between 0 and 1, and a Gaussian distribution for the analysis of ‘duration’ and ‘skewness’. We tested significance with Type II ANOVAs and extracted partial pseudo-R² values by fitting reduced models excluding step-wise predictors and subtracting their R² from the R² of the full model. To investigate whether the relationships between the characteristics of the flowering curves and the competitiveness, stress tolerance or phenological niche of a species were consistent across botanical gardens (Hypothesis 1), we used the accession specific means per species of ‘intensity’, ‘duration’, and ‘skewness’ (hereafter garden mean). We determined these in the same way as before, using LMMs with only year as a random factor. For the following analyses, we restricted ourselves only to the associations that were significant in the previous described quadratic polynomial regression models when comparing the species' means (previous analyses: Stress tolerance vs. ‘skewness’ (1), Competitiveness vs. ‘duration’ (2) and ‘skewness’ (3) and phenological niche vs. ‘intensity’ (4)). For all models we used the glmmTMB() function. For (1)-(3), we considered a gaussian distribution with the garden-means of ‘duration’ and ‘skewness’ as response variables, and each with the accession specific competitiveness or stress tolerance as explanatory variables. For (4), we used a beta distribution with a quadratic polynomial regression model considering the gardens’ means of ‘intensity’ as response variable and the garden mean of the phenological niche (included as a linear and a quadratic predictor) as explanatory variables. To examine the consistency of the associations across botanical gardens, we added the interaction term for the explanatory variables in the four models (explanatory variable X botanic garden). Subsequently, we used model simplification with the comparing models with and without the interaction terms to find the best model considering an AIC criterion. To ensure that a sufficiently large number of species were included in the analysis for each botanical garden, we only considered seven gardens for which measurements of CSR-strategies were available for more than 30 species (see Table 1). For all models, we visually checked simulated residuals for normality and heteroscedasticity using the DHARMa package (Hartig 2024). All analyses were done in R (Version 4.1.1). The analysis steps performed can then be found in script 2_Flowering_Curves_Analyses.R. To create a graph in the appendix of the article, you can use 3_Analyse_Flowering_Curves_ExampleCurves.R. Description of data files: rawdata_PhenObs_20xx.csv: These files contain the raw data used to create the flowering curves. This data can also be accessed via https://www.idiv.de/research/projects/phenobs/data-access/. CS_values.csv: The data from this file is used as predictors in the models (competition and stress-tolerating score per species). To characterise the species according to their competitiveness and stress tolerance, we used Grime's CSR strategies (Grime 1974). The globally calibrated CSR analysis tool StrateFy (Pierce et al. 2017) combines three functional leaf traits (leaf area [mm2] of one leaf, leaf dry matter content (LDMC) as the dry weight of a leaf divided by the fresh weight [%] and specific-leaf-area (SLA) [mm2*g-1]) to extract percentage proportions of the competitiveness, stress tolerance, and ruderal component of the species. Phenology_stages_2019_2024: The data from this file is used as response variable (FlDu = Flowering duration) and as predictor variable (FlOn = Flowering onset, i.e. phenological niche) Example curves.csv: The data is necessary to create Figure S5. Metadata.xlsx: This file contains the metadata for CS_values.csv, rawdata_PhenObs_20xx.csv, phenology_stages_2019-2023.csv and Example curves.csv in separate data sheets. References: Brooks, M.E., K. Kristensen, K.J. van Benthem, et al. (2017) glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal 9, 378–400. Hartig, F. (2024) _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. R, package version 0.4.7 Pierce, S., D. Negreiros, B.E. Cerabolini, et al. 2017. A global method for calculating plant CSR ecological strategies applied across biomes world‐wide. Functional Ecology 31: 444–457.