Data sets and R scripts used for the publication "Species' competitive ability and phenology niche affects the flowering intensity 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 2023 in 14 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 263 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_2023.csv. Data analyses:For the analyses we divided the species into two groups: species that only had one flowering peak, and species that had more than one flowering peak. 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 determine the species' mean values for ‘intensity’, ‘duration’, and ‘skewness’ we first calculated the mean value of each accession across the years within one garden. Using these mean values, we then calculated the mean of every species. 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 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 the three remaining characteristics of the flowering curves (‘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 quadratic 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. To test for possible nonlinear relationships of the phenological niche, we included this variable in the model both as a linear and as a quadratic predictor. For all models, we tested for the variance inflation factors of the explanatory variables which were below two, expressing low risk for collinearity. We further extracted partial R2s for the four explanatory variables by using the partial_r2() function. 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, we used the accession specific means per species of ‘intensity’, ‘duration’, and ‘skewness’ (hereafter garden mean). 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: Competitiveness vs. ‘duration’ (1) and ‘skewness’ (2) and flowering onset vs. ‘intensity’ (3)) For (1) and (2), we used linear models with the garden-means of ‘duration’ and ‘skewness’ as response variables, and each with the accession specific competitiveness as explanatory variables. For (3), we used quadratic polynomial regression models 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 three models (explanatory variable X botanic garden). To run the models, we used the lm() function and Type III ANOVAs for significance testing. Subsequently, we used model simplification with the stepAIC() function to determine whether the interactions were significant or not. 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. For all models, we visually checked the residuals for normality and heteroscedasticity. 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_2023: 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: 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.