Sex-dependent phenological responses to climate vary across species’ ranges

The genus Populus L. (Salicaceae) consists of deciduous and fast-growing forest tree species that are widely distributed in the Northern Hemisphere. Most species from this genus are dioecious and bloom in early spring prior to leaf emergence. Both male and female flowers are borne on pendant catkins (39). They are wind-pollinated and produce large amounts of wind-dispersed pollen and tiny cottony seeds. Reported sex ratios for natural populations of various Populus species range from female-biased to balanced (approximately 1:1) and male-biased (40, 41). Site-specific biases in sex ratio are hypothesized to be associated with environmental conditions, though more studies are needed (39, 41). Individual trees flower for 1 to 2 wk (42).

In this study, we focus on eight wind-pollinated dioecious species from the genus, mainly distributed in North America: Populus alba L., Populus angustifolia James, Populus balsamifera L., Populus deltoides W. Bartram ex Marshall, Populus fremontii S. Watson, Populus grandidentata Michx., Populus tremuloides Michx, and Populus trichocarpa Torr. and A.Gray ex Hook. We selected these species as they i) are wind-pollinated and dioecious; ii) have a large geographic distribution in North America; and iii) have a large number of digitized specimen images available (>500).

We obtained images of digitized herbarium specimens for eight Populus species collected across 119 y (1901 to 2019) in the United States and Canada. Images were downloaded from the Integrated Digitized Biocollections database (iDigBio; (43)) and the Global Biodiversity Information Facility (GBIF; (44, 45)). These records comprised a total of 13,761 specimens collected between 1901 and 2019. In addition, we manually examined and annotated all 120 physical specimens of these species in the Kriebel Herbarium (PUL) at Purdue University. We also obtained phenological determinations of 29 specimens of these species from the Consortium of California Herbaria, which provides phenological information gathered as part of the California Phenology Thematic Collections Network.

