Assessing seasonal richness of active flowers throughout UC Reserve sites in the 20th Century

Plant species are well documented to alter both the timing and duration of their flowering in response to changing climate.  Plant species often exhibit different magnitudes or directions of phenological responses to climate changes from each other.  These shifts may have cumulative effects on the diversity of species in flower throughout a given flowering season, resulting in periods of high or low species richness of actively flowering community members that differ from those that occurred under historical conditions.  

In this study we model the effects of warming throughout the past century on the daily species richness of actively flowering species by developing species-specific phenoclimate models for 1,848 plant species documented to inhabit 16 well documented plant communities across California.  These communities encompassed a variety of distinct vegetation types, ranging from coastal marshes and grasslands to chaparral shrublands and mountainous conifer forests.  

By examining consistent patterns in the resultant modeled community-level flowering displays, we demonstrate that recent warming is likely to have consistently shortened the period in which many species flower concurrently, and that the bloom season has advanced by nearly 5 days on average.  Accordingly, within every flora, recent warming was predicted to increase the daily species richness of active flowers early in the local growing season, with corresponding reductions in species richness of active flowers later in the growing season.  Notably, patterns of change in community-level bloom displays were driven primarily by differences among species in the timing of flowering onset, as termination dates tended to advance in unison with onset dates, resulting in minor changes to flowering duration among species. 

