Temperature controls phenology in continuously flowering Protea species of subtropical Africa.

PREMISE OF THE STUDY: Herbarium specimens are increasingly used as records of plant flowering phenology. However, most herbarium-based studies on plant phenology focus on taxa from temperate regions. Here, we explore flowering phenologic responses to climate in the subtropical plant genus Protea (Proteaceae), an iconic group of plants that flower year-round and are endemic to subtropical Africa.
METHODS: We present a novel, circular sliding window approach to investigate phenological patterns developed for species with year-round flowering. We employ our method to evaluate the extent to which site-to-site and year-to-year variation in temperature and precipitation affect flowering dates using a database of 1727 herbarium records of 25 Protea species. We also explore phylogenetic conservatism in flowering phenology.
RESULTS: We show that herbarium data combined with our sliding window approach successfully captured independently reported flowering phenology patterns (r = 0.93). Both warmer sites and warmer years were associated with earlier flowering of 3-5 days/°C, whereas precipitation variation had no significant effect on flowering phenology. Although species vary widely in phenological responsiveness, responses are phylogenetically conserved, with closely related species tending to shift flowering similarly with increasing temperature. DISCUSSION: Our results point to climate-responsive phenology for this important plant genus and indicate that the subtropical, aseasonally flowering genus Protea has temperature-driven flowering responses that are remarkably similar to those of better-studied northern temperate plant species, suggesting a generality across biomes that has not been described elsewhere.

Phenology is the study of periodic life cycle stages, especially how these are influenced by weather and climate (Schwartz, 2013). Phenological change in response to the rapid shifts in climate induced by human activity is a topic of growing interest in ecology as well as conservation and evolutionary biology. In particular, changes in flowering phenology are likely to have downstream effects on important ecosystem processes and biotic interactions, including impacts on the various animal taxa that depend on plants for pollen, nectar, or fruit/seed (Fitter and Fitter, 2002; Parmesan and Yohe, 2003; Visser and Both, 2005; Post and Inouye, 2008).

Whereas some plant species have demonstrated trends toward earlier phenological events with warmer temperatures across space and/or time (Hart et al., 2014), others have shown either a delay in the onset of flowering in response to warming (Fitter and Fitter, 2002) or no significant change in phenological events as the climate has warmed (Keatley et al., 2002). In part, this is thought to be because different plant species use different cues to time phenological events. For example, in many plants, day length (photoperiod) has traditionally been considered more important than low winter temperature (vernalization) in synchronizing flowering to the changing seasons (e.g., Yanovsky and Kay, 2003; Saikkonen et al., 2012), but an interaction between both (and other) cues has been demonstrated in Arabidopsis thaliana (L.) Heynh. (Andrés and Coupland, 2012).

Perhaps as a result of related species using similar cues, phenological shifts are seemingly non‐random across lineages (Willis et al., 2008; Davis et al., 2010; Davies et al., 2013), emphasizing the need to explore phenological change within a phylogenetic framework. However, the phylogenetic conservatism of phenological response has only been tested on a small subset of species (Davies et al., 2013), and has not been explored for entire plant communities with fine‐scale phylogenetic resolution, nor across the broad distributional ranges of numerous co‐occurring species. If phenological responsiveness to climate is phylogenetically patterned within or between lineages, phylogenetic information may have value for understanding phenological cueing mechanisms and forecasting future responses to climate change.

However, data for assessing patterns and processes of phenological change are sparse. Long‐term observational data on flowering, leaf‐out, and fruiting are limited across space, time, and clades, and short‐term warming experiments do not reliably reproduce the effects of long‐term climate change (Wolkovich et al., 2012). A critical bias in long‐term phenology data is that they are available primarily for temperate regions and only in rare cases for the tropics, where most plant diversity occurs. One potential way to overcome the constraints of long‐term field observational data on phenophases is by using historical records in herbaria and museums (Davis et al., 2015; Meineke et al., 2018a, 2019). Although such records have not necessarily been collected expressly for phenological investigations, and therefore present their own biases (Daru et al., 2018; Panchen et al., 2019), a significant body of literature now exists in which historical records have potential for investigating climate‐related phenological trends across plant species (Primack et al., 2004; Bolmgren and Lonnberg, 2005; Coleman and Brawley, 2005; Lavoie and Lachance, 2006; Miller‐Rushing et al., 2006; Bowers, 2007; Houle, 2007; Kauserud et al., 2008; Gallagher et al., 2009; Neil et al., 2010; Park and Mazer, 2018). Like observational studies, however, most herbarium‐based studies on plant phenology are based on collections from temperate parts of the northern and southern hemispheres (see Willis et al., 2017 and references therein).