Records of herbarium specimens with incomplete or incorrect date information were removed. Dates were determined as incorrect when the month of collection was not within the range of 1 to 12 or the day of collection was not within the range of 1 to 31. We selected records with either geographic coordinates or at least county-level locality information. For each specimen, collection date, year of collection, elevation, latitude, longitude, country, state, and county were extracted from metadata. The presence of flowers and fruits was manually identified for all specimen images at maximum available resolution using ImageAnt (https://gitlab.com/stuckyb/imageant). In addition, floral sex (male/female) was also identified for each specimen image. The small number of specimens that comprised branches from both male and female plants were counted separately for each sex—i.e., we interpreted them as indicating the presence of both flowering males and females at the time and location of collection.

Duplicated data, such as multiple photos of the same specimen, or records for the same species collected on the same date in the same county were removed. Outliers in collection date (day of year) were detected and removed using the interquartile range (IQR) criterion for each species, which is a common approach in detecting and removing outliers in datasets in ecology and other fields (46–48). Using this approach, data points above q0.75 + 1.5 × IQR or below q0.25 − 1.5 × IQR were determined as outliers, where q0.25 and q0.75 represent the first and third quartiles, respectively. The final dataset comprised 2,765 specimens with flowers and determinable sex across eight species from 1901 to 2019 (SI Appendix, Fig. S1).

We also obtained phenological observation data for six of our focal species from the National Phenology Network (NPN), generated by volunteer community scientists in the United States (49). We focused on two phenophases; “leaves” and “open flower” (49). Phenological determinations made by nonexpert volunteers using NPN definitions show high interobserver agreement and are highly congruent with expert assessments (50). To determine the first leaf out and first open flower dates, the first “yes” record that was observed ≤10 d after a “no” record was used (51). To further minimize the sampling bias in observations, only sites with observations for both phenophases for each species were considered. Outliers for each phenophase and species were removed using the same IQR criterion approach applied to the specimen data above. After the data cleaning, we excluded P. angustifolia and P. trichocarpa due to low sample size (<20) and were left with four species (P. balsamifera, P. deltoides, P. grandidentata, and P. tremuloides) for which we had sufficient data (SI Appendix, Fig. S1). The final NPN dataset comprised 1,211 phenological observation collected from 2007 to 2021, including 758 first leaf out dates and 453 first open flower dates across four species.

As precise locality data are not available for the majority of historic specimen records (52, 53), we used county as our geographical unit of analysis. For each county and year, we calculated the mean of each climatic variable and assigned these values to each specimen (4, 9). We used Climatologies at High resolution for the Earth’s Land Surface Areas (CHELSA) monthly climate data with 1-km spatial resolution including maximum and minimum temperature and total precipitation in each month from CHELSAcruts 1901 to 2016 (54) and CHELSA v2.1 2017 to 2019 (55). Mean maximum and minimum temperature and total precipitation in winter and spring were calculated for each county in the United States and Canada. For the NPN observations, which have accurate coordinate data, we used the geographic coordinates and elevation of NPN sites to extract climatic data including mean annual temperature, mean annual precipitation, mean minimum, average, and maximum temperature and total precipitation in winter and spring for each observation site over 121 y (1901 to 2021) available from ClimateNA v. 7.3, which downscales gridded climate data to scale-free point locations (56). Temperature variables were in degrees Celsius (°C) and precipitation variables in millimeters (mm). Winter was defined as a 3-mo period from December to February, and spring was defined as a 3-mo period from March to May, following previous phenological studies from the region (4, 9, 57). A large pool of research suggests that winter and spring temperature are the primary climatic drivers of both the vegetative and reproductive phenology of woody plant species in temperate regions (7, 9, 38, 58, 59).

In order to partition geographic and interannual variation in climate, we calculated the long-term average values (1901 to 2019) of mean maximum and minimum temperature and total precipitation in winter and spring for each county and their temporal anomalies (i.e., the deviation between the annual and the long-term climatic condition) in the collection year for each county (11, 60). As local precipitation conditions and the biological effects of precipitation anomalies of equal magnitude can vary greatly across the geographic range of species, we calculated precipitation anomalies as the proportional value to the long-term mean as suggested by Pearson et al. (11), so the precipitation effects across sites were standardized. The same approach was applied to the climate variables from ClimateNA in calculating the long-term average and climate anomalies associated with NPN observations. Temperature anomalies were estimated in degree Celsius (°C) and precipitation anomalies in unitless proportions. In addition, we used global elevation data (GTOPO30, DOI: https://doi.org/10.5066/F7DF6PQS) with 1-km spatial resolution to calculate county average elevation.

We used linear mixed-effect models to examine the inter- and intraspecies variation in phenological sensitivities to climate change between sex and among ecoregions across species’ ranges. Prior studies and our preliminary analyses suggested that the phenology of P. fremontii populations strongly associate with winter and spring precipitation (61, 62), in contrast to the other seven species in our study whose phenologies are mainly affected by temperature (63, 64). Thus, we conducted separate analyses on P. fremontti.

For the other seven species, we fitted a linear mixed-effect model that estimated the phenological sensitivity to warming between sex, among species, and across the landscape within species’ ranges. Preliminary modeling suggested the mean maximum temperature in spring was the best climate variable affecting flowering phenology of these species by comparing the mean minimum, maximum, and average temperatures in spring and winter. Thus, we included climate anomalies of mean maximum spring temperature and spring precipitation and their interactions in the model representing effects of interannual variations in climate. Additionally, latitude and elevation were included as fixed effect components of the model representing the effects of geographic variation. As we hypothesized that phenological sensitivities to climate may vary among species, between sex, and among populations within species, we included species, sex, and ecoregion as nested random components in the model; specifically, a random intercept term and a random slope term of temperature anomaly for each species, between sex, and among ecoregions. The random intercept allowed the model to incorporate the phenological variation among species in specific sex and ecoregion, while the random slope estimated the phenological sensitivities to temperature anomaly change (i.e., interannual variation) for each group. We used ecoregions classified by the United States Environmental Protection Agency Level II delimitations for North America (65). Only ecoregions with at least 3 observations for both sexes of each species were included in the model analysis to reduce potential bias from small sample sizes. The best model explained 49% of the total variance in Populus flowering phenology (conditional R² = 0.49). A linear mixed-effect model was fitted for P. fremontii separately. Based on preliminary modeling of P. fremontii, we selected the anomalies of winter temperature and precipitation and their interaction in the model to represent the interannual variation of climate. Latitude and elevation were included representing the effects of geographic variation. The random effect of the model included a random intercept term and a random slope term of temperature anomaly between sex and among ecoregions. The conditional R² value of the best model for P. fremontii was 0.47.

To estimate the variation in phenological sensitivity between first leaf unfolding and first open flower dates within each species, we fitted a separate linear mixed-effect model to the NPN observational data. The fixed component of the initial model (i.e., full model) included latitude, elevation, the climate anomalies of mean maximum spring temperature and precipitation in spring, and the interactions between the anomalies. The random component in the model included a random intercept term and a random slope term of temperature anomaly for each species and phenophase. The conditional R² value of the best model of NPN phenological data was 0.58.

Akaike Information Criterion (AIC) values of the resulting models were evaluated to further optimize the fixed component in the model. Models with the lowest AIC value and simplest structure (i.e., smaller number of variables) were selected as the final (best) model. All explanatory variables were centered and scaled to a mean value of zero and a SD of one to avoid bias in model development with variables on different scales. The absolute correlation coefficient values among all predictors in the best models ranged from 0.003 to 0.54, limiting the effect of multicollinearity and overfitting (66). The variance inflation factor values of all variables in the best models ranged from 1 to 1.52, lower than the stringent threshold 3 defined in ref. 67 for avoiding issues of collinearity among variables.

To increase our understanding of how continuous climate change in the future will alter the phenological shifts of the Populus spp. between sex and among regions, we used the best models to predict their flowering phenology across two 30-y time periods (2041 to 2070 and 2071 to 2100) and two shared socioeconomic pathway (SSPs) scenarios (SSP2-4.5 and SSP5-8.5). SSP2-4.5 assumes middle of the road development with medium challenges to mitigation and adaptation, which may be the most realistic development trajectory, while SSP5-8.5 represents fossil-fueled development with high challenges to mitigation and low challenges to adaptation, as the low-effort baseline (68, 69). Future predictions under these two scenarios offer a range of references for climate impact assessment on species interactions. Future climate projections data from 8 CMIP6 GCM ensembles for two 30-y time periods under the two SSPs scenarios were obtained from ClimateNA. To evaluate how phenology will shift in the future, we used the best models fitted above to estimate the current female and male flowering phenology of seven species as a baseline using 30-y historical normals (1981 to 2010) (70, 71). Differences in phenology between current and future predictions were calculated for each species by subtracting the estimated dates from the predicted dates at each ecoregion. The temporal gap between the sexes of each species was calculated under both current climatic conditions and future conditions under two different socioeconomic pathways. Changes in interphase temporal gaps were calculated by subtracting the gaps at current climate normal from the predicted gaps in the future. As the historical normal gap size can vary among ecoregions, we calculated the proportion of changes in the temporal gaps between sex for each species among ecoregions, to facilitate comparisons with future predictions. The same approach was used to predict the future changes of two phenophases (first open flower and first leaf out dates) of four species and the gaps between phenophases for the same two 30-y time periods under two climate change scenarios.

Male flowering phenology was consistently earlier than that of females across all eight species (Fig. 2). With the exception of P. fremontii, the timing of flowering was significantly affected by latitude, elevation, temperature, and precipitation in spring. Specifically, earlier flowering was associated with lower latitude, lower elevation, higher maximum temperature, and lower precipitation in spring (Table 1). Male flowering phenology exhibited stronger sensitivity to changes in temperature than female flowering phenology across all species (Fig. 3A) except P. fremontii as explained below. On average, the phenological sensitivity of male and female flowers to temperature was –2.2 d/°C and –1.6 d/°C, respectively. Phenological sensitivity to changes in spring temperature also varied among ecoregions within the species’ geographic ranges (SI Appendix, Fig. S2). The degree of variation in phenological sensitivity among species was not significantly larger [t (14) = –0.72, P = 0.48, 95% CI (–0.40, 0.20)] than the variation among ecoregions within each species’ ranges (Fig. 3B).

In the case of P. fremontii, the best model included significant effects of elevation and winter precipitation on flowering phenology (SI Appendix, Table S1). Advanced flowering time in P. fremontii was associated with lower elevation and higher degrees of winter precipitation. In contrast to other species, we did not find significant variations in phenological sensitivity between male and female flowers or across the geographic range of P. fremontii.

To evaluate how future climate change may alter the relative flowering times of male and female plants, we made phenological predictions for the seven species with varied phenological sensitivities to warming between male and female flowering. Increased temporal gaps between male and female flowering time were predicted for most ecoregions and species for two 30-y time periods (2041 to 2070 and 2071 to 2100) under two climate change scenarios (SSP2-4.5 and SSP5-8.5; SI Appendix, Figs. S3 and S4). The specific degrees of change in flowering time gaps between males and females varied across species and ecoregions, ranging from −5 to +11 d with the average changes ranging from an increase of 2 to 6 d under the highest warming scenario (SSP5-8.5), compared to the model predicted baseline (1981 to 2010). Greater phenological shifts were predicted under the scenario with greater warming (SSP5-8.5; SI Appendix, Fig. S5). On average, the time between male and female flowering was predicted to increase by 10% across the seven species’ geographic ranges with 90% of predictions ranging from 8 to 28% increases (Fig. 4).

Phenological observations made by community scientists showed that flowering consistently started earlier than leaf out across the four species for which we had data (Fig. 5A). Both the onset of flowering and leaf out were significantly affected by latitude, elevation, maximum spring temperature anomaly, and spring precipitation anomaly (SI Appendix, Table S2). Advanced first open flower and first leaf out dates were associated with lower latitude and elevation, higher maximum spring temperature, and higher spring precipitation. In addition, phenological sensitivity to maximum spring temperature varied between the two phenophases and among four species. First open flower date was more sensitive to change in maximum spring temperature than first leaf out date in all four species (Fig. 5B).

Increases in the temporal gaps between the first open flower date and first leaf out date were predicted for all four species for the future time periods 2041 to 2070 and 2071 to 2100 (SI Appendix, Figs. S6 and S7). The gap between the two phenophases was predicted to increase by 0.4 to 10 d (SI Appendix, Fig. S7) across these species, corresponding to increases of 4 to 52% (Fig. 6). Greater increases were predicted under the scenario with stronger warming. Under the highest emission scenario (SSP5-8.5), the gap was predicted to increase from 1 to 10 d across four species (SI Appendix, Fig. S7). P. balsamifera was estimated to have the greatest increase in the gap between first open flower and first leaf out dates among four species in the future, while the smallest increase was predicted for P. deltoides.

We found similar degrees of inter- and intra-specific variation in the phenological sensitivity of flowering, which is consistent with the findings from prior studies focused on biotically pollinated plants (9). Wind-pollinated species tend to exhibit shorter, more locally synchronous flowering times (24, 25). To overcome the particular problems associated with anemophily, it has been suggested that wind-pollinated plants may benefit from the mass release of pollen on specific environmental cues (21, 75). Interspecific overlap in flowering time may be more common in wind-pollinated species as they need to track environmental conditions facultative to pollen dispersal (e.g., windy, dry) but do not compete for pollination services (76). Further, there is pressure to release pollen when physical obstacles are limited—i.e., before leaf out and canopy closure (21, 77). Therefore, flowering phenology in wind-pollinated plants has been expected to show locally adapted intrapopulation synchronization (21, 24, 25). Our results demonstrate the existence of large variation in phenological sensitivity to climate across species’ ranges, which indicates a high degree of specificity in how the time of flowering responds to climate. It has been suggested that the sensitivity of flowering phenology correlates with long-term climatic conditions, with some studies arguing that sensitivity tends to be greater in colder areas where growing seasons are shorter and thus advantageous to reproduce as soon as possible (78), and others finding greater sensitivity in warmer, more humid regions where climatic cues are more reliable and the risks associated with flowering early are lower (9). We did not, however, find any statistically significant associations between the degree of phenological sensitivity to temperature and the regional climate within the species’ ranges. Further exploration of environmental factors related to wind pollination may help us better understand the distribution of phenological sensitivity across wind-pollinated species’ ranges.

Our results suggest that Populus males generally flower earlier and are more phenological sensitive to changes in climate. Though our specimen data are more likely to represent collections during peak flowering (79), the timing of peak flowering has been demonstrated to correlate strongly with the beginning of flowering and have significant effects on reproductive interactions (80, 81). Also, leaf out preceded flowering in only 19% of the 437 instances of NPN records where flowering and leaf out were observed for the same Populus individual in the same year. Further, male-first flowering in Populus spp. has been frequently documented in the field (42, 82–86). The earlier flowering of males has been associated with several factors, including competition among males for pollination opportunities (87); earlier seed emergence (88); and faster growth rates (89, 90); though these studies were mostly focused on animal-pollinated dioecious species. In the reproduction of dioecious species such as Populus spp., the investment gain of fitness should be mate-limited in males and resource-limited in females (91–93). This is because the reproductive investment of males is mainly restricted to flowering, whereas that of females extends to embryo development and dispersal. Given this disparity in reproductive investment, we may expect males of dioecious plants to be more abundant than females and bloom earlier to maximize their chances of successful reproduction (94–96). Females, on the other hand, may be more averse to risk and adapted to flower once i) climatic conditions for reproduction are safely met and ii) pollen availability is not of concern. Higher phenological sensitivity to temperature in female flowers can be detrimental to reproductive success, as i) warm periods may often be followed by periods of dramatic chilling in temperate regions which can cause floral damage (97), and flowering times can advance faster than the end of the frost season in response to warming (98–101); and ii) females may flower ahead of males when pollen is not available. Alternatively, in the face of rapid climate change, it is possible that females' lesser ability to track climatic variability could result in lowered fitness, as species that do not closely track changes in climate have been suggested to exhibit lower reproductive success, growth, productivity, and abundance compared to those that do (102, 103). We did not, however, find significant responses to temperature or differences in phenological sensitivity among male and female flowers of P. fremontti. Instead, winter precipitation was the climatic factor with the largest effect on the flowering time of this species, which supports findings from previous studies (61, 62). This may be related to the intensive drought and high temperatures prevalent in the species’ range in the southwestern United States resulting in a shorter window of time for reproduction, which may select strongly for more synchronous flowering among sexes during this limited time (61, 104, 105). Interestingly, the specimens we examined comprised significantly more females than males. This may be due to collection bias in herbarium specimens that favors female individuals (14). Along these lines, the targeted collection/observation of male individuals can help reduce this bias and further our understanding of intersexual variation in phenological responses in wind-pollinated plants.

