Flowering time in plants is a highly variable trait that influences species' resource use and exchange of pollen with con- and heterospecifics. Levin (2009) suggested that habitat shifts within species might cause plastic shifts in flowering phenology, reducing pollen exchange across habitats. Coupled with divergent selection across habitats, diverged flowering time might thus pave the way towards ecological speciation. Some of these ideas may apply across species as well. If close heterospecific relatives share phylogenetically conserved flowering times and negatively affect each other's fitness, habitat shifts to microallopatry might provide a means for local coexistence by close relatives by reducing resource competition, shared enemies, or negative interactions via pollination. Habitat shifts might also select for diverged flowering time, or cause flowering time divergence, if phylogenetically conserved cues arrive at different times across habitats. Here, we ask if flowering phenology is phylogenetically conserved for 208 species at our coastal field site in northern California, whether flowering phenology differs systematically across habitat types, and whether habitat shifts are associated with phenological separation, especially in congeners. Because annuality and perenniality have been shown to be associated with habitat traits and flowering time, we included life history in our analyses as well. We also explore the frequency of habitat shifts between congener and noncongener pairs. We use both field observations and data from Jepson eFlora/Jepson Manual 2 (Baldwin et al. 2012) to explore patterns in flowering phenology. The two data sources were well-correlated across 59 species. Phylogeny, habitat, and life history all influenced flowering time, and habitat and life history were also phylogenetically conserved across 208 spp. Congeners differed in habitat more often than noncongener pairs, and also overlapped more in flowering time. Habitat shifts were not associated with shifts in flowering time in congeners, despite mean peak flowering time differences across habitats, and phylogenetic conservatism in habitat use. Congeners that differed in both habitat use and life history, however, did have the greatest difference in peak flowering dates. Habitat shifts likely play a role in local coexistence of close relatives, but our data do not support habitat-mediated changes in phenology as a possible mechanism. Experimental approaches may elucidate the role of phenology, resource competition, pollinators, and other associates in mediating coexistence of congeners at our coastal California field site.
Understanding the drivers of plant flowering time provides insights into the ecology and evolution of plants (Sargent and Ackerly 2008) and their conservation (Sherry et al. 2007; Davis et al. 2010). Flowering phenology is a highly variable trait that determines the timing of resource use, the shared use of, or competition for, pollinators, as well as reflects plant responses to changing environments (e.g., Morales and Traveset 2008; Davis et al 2010). In 2009, Levin proposed that habitat shifts should pave the way for ecological speciation. New habitats might cause phenological divergence in conspecifics through plastic responses to different conditions, and diverged phenologies across habitats – with divergent selection across habitats – could reduce gene flow across habitats and thus promote speciation (Winterer and Weis 2004; Antonovics 2006; Savolainen et al. 2006).
The same effects of habitat shifts might cause flowering time divergence in heterospecifics that could influence speciation or local coexistence. Levin (2006) suggested that habitat shifts could be associated with increased rates of lineage diversification. If heterospecifics, especially close relatives, share cues that induce flowering, but cues come at different times seasonally across habitats, then habitat shifts may reduce co-flowering in heterospecifics without changes in flowering cue traits (Davies et al. 2013). For example, plants often start flowering in response to drying conditions (Rathcke and Lacey 1985; Warren et al. 2011; Ivey and Carr 2012; Li et al. 2016); thus, plants in drier habitats could be expected to flower earlier than close relatives in wetter habitats. Moreover, annuals often occupy drier habitats than perennials and can differ in flowering time (e.g., Hall and Willis 2006), thus we also consider life history effects on flowering time in addition to habitat effects.
Co-flowering species may interact if they facilitate or compete for pollinator services (e.g., Moeller 2004; Agrawal and Fishbein 2008; Morales and Traveset 2008; Tur et al. 2016), if they hybridize, share floral predators, suffer reduced seed set from heterospecific pollen receipt or loss of pollen to heterospecifics (Arceo-Gomez and Ashman 2014, 2016; Toll and Willis 2018; Christie and Strauss 2020) or if they have peak resource use at similar times (Jensen et al. 2019). Moreover, co-flowering in (congeneric) close relatives in particular has been shown to result in greater rates of heterospecific pollen transfer and reproductive interference (Arceo-Gomez and Ashman 2016; Christie and Strauss 2020). Many of these effects would select for mechanisms that reduce co-flowering.
Here, we ask: what are the relative contributions of habitat, life history, and phylogeny as predictors of flowering time overlap between species, and do habitat shifts increase phenological divergence, as suggested by Levin (2009)? To address the specific challenges of close relatives, we also examine: 1) whether habitat and life-history shifts in congeners result in divergence in flowering time, potentially facilitating local coexistence; and 2) whether habitat shifts are more likely to occur in congeneric close relatives.
Habitat shifts may potentially be a mechanism that could result in changes in flowering time (Mallet et al. 2014). Variation in temperature, photoperiod, and moisture are all important cues for flowering time (Rathcke and Lacey 1985; Eckhart et al. 2004; Marques et al. 2004), and may be phylogenetically conserved (Davies et al. 2013). Co-flowering may be triggered by shared environmental cues like moisture availability or temperature (e.g., Diekmann 1996; Pau et al. 2011), which can differ across habitats (Franks et al. 2007; Jentsch et al. 2009; Levin 2009; Ivey and Carr 2012; Jordan et al. 2015; Anacker and Strauss 2016). Alternatively, co-flowering may be hard-wired, if species share flowering cues like daylength (Marques et al. 2004; Li et al. 2016), which is relatively insensitive to environmental conditions. Thus, for species using photoperiod as the predominant flowering cue, habitat shifts would not be associated with phenological divergence.
Phylogenetic conservatism in cues for flowering may also determine flowering time and flowering time overlap, especially in close relatives (Davies et al. 2013; Li et al. 2016; Lessard-Therrien et al. 2014; Anacker and Strauss 2014). Several studies find phylogenetic signal in plant phenology among large groups of plant species (Davies et al. 2013; Li et al. 2016; Lessard-Therrien et al. 2014), meaning that constraints to phenological plasticity and adaptations may limit phenological divergence among relatives (Willis et al. 2008; Davis et al. 2010). Moreover, beyond phylogenetic signal in flowering time across diverse species, sister and closely related species often have similar flowering periods (Mosseler and Papadopol 1989; Ferguson and Jansen 2002; Debussche et al. 2004; Anacker and Strauss 2014; Li et al. 2016).
Using field-collected phenology data and flora databases of the Jepson Manual 2 and Jepson eFlora (Jepson Flora Project 2021), we explore the association of habitat and phylogeny with flowering phenology overlap at our field site, with a specific focus on congeners.
Materials and Methods
The Bodega Marine Reserve (BMR) is located in Sonoma County (38.3070°N, 123.0660°W) and is part of the University of California Natural Reserve System. BMR covers 362 acres along and adjacent to the California coast north of San Francisco. Because two continental plates meet along the San Andreas Fault through the reserve, BMR is also rich in discrete habitat types, despite being a small reserve. There are sand dunes, coastal grasslands/prairie, wetlands, and rocky coastal bluffs that receive salt spray. Some of these habitats differ in parent material and in moisture availability (Anacker and Strauss 2016).
We used two flowering time datasets: first, the phenology of 59 native forb species in the field (hereafter, ‘field data’) and second, we collected flowering time data for 208 forb species that occur at the field site from the Jepson Manual 2, also available on the Jepson Flora Project (2021). The latter flowering periods are presented by month and are based on a combination of herbarium records and knowledge of the author of the species description (Baldwin, Jepson Herbarium, personal communication). These 208 species consisted of the 59 observed in the field plus 149 more present at the reserve; hereafter, ‘Jepson data’. Using the 59 species observed in the field and for which we also had Jepson data, we asked how well our field observations were correlated with reported phenology from the Jepson data.
Field-based flowering time was censused from January to September 2011, covering the Mediterranean climate flowering season. Observations were collected in biweekly surveys of 59 native species of flowering plants across the reserve, including the most abundant flowering plants that year. We recorded the Julian date of flowering onset and then the number of flowers, fruits, and buds on up to 30 individuals of each species at each census period. We determined the date of peak flowering by selecting the date when the most flowers in the population were open. Peak flowering can be less erratic between years than beginning or end of flowering (e.g., CaraDonna et al. 2014), and reflects the maximum number of flowers open for pollen donation and receipt. We used Julian dates for analyses, as all species at BMR start flowering after January 1 and finish flowering before the end of the long dry Mediterranean summer, well before December 31; thus, circular statistics (Morellato et al. 2010) for flowering dates were not required.
Our data reflect what species are doing at our field site, but 59 species is not a large sample size for the questions we wished to ask, especially with regard to congeners. At the time of this study in 2012, flora-based flowering time data by month were available for 208 of the 295 native forb species at BMR (71.5%) from the Jepson Manual 2; they thus reflect a much larger phylogenetic temporal and spatial scale than our data collections. To estimate peak flowering from the Jepson data, we used the midpoint of recorded flowering time. All flowering time data were converted to Julian dates.
One might expect large discrepancies between our temporally and spatially limited observations and the coarse range-wide data (see Discussion in de Keyzer et al. 2017). However, we found that for the 59 species for which we had both field and Jepson phenology estimates, field peak flowering time was significantly and biologically meaningfully correlated with the Jepson data (r = 0.68, P < 0.001). This result gave us some confidence that estimates from the Jepson could be used to address our questions, and provided a larger sample size encompassing the majority of the native forb community at the BMR (208 of 259 species). Moreover, our results (below) are always in accordance across the field and Jepson data sets, though the field data suffer from a lack of power in some cases.
Habitat Use and Life History
The BMR contains discrete habitat types – coastal grassland and sand dunes – that flank each side of the San Andreas fault that runs through the reserve. In addition, there are seeps, and a freshwater marsh, as well as rocky outcrops and coastal bluffs close to the ocean exposed to extensive salt spray. In a previous study, we found that some of these habitats differ in soil moisture levels (e.g., rocky outcrops are much drier than marsh or bluff habitats; Anacker and Strauss 2016), an important flowering cue. Each of the 208 native forb species was assigned to one of the habitat types based on their occurrence at BMR by botanical experts and reserve managers at that time – Peter Connors and Jackie Sones. Species that occupied more than one habitat were assigned to the habitat in which they were most common; because these habitats are quite discrete and different, all species studied could be assigned primarily to one habitat. At this relatively small reserve (146 ha), habitats were within the foraging distances of many abundant pollinators, notably those recorded for Bombus spp. (e.g., Jha and Kremen 2013).
We built a phylogeny for the 208 native forb species based on molecular sequences for three genes (ITS, matK, and rbcL), downloaded from GenBank (Benson et al. 2012). We supplemented missing sequence data with sequences taken from congeners. The gene by species matrix is available from the Dryad digital repository ( https://doi.org/10.5061/dryad.qfttdz0hqand). Sequences were then aligned using MUSCLE.
For the phylogeny, we first used the software program Phylomatic (Webb and Donoghue 2005) to generate a partially resolved topology that was used subsequently as a topological constraint tree. This phylomatic tree was based on a recent Angiosperm Phylogeny Group tree (R20100428). We then conducted a maximum likelihood analysis in RAxML (Stamatakis 2006), using the phylomatic tree as a topological constraint, a GTRCAT model, and 100 bootstrap replicates. The resulting RAxML tree was fully dichotomous with branch lengths in substitutions per site. Due to the size of the gene matrix, we used the RAxML tree to fix the topology during divergence time estimation in BEAST (Drummond et al. 2012). We constrained several nodes using fossil calibrations from Bell et al. (2010; listed below) with an arbitrary standard deviation of 0.1 Ma. We ran a single MCMC chain for 10 million generations, sampling every 1000 generations. We repeated the analyses twice and combined the resulting posteriors to assure convergence of the posterior distribution. From the combined BEAST posterior, a maximum clade credibility tree was made and uploaded to the open access repository Figshare ( https://doi.org/10.6084/m9.figshare.15135450.v1).
Seedplant, 325; Apiaceae, 33; Asteraceae, 44; Boraginaceae, 59; Brassicaceae, 24; Convolvulaceae, 20; Crassulaceae, 41; Cucurbitaceae, 20; Fabaceae, 56; Gentianales, 71; Iridaceae, 32; Liliaceae, 48; Montiaceae, 83; Nyctaginaceae, 13; Onagraceae, 20; Papaveraceae, 112; Plumbaginaceae, 27; Polemoniaceae, 35; Ranunculaceae, 65; Rubiaceae, 56.
To test for phylogenetic signal in peak flowering time, we calculated Pagel's λ, which provides a more robust measurement than other metrics (Münchmueller 2013), using “fitcontinuous()” function in the R package geiger (Harmon et al. 2008). We also estimated phylogenetic signal in habitat use and life history, using Pagel's λ for discrete traits using “fitDiscrete()” function in the same package; lambda values near zero indicate no phylogenetic signal, whereas values near one indicate strong phylogenetic signal. All analyses above and below were repeated for the field data and the Jepson data sets.
We next tested if phenology was correlated with habitat use and life history, while accounting for phylogenetic non-independence using the phylANOVA() function from the phytools package (Revell 2012). To disentangle the contributions of phylogeny, habitat, and life history on peak flowering time, we also used variance partitioning (Desdevises et al. 2003; Peres-Neto et al. 2006; Gonçalves-Souza et al. 2014). We fit an ecological trait model (habitat + life history), a phylogeny model, and a trait + phylogeny model. To represent the phylogeny in linear terms, we decomposed the phylogeny into a set of principal coordinates (PCs) using the “PVRdecomp” function of the PVR package. We regressed each PC against peak flowering time, retaining those that were significantly related at α = 0.05. For the 59-tip phylogeny, we retained just one PC; for the 208-tip phylogeny, we retained 10 PCs. Next, we derived the four constituent components via subtraction of the R2 values from the three regression models, as described by Desdevises et al. (2003).
To address the patterns specifically in close relatives, we conducted a set of analyses on congener pairs; the mean estimated divergence of congeners for the field and Jepson data was 6.9 and 11.7 My, respectively. First, each congener pair was placed into one of four “shift” categories: none (shared habitat and life history), habitat only, life history only, and both (differ in habitat and life history); the counts were compared using a Chi-squared test. We then compared the peak flowering time difference, calculated as abs[peak flowering species A – peak flowering species B], with the type of shift for every congener pair using a one-way ANOVA. Because of the pairwise nature of the data, we reduced our degrees of freedom in the statistical test to equal the number of unique genera, rather than the number of pairs. For the field data, the number of congener pairs was 17, and the number of unique genera was 10; for the Jepson data, the number of congener pairs was 116, and the number of genera was 42. All analyses were performed in R version 3.0.2 (R Foundation for Statistical Computing, Vienna, Austria).
Effects of Phylogeny, Habitat, and Life History on Phenology
Peak flowering time, habitat affinity, and life history all contained moderate to strong phylogenetic signal (Fig. 1, Table 1). Peak flowering time was also related to both habitat and life history for both datasets, even after accounting for phylogenetic non-independence (Fig. 2).
Phylogenetic Patterns in Flowering Time, Habitat-use, and Life History (Annual or Perennial) Using Pagel's λ. *P < 0.05; **P < 0.01; ***P < 0.001.
Phylogeny also explained substantial variation in peak flowering time using variance partitioning methods (Table 2). For the Jepson data, the total variation in peak flowering time was partitioned into a traits-only component (17%), a phylogeny-only component (17%), a shared trait + phylogeny component (16%), and unexplained variation (49%). The 17% shared trait + phylogeny component reflects the correlation of phylogeny and traits (Fig. 2). The 17% trait-only component suggests that the relationship between the traits and peak flowering times remained after accounting for phylogenetic autocorrelation, consistent with the results of the phylANOVA tests. The 17% phylogeny-only component represents effects of shared evolutionary history. In total, we explained up to 51% of the total variation in peak flowering time using habitat, life history and phylogeny as predictors, and all of these terms were included in the best model to explain peak flowering time.
Plant Flowering Time Partitioned Among Habitat and Phylogenetic Components. Adj. R2 values for partition of variance (bottom of table) were obtained via subtraction of the Adj. R2 values of the top model. Shared and unexplained components are untestable (Peres-Neto et al. 2006)
Flowering Similarity in Congeners
On average, median divergence in peak flowering time was 20 d less between congeners than noncongeners, based on all pairs of flowering species [congeners: 33 mean, 28 median d divergence; noncongeners: mean 53 d, median 49 d divergence; t = 3.3, 16.8 df; P < 0.01 from the field data; similarly, 31 mean and 30 d median for congeners, vs 48 mean and 45 d median from Jepson data; t = 7.1, 117.6 df; P < 0.001). Importantly, 15% of congeners had complete overlap in peak flowering time, while only 9.8% of noncongeners did, based on the larger Jepson data-set.
Despite more similar peak flowering times overall, congeners still exhibited about a 3 wk divergence in peak flowering. The majority of congener pairs differed in one or both habitat or life history traits (70.5% for field data; 51.7% for Jepson data; Table 3).
Difference in Peak Flowering Date (Flowering Time Distance), Habitat Shifts, and Life History Shifts Among Congener Pairs. The number of congener pairs in each shift category (habitat shift only, life history shift only, habitat and life history shift, no shift) differed from the null expectation based on the relative frequency of habitats and life histories; these patterns were stronger for the Jepson data (P < 0.001), but showed the same trend in the field data (P = 0.12).
The number of congener pairs in each shift category (habitat shift only, life history shift only, habitat and life history shift, no shift) differed from the null expectation of the relative frequency of habitats and life history for the Jepson data (P < 0.001), and showed the same trend in the field data (P = 0.12; n = 17 pairs for the field data vs. n = 116 for the Jepson data). Congeners were more likely to partition habitats than noncongeners. Habitat shifts were ∼2 times more common than life history shifts, but alone did not affect peak flowering (Table 3; 29.6 vs 29.3 d for field data and 23.6 vs 23.2 for Jepson data), counter to Levin's (2009) hypothesis.
Plant species that flower at the same time overlap in use of resources and pollinators, and shared floral predators. If flowering time is phylogenetically constrained, either by shared cues or indirectly through shared habitats, then shared flowering time can 1) increase resource competition, 2) facilitate pollination or increase competition for pollinators, and 3) increase opportunities for heterospecific pollen transfer. Levin (2006, 2009) suggested that habitat shifts might cause plants to diverge rapidly in flowering time, assortatively mate, and thereby become reproductively isolated through ecological speciation (e.g., Savolainen et al. 2006; Osborne et al. 2019); he provided a convincing review of divergence in phenology with habitat shifts within species (see also Stam 1983; Winterer and Weis 2004; Gavrilets and Vose 2005, 2007). We explored whether flowering time was phylogenetically conserved, and whether habitat shifts might provide a way for heterospecifics to partition flowering time, potentially reducing negative effects on each other's fitness resulting from reproductive and competitive interactions.
We found evidence for phylogenetic conservatism in habitat use, as many other studies have found (e.g., Cavender-Bares et al. 2006; Kraft and Ackerly 2010), and that species occupying different habitats differed in mean flowering time, after taking phylogeny into account. Thus, both phylogeny and habitat influence flowering time. The predictive power of phylogeny, both in terms of phylogenetic signal in habitat use and life history and on its own, was quite large. Peak flowering time had phylogenetic signal of a magnitude consistent with the findings of previous studies of plant phenology (Staggemeier et al. 2010; Davies et al. 2013; Seger et al. 2013; Lessard-Therrien et al. 2014; Li et al. 2016; Fig. 1; Table 1).
Despite the fact that we found divergence in mean peak flowering time across habitats, consistent with Levin's hypothesis, an analysis focusing only on congeners at BMR found that congeners differing in habitat alone did not diverge more in peak flowering time than those sharing the same habitat type. Congeners were more likely to completely overlap in peak flowering time than were noncongener pairs, and their peak flowering times were less diverged than those of non-congeners. Thus, it is likely that conservatism in flowering times outweighed habitat effects on phenology at the congener level.
Phylogenetic signal in flowering time may arise through several pathways: first, flowering cues may be conserved (Rathcke and Lacey 1985); if based on daylength, they may be largely invariant across a range of habitats (Marques et al. 2004; Li et al. 2016). Initiation of flowering in Mediterranean ecosystems may be especially cued on daylength, as other systems with less predictable rainfall show the importance of daylength in determining flowering time (e.g., Cortes-Flores et al. 2017). We did not find evidence for a second pathway to conserved flowering time in congeners: phylogenetic signal in the types of habitats where species grow, and habitat-specific flowering times. Plants often initiate flowering in response to abiotic cues, such as drydown (Rathcke and Lacey 1985; Warren et al. 2011; Ivey and Carr 2012; Li et al 2016), and we had data showing that habitats differed in water content of soils at BMR (Anacker and Strauss 2016); thus, plants in drier habitats could be expected to flower earlier than close relatives in wetter habitats (Eckhart et al. 2004; Jentsch et al., 2009, Mazer et al. 2021, this issue). However, we did not find that close relatives diverged in peak flowering time when they occupied different habitats. A caveat is that our sampling intervals were coarse, and we would not have detected the 3-d divergence of two Clarkia congeners occupying slightly different habitats (Mazer et al 2021). In our study, divergence in flowering time was greatest in congener pairs differing in both habitat and life history.
Another caveat to our study is that we could not take into account the influence of time since divergence between congeners on their flowering time divergence because phylogenetic relationships/sequence data for congeners at Bodega Bay were not available for many species. When these relationships become known, an analysis taking into account the effects of time since divergence between congeners on the magnitude of phenological shifts would be very informative.
Phenological overlap and the role of phylogeny in constraining peak flowering time may impose challenges to or facilitate coexistence, especially for close relatives (Runquist and Stanton 2013; McEwen and Vamosi 2010; Weber and Strauss 2016). For congeners, challenges to coexistence lie in how resources are partitioned, how much pollen is exchanged, and how enemies are shared between ecologically similar species with similar habitat preferences. Our prior work has shown that congeners compete more intensely at the BMR field site than do less closely related taxa, and also that congeners exhibit spatial overdispersion at the field site (Anacker and Strauss 2014, 2016). Challenges to co-flowering could entail in increased competition for pollinators (Morales and Traveset 2008; Albrecht et al. 2016) or greater heterospecific pollen transfer that reduces fitness (Grossenbacher and Stanton 2014; Runquist and Stanton 2013; Arceo-Gomez and Ashman 2016; Christie and Strauss 2020). Thus, phenological overlap in close relatives may favor habitat shifts to reduce any of these types of negative interactions. Spatial sorting would be the mechanism through which habitat divergence mitigates the negative effects of flowering overlap. We found that congeners were more likely to partition habitats than noncongeners in this study. Only additional field experiments can further elucidate the underlying mechanisms driving habitat partitioning among congeners at BMR.
Our results are also important from a methodological perspective. We had good agreement of field-collected phenological data for 59 species from diverse families taken at a local field site –where character displacement in peak flowering time might be occurring at microallopatric scales– with phenological records for the same species from the Jepson Manual 2/Jepson eflora summarized over the range of the species. This study reinforces the conclusions of a growing body of research demonstrating the value of survey resources like field guides and herbarium records in the study of plant phenology (e.g., Davis et al 2015; Willis et al. 2017; Love et al. 2019).
Another caveat in interpreting our results is that total flowering overlap may not coincide with estimated peak flowering overlap. That said, we chose peak flowering because it was more consistent from year to year than other aspects of phenology in previous long-term phenological studies (CaraDonna et al. 2014, and it represents when most of the flowers in the population are open.
Levin (2009) suggested that habitat shifts might drive ecological speciation through a reduction in gene flow across habitats. He did not find support for this hypothesis across species; habitat shifts were not associated with lineage diversification on islands or mainlands in a literature review (Levin 2006). We find that flowering time is influenced by phylogeny, habitat use and life-history across 208 native forb species occupying our field site. However, we did not find evidence that habitat shifts caused phenological divergence in close relatives, counter to the initial expectations. Habitat shifts were, however, more prevalent in congener pairs than noncongener pairs at BMR, despite overall phylogenetic signal in habitat use. These results suggest that microsympatric congeners interfere with each other. Because congeners share flowering times, flower color and morphology, heterospecific pollen deposition and reproductive interference could be an underappreciated mechanism contributing to selection for habitat divergence in close relatives, in addition to traditionally considered forces of resource competition.
Field data were collected by Anna M. Truszczinski. Support was provided by NSF DEB 11-20387 to S.Y.S., the California Agricultural Experiment Station, and an NSF GRFP to A.M.T. We thank Peter Connors and Jackie Sones for assistance with fieldwork and habitat sources, Jean Burns for discussion, Pedro Peres-Neto for statistical discussions and an anonymous reviewer and Marjorie Weber for comments on the manuscript.