Here, we investigate phenological drivers for Protea L. (Proteaceae), commonly known as proteas or sugarbushes, an iconic flowering plant genus endemic to sub‐Saharan Africa, with its center of diversity in southern Africa (Rourke, 1982). The genus Protea comprises about 115 currently recognized species, the bulk (ca. 70) of which are endemic or near‐endemic to the Cape Floristic Region, and is a flagship taxon for this world biodiversity hotspot (Vogts and Patterson‐Jones, 1982; Manning and Goldblatt, 2012; Daru et al., 2015, 2019; Gollnow and Gerber, 2015). Widely used in horticulture (mainly as cut flowers), Protea is composed of woody shrublets, shrubs, or trees displaying a suite of floral adaptations (Fig. 1) for a diverse array of pollinators including beetles, birds, and small mammals (Collins and Rebelo, 1987; Wright et al., 1991). These pollinator resources are considered ecologically important—for arthropods, for instance, a phenology‐based study involving four species of Protea showed that seasonal community patterns were significantly influenced by infructescence phenology (Roets et al., 2006).

Figure 1

Representatives of Protea (Proteaceae) showing variation in flower head and growth form. (A) Flower head of Protea burchellii, Caledon, South Africa; (B) P. acaulos, Sir Lowry's Pass, Somerset West; (C) P. aristata flower head; (D) P. caffra subsp. caffra, Fairy Glen Nature Reserve, Pretoria; (E) P. angustata, Kleinmond, Western Cape; (F) green‐leaved form of P. roupelliae subsp. roupelliae tree in flower, Mtentu River Gorge, Eastern Cape; (G) flower head of P. cynaroides; (H) P. scolymocephala; (I) P. nana, Paarlberg Nature Reserve, Paarl; (J) P. magnifica; and (K) P. mucronifolia. Pictures A–K courtesy of SAplants [CC BY‐SA 4.0], Wikimedia Commons.

Climatically, the Cape Floristic Region is temperate to subtropical (Van Wyk and Smith, 2001), with mean annual temperatures varying from 15–16°C at the coast to 17–18°C in interior areas, but lower than 13°C at higher altitudes. Frost is restricted to the inland valleys, and snow falls on the higher mountains. The western part of the Cape Floristic Region receives most of its rainfall in winter (so‐called Mediterranean climate), but to the east the rainfall is more evenly distributed throughout the year. Average annual rainfall is mostly between 300 and 2000 mm, but it is estimated to be as high as 5000 mm on some of the mountain peaks. Elsewhere in southern Africa, the geographic range of Protea falls mainly in areas with summer rainfall (±14 species), whereas about 35 species are dispersed farther north in central and tropical Africa (Rourke, 1982).

Despite its considerable ecological significance and cultural importance, the flowering phenology of Protea remains poorly studied (Pierce, 1984; Johnson, 1993). Even in the case of commercially grown plants, the factors that trigger flowering onset are poorly understood. The great variation that exists in Protea flowering times and apparent flowering prerequisites suggests multi‐factorial control of the Protea flowering cycle (Hoffman, 2006). In particular, for some cultivated members of Protea (e.g., Protea cynaroides (L.) L. [Fig. 1G]), flowering may not depend strictly on photoperiod, and a cold treatment during winter may be required by some species to trigger flowering (Hoffman, 2006, and references therein), whereas in the case of resprouters, the age of vegetative growth since the last fire may also play a role (Rebelo, 2008).