Methods
The University of California Reserve sites used in this study consisted of all UC reserves for which plant lists were available, and for which the local floras included at least 100 angiosperm taxa. Plant lists used in this study were posted on-line by each reserve and assembled by Brian Haggerty and Susan J. Mazer (https://ucnrs.org/a-flora-for-the-nrs/z). To ensure that each site represented a community of plants located within a small enough area that they might reasonably be considered sympatric or be accessible to shared pollinators, we further excluded all reserve sites covering areas exceeding 3,000 ha, or in which the range of elevations within the site exceeded 250 meters. The remaining sites included 16 distinct locations distributed across California, including both coastal, inland, and mountainous sites, and represented a range of ecoregions and vegetation classes.

Herbarium Data

Records of flowering phenology used in this study were drawn from 9,216,145 digital specimen records acquired from the digital archives of 440 herbaria throughout North America (Park et al. 2023). In order to ensure the quality of the dataset under examination, we excluded from the data set analyzed here: all specimens not recorded as being in flower at the time of collection; duplicate specimens of a given species (i.e., specimens collected on the same DOY (day of year), in the same year, and at the same location) and specimens for which the date of collection, the latitude and longitude of the collection site, or the species name was not available. As taxonomic nomenclature has varied both over time and regionally across the historical record represented by this dataset, we were concerned that changes in taxonomic nomenclature generated ambiguity in the species’ identifications represented in the data. To rectify this, we standardized the taxonomic nomenclature used to describe each specimen using taxonomic identifications provided by the Taxonomic Name Resolution Service iPlant Collaborative, Version 4.0 (Boyle et al., 2013, Accessed: 30 August 2021; https://tnrs.biendata.org/), which matched outdated or ambiguous identifications to a standardized taxonomic nomenclature.

To ensure that a sufficient number of specimens were observed for each species to model the timing of the bloom display for that species, we also eliminated all species represented by fewer than 100 unique records.

We then identified all species that were reported within at least one of the University of California Reserves based on the species lists mentioned above (https://ucnrs.org/plant-list/). Finally, we eliminated all specimen records of species not present in the species list of at least one reserve. The remaining data encompassed 1,908,706 flowering specimens representing 1,848 species in 234 angiosperm families distributed across North America. Detailed methods and code used to prepare this dataset for the analyses presented here can be accessed via DRYAD at https://doi.org/10.5061/dryad.0gb5mkm9d.

Climate data

Climate conditions in the year and location of each collection were then estimated using the ClimateNA v7.21 software package available at http://tinyurl.com/ClimateNA (Hamann et al. 2013). Climatic parameters evaluated here consisted of mean annual temperatures of Winter (January-March), Spring (April-June), Summer (July-September) and Fall (October-December) months. Similarly, long-term normal conditions (represented by average Winter, Spring, Summer, and Fall temperatures at each location over the years 1911-1940 and 1991-2020) were then estimated at the location from which each specimen was collected, as well as the centroid of each UC Natural Reserve using similar methods. To separate the effects of phenotypic differences in flowering phenology among populations occupying different collection sites that experienced different long-term climate conditions from the effects of plastic responses in flowering phenology associated with year-to-year variation in climate, we calculated for each climate parameter the annual anomaly associated with each digital specimen. That is, for each specimen represented in the digital data, we calculated the difference between the annual conditions in the year and location of specimen collection and the 1911-1940 normal conditions at that location. Positive values for this difference indicate that a given specimen was collected in a warmer-than-average or wetter-than-average year. Sites within the UC Natural Reserve System that were used in this study ranged in mean annual temperature from 4.3°C to 16.7° (Table 1). Note that, for this study, the phenological sensitivity of each species to interannual climate variability was considered to be constant throughout the range of each species.

Modeling plant phenology

We developed species-specific phenological models for each species from the available digital specimen records. Using the available data for each species, we regressed the observed DOYs against local historical climate normals (represented by 1911-1940 January - March, April-June, July-September, or October-December mean temperatures) and anomalies (represented by annual departures from 1911-1940 January - March, April-June, July-September, or October-December mean temperatures) using quantile regression at the 10th percentile (representing population-level flowering onset) and 90th percentile (representing population-level flowering termination) for each species. As Winter, Spring, Summer, and Fall temperatures are often highly collinear, each model was allowed to select temperature normals and anomalies from only the single season that best explained observed phenological variation (based on Akaike information criterion of models created using temperatures from each season).

These methods have previously been demonstrated to accurately predict the dates of flowering onset and termination for species’ population-level bloom displays in response to a given set of climate conditions using simulated natural history collections data (Park et al. 2024). The predicted duration of each species’ bloom display at each location in which it occurred under each climate scenario was then estimated as the difference in DOY between the predicted dates of population flowering onset and termination. To avoid conflating phenological shifts that resulted from plastic responses of individual plants with population-level phenotypic differences in phenology that occur along climate gradients (which can be generated by both plasticity and evolutionary change), both seasonal temperature normals and anomalies were evaluated for inclusion in each model. The partial regression coefficients associated with temperature normals are interpreted as a combination of the phenotypic differences in phenological timing among populations inhabiting locations with differing climatic conditions and the plastic responses to spatial variation in local temperature, while the partial regression coefficients associated with the climate anomalies are interpreted as solely plastic responses to inter-annual variation in local temperature.

Assessing patterns of change in species-level flowering time and duration

Mean changes in species-level flowering times relative to historical (1911-1940 normal) climate conditions across the UC reserve sites examined in this study were assessed by calculating the predicted onset (10th percentile) and termination (90th percentile) of each species’ flowering within each site in which it occurred under both historical (1911-1940) and contemporary (1991-2020) normal climate conditions. Changes between historical and contemporary DOYs of flowering onset, DOYs of flowering termination, as well as in the duration of the local flowering period, were then assessed for each species in each UC Reserve site using pairwise means comparisons.

Predicting historical and contemporary flowering periods within reserves

Using the resulting phenoclimate models, the DOYs of flowering onset and termination were then predicted for every species documented to occur within each site using the climate conditions located at the centroid of that site using ClimateNA v7.21. Typical flowering times during the early 20th century were represented by predictions of flowering onset and termination by each species under 1911-1940 normal temperatures at each site in which it occurred, while typical flowering times during recent decades were represented by flowering times of each species under 1991-2020 normal temperatures. All sites examined in this study experienced warming between these two time periods, with increases in mean annual temperature (MAT) ranging from 0.8°C to 1.4°C.

Using the predicted DOYs of flowering onset for each species within each reserve in which it occurred, we then calculated the number of species predicted to be in flower throughout the (broadly defined) flowering season under historical (1911-1940) and modern (1991-2020) temperatures at each site. To compensate for differences in taxonomic diversity among locations, relative species richness of the taxa expected to be flowering at each site and on each DOY was then calculated as the number of species predicted to be in flower on that DOY divided by the total number of species documented to inhabit that site.

Evaluating changes in peak species richness of the bloom season

Peak species richness of active flowers within each site under historical and modern conditions was measured as mean number of species predicted to be in flower within a UC Reserve site and under a given climate scenario during the 15-day period in which floral species richness (the number of active, co-flowering species) was predicted to be the highest (note that the timing of this period could differ between sites as well as between historical and contemporary climate scenarios). We then evaluated whether differences between modern (1991-2020) and historical (1911-1940) climate conditions resulted in systemic changes to peak daily species richness of flowers across the UC reserve sites using pairwise t-tests. We further evaluated whether the greater temperature changes over the study period were associated with greater shifts in peak daily richness of flowers by regressing the observed shifts in peak daily species richness of flowers within each UC reserve against the observed change in mean annual temperature within that site.

Evaluating changes in timing and duration of the bloom season

The bloom season was defined in this study as the period of time between the first day of the year on which species richness of active flowers exceeded 20% of the local historical maximum, and the last day of the year on which species richness of active flowers exceeded 20% of the local historical maximum. The duration of the bloom season (i.e., that portion of the year in which the majority of flowering occurs) was measured as the number of days during which species richness of blooming taxa remained above 20% of the historical maximum for that site. This duration was calculated separately for each site under both historical and modern climate conditions. We then tested whether systemic changes to the onset date, termination date, and duration of the bloom season had occurred in response to recent climate changes by comparing durations of the flowering season under historical and modern climate conditions using pairwise t-tests. These analyses were conducted to evaluate whether greater site-specific temperature changes over the study period were associated with greater shifts in onset date, termination date, and duration of the bloom season by regressing the observed change of each of these metrics within each UC reserve against the observed change in mean annual temperature at that site. In order to evaluate whether recent climate change had impacted the duration of high species richness of flowers across these sites, we also tested whether the shifts in the duration of the peak bloom season (defined as the number of days in which each site exhibited >75% of its peak historical maximum richness) had occurred in response to recent warming using similar methods.

Climate warming might also cause changes in the patterns of daily species richness of actively flowering species throughout the bloom season that would not be captured through examinations of changes in the duration, timing, or peak richness of active flowers within the bloom season at each site. For example, warming conditions could result in decreased daily species richness of flowers during some portions of the season. Such changes could be evaluated, however, by determining the magnitude of daily changes in species richness of active flowers that occurred within each site as a response to warming conditions. Thus, to evaluate whether larger magnitudes of warming were associated with greater changes in daily species richness of blooming plants throughout the growing season, for each reserve we calculated the mean absolute departures (or mean absolute errors) between daily historical (1911-1940 normal) and contemporary (1991-2020 normal) species richnesses of flowering plants throughout the historical bloom season of the flora at that reserve (defined as the period between the DOY on which that flora’s relative species richness of active flowers reached 20% of its maximum under historical conditions, and the DOY on which it subsequently fell below 20%). Mean absolute departures measure the degree to which warming produced changes to species richness of species in flower within each site. Accordingly, this metric can be thought of as a general estimate of the mean magnitude of change in daily richness of flowering species that occurred throughout the growing season in response to climate warming, irrespective of whether such changes involved increases or decreases in daily species richness of active flowers. In order to determine whether greater warming was associated with greater overall change in patterns of daily species richness of flowers, we then regressed the mean absolute departures in relative daily species richness among sites against the observed differences in mean annual temperature that had been observed between historical (1911-1940) and contemporary (1991-2020) time periods at each site.

Evaluating climate-driven shifts to skewness in seasonal distributions of floral diversity

As individual species collectively shift their flowering phenology in response to changing climate conditions, this may result not only in changes to the timing and diversity of floral resource availability, but also to the shape of the distribution of floral resources across the growing season. Thus, we calculated the skew of the distribution of the species richness of floral resources throughout the growing season under both historical and modern conditions, and then used pairwise t-tests to determine whether phenological responses to modern (1991-2020) conditions exhibited significantly different skew from floral displays under historical (1911-1940) conditions.

Evaluating systemic shifts in species-level phenological timings and durations

Changes in the timing, duration, and peak richness of the community-level flowering display could be caused not only by differences among taxa in the timing of their flowering but by systematic reductions or increases in their flowering durations. For example, earlier onset of the bloom season could be explained by advances in the timing of a small percentage of species, or by systematic advances across the majority of taxa occupying a given site. Similarly, declines in daily species richness could be explained by increased separation in the timing of flowering onset among taxa in response to warming conditions, or by consistent reductions in flowering duration. Thus, in order to determine whether the changes we observed in community-level flowering metrics were due to systematic changes in the timing of flowering or the duration of flowering, we tested whether the plant species examined in this study had collectively exhibited significant shifts between historical (1911-1940) versus modern (1991-2020) conditions in flowering onset, peak flowering, flowering termination, or flowering durations of each species under using pairwise t-tests.