Climate is not predicted to change consistently across species’ ranges, and phenological responses to this change will vary greatly among and within species. For instance, previous studies have suggested that climate change may prolong or shorten the time between leaf-out and flowering (37, 38) and increase the time between the flowering of cooccurring congeneric species (4, 106). Temporal synchrony among sexes is important for effective reproduction, particularly in dioecious, wind-pollinated plant species as they rely upon a nonspecific vector and the probability of successful pollination greatly decreases once leaves emerge (21, 23, 107). Our results suggest that sexual variation in phenological sensitivity will cause increases in the temporal gaps between male and female flowering times of Populus species in response to warming in spring. The predicted divergence in flowering time among sexes are driven by higher degrees of phenological sensitivity to temperature in male flowering. Though temporal overlap in flowering times does not directly translate into successful pollination, reduction in flowering synchrony increases pollen limitation by decreasing pollen availability during the time females are receptive and will likely result in lowered reproductive success, especially given the relatively short pollen longevity and pollination periods of wind-pollinated species such as Populus spp. (22, 85, 108). Further, shifts in the timing of pollen dispersal may also affect other interacting organisms, such as bee species, which have been shown to collect Populus pollen (109, 110), as insects can display different phenological responses to changes in climate from the plants they rely upon (7).

Another factor that influences pollination success is the relative timing of flowering vs. leaf-out. This is especially the case when flowering occurs prior to leaf development, as leaves function as a physical deterrent to pollen dispersal and reception—both those on the same plant and others (21, 111). Our results suggest that the temporal gap between flowering and leaf-out in Populus species will increase with climate change. We expect that this interphase increase itself will not negatively affect the probability of reproductive success, as it will likely translate into longer window of time during which pollen can be dispersed unobstructed. However, the net effect of temporal divergence between sexes and phenophases on pollination efficacy is uncertain. Along these lines, other factors, including flowering duration, pollen viability, intrapopulation phenological variability, and their responses to climate must be considered to ascertain the true reproductive cost of phenological divergence among sexes and phenophases (94). Further, studies of sexual differences in dioecious plants should also consider resource investment across the entire reproductive season, and not just flowering (112).