The genus as a whole can be described as having a year‐round flowering phenology. Each month of the year has multiple Protea species that are described as being in full flower, with spring and summer being most common (Rebelo, 2001), and while individual species tend to exhibit flowering peaks in certain seasons, most species are known to flower year‐round. Unlike highly seasonal temperate floras with discrete flowering windows, the relatively aseasonal cycle and year‐round flowering for most angiosperms in subtropical systems poses analytic challenges for statistical models of flowering date. Whereas prior studies have addressed the circularity of phenological data sets using circular statistics that convert dates to angles and analyze them in a trigonometric framework (Morellato et al., 2010), formal circular statistics have not, to our knowledge, been integrated with modern inferential techniques, such as hierarchical mixed effects modeling, that are now widely used in ecology. In the context of herbarium data, the use of hierarchical models is advantageous as they can allow data that are sparsely sampled across space and time to be pooled across closely related species to gain inferential power. To take advantage of linear mixed effects modeling while working with these year‐round data, we here develop a new approach that uses sliding windows to linearize the calendar year separately for each species based on its period of minimum flowering.

Here, we use a database of 1727 carefully vetted herbarium specimens representing 25 species collected between 1950 and 2011 to explore flowering phenology across time and space for Protea. Specifically, we: (1) characterize seasonal and geographic flowering phenology patterns across Protea species, (2) investigate how site‐to‐site and year‐to‐year variation in temperature and precipitation influence Protea flowering phenology, and (3) test for phylogenetic conservatism in these climatic effects on phenology. Our study reveals how an iconic plant genus in an understudied part of the world responds to climate variation, and we discuss how this might cascade to affect its native ecosystems.

MATERIALS AND METHODS

Herbarium data

We compiled information from 7770 herbarium specimens for 87 species of Protea. The data were collected from herbaria across South Africa and are archived at the National Herbarium of the South African National Biodiversity Institute in Pretoria, and include records from Eswatini and Lesotho as well as South Africa. From each of these specimens, we extracted four types of data: species identity, collection date (year, month, day), geographic coordinates (often represented as quarter degree grid cells [QDGC]), and flowering status (whether or not the specimen was in full flower). Each QDGC was converted to decimal degrees following Larsen et al. (2009). For example, a QDGC of 3419AD = longitude 19.375 and latitude −34.375, or 3318CD = longitude 18.375 and latitude −33.875, and so on (Larsen et al., 2009). For each specimen, full flowering was assessed using the first‐ and second‐order phenological scoring protocol of Yost et al. (2018) by first examining whether any reproductive structure is present, before determining whether flowers were in anthesis (see also Primack et al., 2004; Hart et al., 2014). Specimens not in full flower were removed from the analysis. We also removed records falling farther north than the northernmost point in South Africa, records from before 1950, and records listed as collected on the last day of any month (these dates were dramatically overrepresented in the data, suggesting that they were often used as a default when the true collection day was unknown). Lastly, species with fewer than 50 remaining records were removed. After data cleaning, the final data set included 1727 records representing 25 species of Protea.

Phenological patterns

For each species, we used the frequency distribution of specimen collection dates as a proxy for flowering phenology (Panchen et al., 2012). We converted dates on specimen labels to Julian Day of Year (DOY; where January 1 = 1 DOY, February 1 = 32 DOY, and so on). To characterize the flowering phenology for each species, we defined the “peak season” as the center of the sliding three‐month window with the largest number of occurrences, “low season” as the center of the sliding six‐month window with the fewest number of occurrences, and “aseasonality” as the ratio of the numbers of low season to peak season occurrences. These sliding windows were wrapped to reflect the circular nature of the calendar year, in which December 31 is adjacent to January 1. For each species, we calculated an adjusted version of the Julian DOY of collection for each herbarium specimen, day of flowering year (DOFY), measured as the number of days since the middle of a species' low season. Unlike DOY, which has a disjunction between December 31 and January 1, DOFY can be treated as a linear variable.

To test the assumption that the seasonal distribution of herbarium specimens with 50% open flowers is a good representation of actual peak flowering in the field, we compared the center of each species' peak flowering season derived from our data set to phenograms provided in the Protea Atlas of Rebelo (2001), which represents the most comprehensive treatment of all described Protea species, including relevant information on the ecology, spatial distributions, and species abundance.

Climate analyses

Studies have shown that the distribution of plant species richness in southern Africa is driven largely by shifts in rainfall and temperature regimes (O'Brien, 1993; O'Brien et al., 1998, 2000). We used temperature and precipitation data from the University of Delaware Air Temperature and Precipitation dataset version 3.01 (Willmott and Matsuura, 2001), which includes 0.5‐degree gridded monthly mean temperature and precipitation for every month of every year from 1950–2010. We extracted the full climate time series (2 variables × 12 months × 61 years) at the collection location of every herbarium specimen. Precipitation values were log‐transformed for normality and for ecological realism.

We used these data to explore the roles of spatial and temporal climate variation in driving Protea flowering dates. For each specimen record, a spatial and a temporal climate anomaly were calculated for temperature and precipitation for four predictor variables. Spatial anomalies were calculated as the difference between a location's long‐term mean climate (across all months and years in the climate data set) and the species‐wide average long‐term mean climate across all specimen locations, with positive values representing specimens from warmer or wetter parts of a species range and negative values representing specimens from cooler or drier parts of the range. Temporal anomalies were calculated as the difference between annual temperature or precipitation at the collection location in the year the specimen was collected and the long‐term mean climate at that location, with positive values representing collections in years that were locally warmer or wetter than average and negative values in years that were locally colder or drier than average. Because species flower at different times of year, these annual climate anomalies were defined differently for each species, calculated as the average across the 12‐month period from eight months before through three months after a species' peak flowering month; this asymmetrical window was chosen in order to encompass climate during a given peak flowering season and the preceding low season, which together are likely to influence flowering phenology. DOFY anomaly, the dependent variable, was calculated as the difference between each specimen's DOFY as described above and the average DOFY for that species.

We fit a single mixed effects model using the R package lme4 (Bates et al., 2015), predicting DOFY anomaly as a function of these four climate variables (spatial temperature anomaly, spatial precipitation anomaly, temporal temperature anomaly, and temporal precipitation anomaly), with random effects of species on slopes but not intercepts (because intercepts are by definition zero as a result of the de‐meaning described above). This hierarchical modeling approach allows the simultaneous estimation of each climatic predictor on Protea flowering phenology both overall and for each individual species, both levels of which are of interest in this study. The maximum likelihood optimization criterion was used over restricted maximum likelihood, to allow significance testing via model comparison. To test for the statistical significance of each of the predictors, the full model was compared to four reduced models using likelihood ratio tests, each with one of the four predictors removed.

Phylogenetic analysis

Phylogenetic relationships among Protea species were reconstructed using DNA sequences from four plastid (trnL, trnL‐trnF, rps16, and atpB‐rbcL) and two nuclear regions (ITS and ncpGS) available from GenBank. The sequences were aligned using SeaView version 4.5.4 (Gouy et al., 2010) and manually adjusted using Mesquite version 2.5 (Maddison and Maddison, 2008). The combined sequence data set comprised 3386 loci.

Next, we reconstructed phylogenetic relationships using maximum likelihood (Stamatakis et al., 2008) via the CIPRES gateway (Miller et al., 2009). Branch lengths were transformed to millions of years by enforcing topological constraints assuming the APG III backbone from Phylomatic version 3 (Webb and Donoghue, 2005). We then used BEAUti version 1.7.5 (Drummond and Rambaut, 2007) to generate the dated phylogenetic tree using Bayesian inference and one independent fossil calibration with normal prior distribution as follows: Protea root node (28.4 Ma, SD 2 Ma). We carried out a Bayesian Markov chain Monte Carlo analysis by running four chains simultaneously for 2 million generations and discarding the first 20% of trees as burn‐in. The distribution of posterior probabilities from the different chains was assessed by constructing a 50% majority rule consensus tree for further analysis.

We then used this dated phylogenetic tree to estimate phylogenetic signal on five aspects of Protea flowering phenology: aseasonality index, spatial temperature anomaly, spatial precipitation anomaly, temporal temperature anomaly, and temporal precipitation anomaly. We tested whether closely related species tend to exhibit similar phenologies or diverge in the timing of reproductive events more or less than expected by chance based on these facets of Protea flowering phenology. We used Abouheif's Cmean statistic (Abouheif, 1999), Blomberg's K (Blomberg et al., 2003), and Pagel's lambda (λ) (Pagel, 1999). Significance was assessed by shuffling the trait values 1000 times across the tips of the phylogeny and comparing it to expectations by random models. Values of Blomberg's K, Abouheif's Cmean, and Pagel's λ > 1 indicate high phylogenetic signal, i.e., closely related species share more similar traits than expected by chance. Both Blomberg's K and Pagel's λ were calculated using the R package phytools (Revell, 2012), whereas Abouheif's Cmean was calculated using adephylo (Jombart and Dray, 2008).