Finally, as these shifts in phenology are also predicted to vary across species, they may result in altered degrees of inter-specific competition and gene flow (31, 90, 113). Many Populus species have overlapping ranges and can hybridize (114). Given their reliance upon a nonspecific pollination vector, interspecific variation in their reproductive phenology may be critical for maintaining species boundaries. Along these lines, it has been suggested that flowering time displacement could be common in wind-pollinated plants (115–117). Heterogenous phenological shifts across species and their ranges will likely result in altered patterns of phenological divergence/convergence among cooccurring congeneric taxa and hybridization rates.

As our study builds upon herbarium specimens and community science observations that were not originally intended for phenological research, our results may be affected by the gaps and biases that these data inherently comprise (118, 119). For instance, herbarium collections can be geographically biased towards population centers and the sex ratios represented within herbarium collections may not accurately reflect those in the field (14, 118). In general, we had fewer observations of male flowering, and sampling efforts were likely not equal across the entirety of species’ ranges. While we limited our analyses to well-collected species for which large amounts of data (n > 500) were available and removed duplicates and outliers to minimize such issues, our results are only as reliable as the data they are based upon. In addition, it is possible that biases in the collection of male vs female individuals due to differences in abundance, longevity, conspicuousness, or collector preference, may have affected our results (14). For example, the persistence and collection of unreceptive female flowers that have not yet visibly started the transition to fruit can potentially influence estimates of flowering phenology, especially when focusing on the period during which pollination is possible. Further study is needed to investigate the effects of such collecting biases on phenological estimations. Community science observations can harbor biases as well. For instance, observation frequency can be variable and biased toward weekends, and the duration of participation can vary greatly across volunteers (120, 121). NPN provides extensive training materials and specific guidelines for making phenological observations, and local phenology projects often provide additional hands-on training and support. Still, there can be some variation in how different community scientists perceive and document phenological events (122, 123). Thus, with additional data, it is possible that our results may change. Along these lines, ensuring that individuals and organs representing both sexes are collected/surveyed and labelled accordingly would greatly facilitate future research into sexual variation in phenological responses. Nonetheless, previous studies have demonstrated that specimen-derived phenological trends are not likely to be systematically biased and are compatible with those derived from field surveys (79, 124). Further, the consistent patterns we observe across species and their ranges suggest that these patterns are real.

Harnessing the power of natural history collections and community science, we show that plant phenology and its sensitivity to environmental cues can vary across and within species and between sex and phenophases in wind-pollinated dioecious species. This variation will likely manifest as increased temporal gaps between male and female flowering and flowering and leafing, the magnitude of which will vary across species’ ranges. Our study demonstrates the importance of considering multiple dimensions of inter- and intra-specific variation when predicting phenological responses to climate change.

Y.X. and D.S.P. designed research; Y.X., H.T.T., L.G.B., C.H.H., and D.S.P. performed research; Y.X. and D.S.P. analyzed data; H.T.T., L.G.B., and C.H.H. collected data; and Y.X. and D.S.P. wrote the paper.

The authors declare no competing interest.