RESULTS

Peak flowering over time and across geographic space

We found wide variation in flowering phenology across geographic space, climate, and time. Specimen records were not evenly distributed across years; earlier years showed sparser records, and a high density of collecting was seen between the 1960s and 1980s (Fig. 2). Our cleaned and validated data set confirmed the year‐round flowering phenology of Protea: across the 25 species, 11 of 12 months contained the center of the peak flowering window of at least one species, more than half of species had flowering observations in all 12 months of the year, and no species had flowering records in fewer than eight months. Each species exhibited a period of peak flowering, with the timing of these peaks varying widely across the year among species (Appendix S1). Importantly, we found a strong correlation between peak flowering date derived from herbarium specimens and flowering time as recorded in the literature (Rebelo, 2001; r = 0.93, Appendix S2), supporting our claim that collection dates on herbarium specimens can serve as surrogates for flowering dates in Protea, given careful data cleaning and validation of specimen flowering status.

Figure 2

Temporal distribution of collection date of herbarium specimens of Protea species in South Africa, spanning 1950 to 2011 (n = 1727 specimens representing 25 species).

Geographically, we found a spatial gradient in peak flowering season of Protea species. In the Mediterranean‐type Cape Floristic Region with wet winters and dry summers, flowering time for most species tended to peak in the winter, whereas the non‐Mediterranean regions with wet summers and dry winters extending from Mthatha in Eastern Cape Province to Limpopo in the north showed more flowering during summer (Fig. 3A, C). Orthogonally to this east–west gradient, a coast–inland gradient in aseasonality of flowering was also apparent, with inland species tending to exhibit larger differences in flowering between the peak and low seasons (Fig. 3B).

Figure 3

Geographic and climatic patterns of flowering times of Protea species in South Africa derived from preserved herbarium specimens. (A) Day of specimen collection represented in Julian day of year (DOY; where Jan 1 = 1 DOY and Feb 1 = 32 DOY, and so on). (B) Species aseasonality of flowering, i.e., the ratio of low season to peak season specimen counts, with low values representing species with most flowering concentrated during peak season and high values representing species with a significant portion of total flowering occurring in low season. In panels A and B, a small amount of noise was added to geographic coordinates to prevent overplotting of points at the centroids of quarter degree grid cells. (C) Monthly mean temperature and precipitation and peak and low flowering months for the 25 Protea species included in the analysis; species are ordered by the slope of the monthly relationship between temperature and precipitation, roughly translating to the gradient from wet‐winter Mediterranean‐type climates in the Cape region to wet‐summer climates in the northeast. N = 1727 specimens for 25 species.

Effects of climate variation on flowering phenology

For temperature, model selection based on likelihood ratio tests identified highly to marginally significant effects of both spatial variation (χ = 14.45, df = 5, P = 0.013) and inter‐annual variation (χ = 10.67, df = 5, P = 0.058) on specimen collection dates. Neither spatial nor temporal variability in precipitation had significant effects (P = 0.93 and P = 0.75, respectively). The fixed effect coefficients were −3.83 days/°C for spatial temperature variation and −5.18 days/°C for year‐to‐year climate variation (Appendix S3), indicating that both dimensions of climatic variability similarly cue the timing of peak flowering season in Protea (Fig. 4).

Figure 4

Fitted model predictions of phenological responsiveness of Protea species across South Africa as a function of spatial and temporal variation in temperature. Larger x‐axis values represent warmer parts of a species' geographic range, while larger y‐axis values represent warmer years at a given location. Each point represents a Protea specimen in the data set, colored according to its predicted flowering date relative to the species mean (i.e., DOFY anomaly). Variability among nearby points is caused by model variables not pictured in the chart, including precipitation variation as well as species‐specific responses to temperature. Diagonal contour lines represent the overall fixed effects of the two temperature variables on flowering phenology, and are spaced at 1‐day intervals.

In addition to the overall coefficients describing climate effects on phenology, the hierarchical model also includes individual species‐level coefficients. Eighty‐eight percent of species exhibited advancement in flowering phenology in warmer locations within their ranges whereas 56% exhibited advancement in warmer years, with well‐known species such as P. cynaroides and P. scolopendriifolia Rourke both showing greatest advancements of −9 days/°C each (Appendices [Link], [Link]). These sensitivities were not significantly predicted by the aseasonality of species' flowering phenologies (Appendix S6).

Phylogenetic signal in flowering phenology

We tested the hypothesis that closely related species shift flowering phenology more similarly (either toward early or late flowering) than expected by chance. Using various dimensions of Protea flowering phenology in relation to climate, we found significant, but weak phylogenetic signal in aseasonality, temporal temperature, and temporal precipitation anomalies (Abouheif's Cmean = −0.017, 0.22, and 0.15, respectively, all P < 0.01; but both λ and K = [ns]; Appendix S7), showing that species within lineages shift flowering time more similarly with these independent facets of climate variation than expected by chance (Fig. 5).

Figure 5

Phylogenetic conservatism of flowering phenology in relation to temporal temperature anomaly in Protea, i.e., the tendency of closely related species to change flowering time similarly under given temporal changes in temperature. The color scales correspond to mixed effects model coefficient for effect of temporal temperature anomaly on flowering date with warming temperature, with shifts toward early flowering indicated in red and late flowering indicated in blue.

DISCUSSION

In this study, we explored shifts in flowering phenology using herbarium specimens of Protea species in southern Africa, providing the first assessment of phenological responses to climate in Africa, within an area unrepresented by historic observational data. We show that across temporal and spatial gradients in climate, Protea species display advanced flowering phenologies in response to warmer temperatures. Although species vary in their phenological responses to climate, these responses are phylogenetically conserved, such that closely related species tended to shift flowering time similarly with temperature.

We found that flowering phenology in Protea species advanced by an average of 3–5 days per degree of temperature across both space and time. Our analysis combined these independent dimensions of climate variation—across sites and across years—and found very similar responses of Protea species to temperature variation along both dimensions. This strengthens our faith in the results, as simple first principles would indeed predict that these two uncorrelated axes of temperature variation would generate similar flowering phenology responses. It also implies that space‐for‐time substitution, a widely practiced but less validated approach for predicting ecological impacts of future climate change such as species geographic range shifts, may have viability in plant phenology research (Pickett, 1989). Our results imply that Protea phenology may be sensitive to ongoing anthropogenic climate change, although our data were not sufficient to address this question directly due to low coverage post‐1980 and high variability in flowering dates in the data set.

Our results are also in agreement with studies from the northern hemisphere showing that plants accelerate their reproductive phenologies with warmer spring temperatures (e.g., Primack et al., 2004; Roberts et al., 2015). These results agree qualitatively but, importantly, are also remarkably similar in effect size. For example, Primack et al. (2004) and Miller‐Rushing and Primack (2008) found slopes of roughly −5 days/°C and −3 days/°C, respectively, for flowering phenology of species in highly seasonal temperate environments of North America. Whereas most such studies have focused on temperate latitudes with seasonal flowering times (Primack et al., 2004; Amano et al., 2010; Panchen et al., 2012; CaraDonna et al., 2014), our study is the first to demonstrate similar variation in phenology for a genus with year‐round flowering. The congruence between our results in subtropical Africa with aseasonal flowering to those from temperate regions with highly seasonal phenologies suggests a general trend across disparate ecosystems and species with different natural histories. However, the herbarium records used here lack precise coordinate data and thus are assigned to centroids of quarter degree cells. Therefore, it is possible that elevational differences might lead to spurious inferences in some cases. Despite the coarseness of the data described here, the congruence of flowering responses to temperature in ours and previous studies lends credence to the use of flowering data from herbarium specimens despite spatial uncertainty inherent in most specimen‐based data sets.

The data preprocessing methods employed here allow for linear regression methods previously used in phenological studies in temperate regions (e.g., Primack et al., 2004; Miller‐Rushing and Primack, 2008) to be applied to more aseasonal systems by using sliding windows to re‐center each species' observations on periods of maximum and minimum flowering activity. When comparing our assessment of flowering phenology to expert‐derived estimates in the literature (Rebelo, 2001), we found very high correlations (r = 0.93), indicating that herbarium records combined with a sliding window approach can indeed capture key phenological patterns. Thus, the methods we detail here may be of broader use across tropical and subtropical regions, which contain most of the world's plant diversity but remain largely unexplored with regard to climate change and phenology.

Phylogenies can provide important insights into species phenological responses to climate change (Davis et al., 2010; Davies et al., 2013). We found that the variability of flowering phenology across species ranges are phylogenetically patterned, such that close relatives have shared phenological responses. Although few studies have investigated phylogenetic conservatism in plant phenology, those that have generally find that phenological responsiveness is phylogenetically conserved (Davis et al., 2010; Jia et al., 2011; Lessard‐Therrien et al., 2013; Davies et al., 2013). For example, Willis et al. (2008) found that phenological responses to climate among plant clades in New England are phylogenetically conserved such that species within less responsive lineages correspond to those in severe decline. We found weak phylogenetic signal, which we believe is a result of limited taxon sampling—25 species (out of about 115 currently recognized species in the genus)—and increasing the taxon sampling might yield stronger phylogenetic pattern. Nonetheless, our results can provide important baselines for more focused investigations of, for example, mechanisms underlying phenological response to climate, and the formation of reproductive barriers that lead to reproductive isolation and sympatric speciation.

In this study, we show that warmer temperatures cue earlier flowering in Protea, but how this phenomenon influences pollinator abundance remains poorly understood. If phenological responsiveness to climate occurs independently among plants and pollinators, some aspects of plant–animal associations including that of pollinators may be modified, potentially leading to phenological mismatch (Kudo and Ida, 2013). We found weak but significant phylogenetic signal in the affinity of closely related species to shift phenology similarly, which can potentially influence co‐evolved mutualists including pollinators. Although we did not explicitly test for phenological mismatches, we believe that herbarium specimens can serve as a critical first step in monitoring species interactions (e.g., Meineke et al., 2018b) and resulting population dynamics, including distributional migrations, expansions, and contractions (Feeley, 2012). In the future, it may even be possible to forecast flowering (Park et al., 2019) and bee pollination to address the conservation challenges posed by changing phenology.

APPENDIX S1. Specimen collection frequency across day of flowering year (DOFY), a normalized version of the Julian day of year. Red vertical dashed lines correspond to January 1.

Click here for additional data file.

APPENDIX S2. Comparison of species peak flowering season (in Julian days) recorded from herbarium specimen records versus the literature (Rebelo, 2001). Rebelo (2001) reports both a “long” season of increased flowering activity and a narrower “short” season of maximal flowering activity for each species, the centers of which are shown here relative to the peak flowering date we calculated from herbarium data as described in the text. Although the y‐axis ranges from 0–365, the x‐axis has a slightly broader range—given the circular nature of the calendar year, a given Julian date can take multiple values (e.g., 10 = 375), and the value that best communicates alignment with the field guide data set is shown.

Click here for additional data file.

APPENDIX S3. Parameters used to characterize phenologic responsiveness to climate in Protea species, estimated from the mixed effects model.

Click here for additional data file.

APPENDIX S4. Changes in flowering times of Protea species across South Africa in relation to anomalies in temperature. Statistical analysis based on mixed effects model using both spatial temperature variation (A) and temporal climate (year‐to‐year temperature variation) (B) as predictors, with species as random effect. Negative slopes indicate advancement of flowering with warming. Lines indicate fitted slopes for individual Protea species. Points indicate input specimen data, and have been truncated for visualization at the extremes of the y‐axis range.

Click here for additional data file.

APPENDIX S5. Species‐specific statistics generated by the sliding window phenology analysis and the mixed effects model (MEM) climate analysis for each of the 25 Protea species.

Click here for additional data file.

APPENDIX S6. Relationship between the aseasonality of species' annual flowering phenology cycles (aseasonality index) and their estimated phenological responses to temperature variation across space and time (coefficients from the linear mixed effects model). Dashed lines show linear regressions with 95% confidence intervals shaded.

Click here for additional data file.

APPENDIX S7. Tests of phylogenetic signal in different dimensions of Protea flowering.

Click here for additional data file.