Diversity of Reproductive Phenology Among Subtropical Grasses Is Constrained by Evolution and Climatic Niche

Reproductive phenology is sensitive to climatic changes and is associated with species functional types, distribution ranges, and their corresponding climatic niches. Phylogenetic niche conservatism in reproductive phenology also constrains its diversity and the distribution of species. Therefore, we assessed the effects of photosynthetic pathway, life history, phylogeny, and climatic niche on reproductive phenology. For 190 Poaceae species in subtropical China, we compiled data on flowering onset and reproductive period, functional type (photosynthetic pathway and life history), and 18 climatic variables across the species’ global distributions and used phylogenetic models to determine associations. We found strong phylogenetic signals in flowering onset but not in reproductive period. Photosynthetic pathway and life history have significant interactive effects on both flowering onset and reproductive period, such that C3 annual grasses flowered the earliest and had the longest reproductive period. We found that species with wider climatic niches would flower earlier and have longer reproductive periods. Specifically, species that experience wider ranges of mean annual precipitation and coldest-month temperatures would flower earlier, and species with higher mean annual temperature and wider ranges of wettest-quarter precipitation have a longer reproductive period. This study finds that the diversity of reproductive phenology among subtropical grasses is constrained by evolution and climatic niche and that photosynthetic pathway and life history have an interactive effect on the timing and the duration of reproduction.


INTRODUCTION
Reproductive phenology is key to the success of individual plant reproduction and population persistence (Rathcke and Lacey, 1985;Elzinga et al., 2007;Craine et al., 2012a). Differences in reproductive phenology among co-occurring species minimize direct competition and are often important in maintaining high species diversity in communities (Craine et al., 2012b). However, the reproductive phenology of species is not fixed and may be sensitive to changes in climate (Anderson et al., 2012). Flowering onset and duration, for example, are determined by a network of signals including day length (Colasanti and Coneva, 2009), light quality (Cerdán and Chory, 2003), temperature (Fitter and Fitter, 2002;Colasanti and Coneva, 2009), atmospheric CO 2 concentration (Springer and Ward, 2007), and drought (Franks, 2011). In temperate grasses, vernalization, or a prolonged exposure to cold temperature, is required to induce the expression of the FLOWERING LOCUS T (FT) gene, a key gene in flower development (Corbesier et al., 2007;Tamaki et al., 2007;Lv et al., 2014). Besides environmental factors, recent studies have found that the reproductive phenology of species is also influenced by life history, photosynthetic pathway, and phylogeny (Davies et al., 2013;Du et al., 2015;Cortés-Flores et al., 2017;Munson and Long, 2017). Despite these various findings, there has not yet been a systematic study on how reproductive phenology is shaped by plant functional types, phylogeny, and climatic niches.
As the fifth largest plant family comprising ∼12,000 species, grasses (members of the Poaceae family) are an essential food source for many species, including our own, and thus any potential phenological shifts under climatic change might have huge ecological and economic consequences (Munson and Long, 2017). Phylogenetic niche conservatism (PNC) is the tendency of closely related species to share similar niches and traits (Wiens, 2004;Wiens and Graham, 2005) and has been found in flowering onset in many species from different families (Davies et al., 2013;Du et al., 2015;Li et al., 2016). Although the phenology of grasses has been studied in a small number of crops and dominant grassland species (Craine et al., 2012a;Li et al., 2016), studies that encompass the great diversity of grass species are still lacking. This may be partly due to difficulties of long-term phenological surveys on many species or distinctions between flower and fruit stages for grasses in the field (Primack and Gallinat, 2017). However, since PNC in flowering phenology has been widely found and a number of grass morphological and physiological traits show PNC (Liu et al., 2012;Liu and Osborne, 2015), we also expect to find PNC in the reproductive phenology of grasses.
Divergent structural and physiological traits among C 3 and C 4 grasses Sage et al., 2012;Taylor et al., 2014;Lundgren et al., 2019) may lead to differences in reproductive phenology. In grasses, C 4 species have higher photosynthetic rates, greater water use efficiencies, and faster growth under warm and dry conditions than their C 3 counterparts (Ehleringer, 1978;Way et al., 2014;Atkinson et al., 2016). Therefore, where C 3 and C 4 grasses coexist, C 3 grasses might flower earlier to gain an initial reproductive advantage before intense competition with C 4 grasses for space and other resources begins. In addition, the reproductive duration of C 3 grasses might be longer than that of C 4 grasses because C 4 grasses are adapted to warmer regions where the warm climate delays the flowering onset and shortens the duration (Sherry et al., 2011). Using herbarium records of 16 grass species, a recent study found that C 3 grasses indeed start to flower earlier than C 4 grasses (Munson and Long, 2017). However, whether this pattern in reproductive phenology between C 3 and C 4 grasses is widespread has not been tested.
Life history may be another important trait in determining phenology due to fundamental differences in resource allocation strategies across life history functional types (Lundgren and Des Marais, 2020). While annuals tend to flower once in a single growing season and die, biennials and triennials flower in their 2nd and 3rd years, respectively. In contrast, perennials persist over multiple growing seasons but may have multiple flowering periods per year (iteroparity) or flower once over multiple years (semelparity; Lundgren and Des Marais, 2020). Therefore, to ensure that their sole reproductive effort is successful and to limit competition with established perennials, annuals may flower earlier and longer than perennials. Munson and Long (2017) found that annual grasses flowered earlier than perennial grasses in a sample of 16 species, but only in regions with a higher mean annual temperature. Given that, across subtropical grasses, life history influences the diversity of several functional traits (Liu et al., 2019b), we expected that reproductive phenology would also differ between annual and perennial grasses. Therefore, as multiple environmental and functional factors will affect the reproductive phenology of grasses, more species across multiple ecosystems need to be tested to determine whether photosynthetic pathway and life history have interactive effects on grass reproductive phenology and whether the patterns depend on phylogenetic background.
Climatic niches depict where species occur over space and time and reflect their adaptive strategies (Soberón, 2007;MacLean and Beissinger, 2017). Across species, variation in climatic niche is closely related to variation in functional traits (Thuiller et al., 2004). Species with wide geographic ranges are expected to have higher adaptability and are more able to change their traits to better fit their environment than species with narrow ranges Slatyer et al., 2013). In terms of reproductive phenology, species with larger ranges are predicted to be more able to alter the timing of reproduction in response to changes in the local environments than narrowranging species. For example, flowering onset is likely to be related to climatic conditions, especially temperature (Fitter and Fitter, 2002). Therefore, in climate change scenarios, knowledge of the relationships between reproductive phenology and their climatic niches is urgently needed for predicting the survival and the distribution of species.
Previous phylogenetic studies on flowering phenology have included species across multiple families (Munguía-Rosas et al., 2011;Du et al., 2015), but few have focused on many species within one family, which can assess phylogenetic effects at a more detailed scale (Staggemeier et al., 2010). Furthermore, tropical, and subtropical regions contain over half of the herbaceous species in Poaceae (except Bambusoideae and Pooideae; Watson et al., 1992;Osborne et al., 2014). Grass species in subtropical China are diverse in life history and photosynthetic pathway and have been well surveyed and studied, providing an ideal basis for further study (Liu et al., 2019b). Here we use data on reproductive phenology, photosynthetic pathway, life history, phylogeny, and global climatic niche for 190 subtropical grasses in order to test the following hypotheses while controlling for phylogenetic effects: (1) There are strong phylogenetic signals in two reproductive phenological traits (flowering onset and reproductive period); (2) C 3 grasses flower earlier and have longer reproductive periods than C 4 grasses, and similarly, annuals flower earlier and have longer reproductive periods than perennials. Photosynthetic pathway and life history may have interactive effects on flowering onset and reproductive period; and (3) Species with wider global climatic niches might flower earlier and have longer reproductive periods, with some key climatic niche variables playing more important roles than others, for example, coldest-month temperature may be more important than mean annual temperature.

Species Sampling and Reproductive Phenology
Subtropical China has a typical East Asian monsoon climate. Its central region, Guangdong Province, has a mean annual temperature of 21.2 • C (13.6 • C in January and 28.9 • C in July) and a mean annual precipitation of ∼1,700 mm, of which 80% occurs in the wet season from April to September.
There are 316 herbaceous Poaceae species in the native flora of subtropical China, according to the "Flora of Guangdong" (South China Botanical Garden, 2009). Here we only considered herbaceous grasses and so excluded the 153 species of woody bamboos in this region to keep the same life form. We assigned photosynthetic pathway (C 3 or C 4 species) for each species based on published literature (Grass Phylogeny Working Group, 2012;Osborne et al., 2014), as well as life history (annuals or perennials) and reproductive phenology from the flora for 290 species (South China Botanical Garden, 2009). Of these species, 9.5% are C 3 annuals, 14.9% are C 3 perennials, 26.6% are C 4 -annuals, and 49.8% are C 4 perennials. Furthermore, after matching species with global occurrence and phylogenetic information, we had a final sample of 190 species, which constituted 60.1% of the total native Poaceae species in subtropical China. This sample also well represented the floral pool, with 10.0% C 3 annuals, 12.1% C 3 perennials, 27.9% C 4 annuals, and 50.0% C 4 perennials (Supplementary Table S1).
Following previous studies on phenology using monthly data (Hart et al., 2016), we assigned the initial flowering month (e.g., May is recorded as 5) as "flowering onset" and calculated "reproductive period" as the difference between "flowering onset" and flowering and fruiting end month. For the 11 species whose reproductive period ends in the following year, we calculated the period directly (e.g., for Heteropogon triticeus, which flowers and fruits from October until March, the flowering onset is 10 and the reproductive period is 5). This method not only enabled the correct calculation of reproductive period but also avoided the circular problem of underestimating phylogenetic signals in plant phenology (Staggemeier et al., 2019).

Climatic Data and Climatic Niche Calculation
For climatic data, we obtained geo-referenced locality data for each species from the Global Biodiversity Information Facility 1 , using the R statistical computing environment v. 3.5.1 (R Core Team, 2018), and the occ_download function in the rgbif package (Chamberlain et al., 2014). We carefully checked each species' occurrence data to ensure that the database is representative of the distributions of the grasses (Supplementary Figure S1). Next, we extracted the climatic data for all occurrence points from the WorldClim database (Hijmans et al., 2005), using the function extract in the R package raster (Hijmans et al., 2015). The WorldClim dataset has climatic data at 1-km 2 resolution for 1950-2000. We included nine climatic variables to calculate climatic niches: mean annual temperature (MAT, bio1), temperature seasonality (Ts, bio4, standard deviation across 12 monthly temperature values), maximum temperature of the warmest month (Tmax, bio5), minimum temperature of the coldest month (Tmin, bio6), mean annual precipitation (MAP, bio12), precipitation seasonality (Ps, bio15, coefficient of variation across 12 monthly precipitation values), precipitation of the wettest and the driest quarters (Pmax and Pmin, bio16, and bio17, respectively), and mean annual solar radiation (sr). We did not use precipitation of the wettest and the driest months because the monthly precipitation values have high inter-annual fluctuations, and thus the precipitation values based on quarters are more biologically meaningful for plants. Across the 190 species, the sample size of occurrence data per species ranged from 1 to 170,647 (mean = 2,612).
We estimated climatic niches for each species using climatic variable values extracted from localities across its geographic range (Quintero and Wiens, 2013). The mean and range (maximum -minimum) of the nine climatic variables were calculated for each species, using mean, max, and min functions in R, giving a total of 18 indices (9 × 2) to define the species' climatic niches (Supplementary Table S1). The indices related to temperature, precipitation, and solar radiation can set the species' climatic range limits and potential climatic tolerances.

Phylogenetic Tree
A phylogenetic tree of the 190 study species was constructed by integrating published phylogenies based on chloroplast genes (Supplementary Figure S2). We first extracted 146 species from a Poaceae super tree ∼3,000 species  after checking synonymies based on the GrassBase-Synonymy database 2 . We then searched published phylogenies for different genera and identified 27 more species as congeners in the super tree. Finally, we assembled the last 17 species as polytomies into the 190-species tree using the function bind.tree in the R package ape (Paradis et al., 2004). When we focused on phylogenetic relationships rather than the exact divergence times among species, polytomies and missing branch length information only had negligible impacts on the phylogenetic signals calculated using phylogenetic generalized least square (PGLS; Münkemüller et al., 2012). Therefore, we transformed the branch lengths of the phylogenetic tree using Grafen's method to avoid the effects of setting dichotomies/polytomies (Grafen, 1989).

Data Analysis
To determine if there are phylogenetic signals in reproductive phenology traits, we estimated Pagel's λ. This indicates the degree to which residual trait variation shows similarity based on phylogenetic relationships across species (Pagel, 1999). The λ values vary between 0 and 1 in which, λ = 0 indicates no phylogenetic signal and thus this trait is independent of phylogeny. In contrast, λ = 1 implies that the distribution of trait values across the phylogeny is as expected under the Brownian model of trait evolution (i.e., trait variation completely depends on phylogeny). We estimated λ for both flowering onset and reproductive period using PGLS models, using the function pgls in the R package caper (Orme and Freckleton, 2018).
To establish how reproductive phenology differed with photosynthetic type (PT) and life history (AP, annuals and perennials) and whether interactive effects existed between PT and AP, we also constructed PGLS models. PGLS performs well irrespective of the intensity of the phylogenetic signal, which is ideal for comparisons across large numbers of traits differing in their associations with phylogeny (Revell, 2010). We compared four alternative models (y ∼ PT × AP, y ∼ PT + AP, y ∼ PT, and y ∼ AP), where the response variable was either flowering onset or reproductive period, and selected the best-fitted model that had the lowest Akaike information criterion (AIC) scores (Anderson, 2007).
How reproductive phenology related to the 18 climatic niches (nine means and ranges of MAT, Ts, Tmax, Tmin, MAP, Ps, Pmax, Pmin, and sr) was analyzed using stepwise regression based on PGLS. As many factors in a stepwise regression will give redundant factors in the selected best model, we constructed four potential models and compared them. We first built a full model (M1, with all climatic variable means and ranges as factors), from which we obtained its best-fitted model (M1best) based on stepwise regressions by using function step in R. Next, we built a mean-value-based model (using the nine mean variables as factors) and obtained its corresponding best-fitted model M2. Similarly, a range-value-based model (nine range variables as factors) produced its corresponding best-fitted model M3. Then, we combined significant factors from M1best, M2, and M3 to construct a more comprehensive model (M4). The interactive factor PT × AP was added in M1 to M4 according to the results of the abovementioned analysis. We did not add interactive factors among climatic niches because similar previous studies had already tested and treated them as additive factors (Atkin et al., 2015). We then compared the four models to find the best-fitted model based on AIC scores and reported percentages of variation explained by each factor. Finally, we analyzed the relationships between significant climatic factors and reproductive phenology traits using standardized major axis (SMA) analysis as well as PGLS to determine the role that phylogeny plays on these relationships. In order to test the interactive effects (PT × AP) on the relationships, we ran SMAs for all data and the four PT-AP groups (C 3 annuals, C 3 perennials, C 4 annuals, and C 4perennials) separately. SMA analysis was performed using sma function in the R package smatr (Warton et al., 2012).

Reproductive Phenology and Climatic Niches Among Photosynthetic and Life History Types
The 190 grass species belonged to five C 3 and four C 4 independent lineages across the phylogenetic tree, while annual and perennial species were dispersed randomly across the tree (Supplementary Figure S2). C 3 annual grasses started to flower one month earlier than C 3 perennials, C 4 annuals, and C 4 perennials ( Figure 1A). C 3 annuals also showed the longest reproductive period and C 3 perennials the shortest (Figure 1B). Further PGLS model comparisons showed that photosynthetic pathway (PT) and life history type (AP) had significant interactive effects on both flowering onset and reproductive period. y ∼ PT × AP was selected as the best-fitted model ( Table 1B). The interaction PT × AP explained 70.1 and 71.1% of total variation of flowering onset and reproductive period, respectively, far more than single factors PT and AP explained (Table 2A).
For climatic niches, C 3 annuals, and C 3 perennials occurred in localities with the lowest MAT, while C 4 perennials occurred in the highest, based on multiple comparisons ( Figure 1E). The other three key climatic niche variables (Tmin range, MAP range, and Pmax range) did not differ among the four groups (Figures 1C,D,F).

Relationships Between Reproductive Phenology Traits and Climatic Niches
Phylogenetic generalized least square stepwise model comparison on the relationships between flowering onset and climatic niches selected the range-based model (M3) as the one with the best fit (Table 1C). Further analyses using this model found that, among multiple factors, Tmin range, MAP range, and PT × AP were the three significant factors that explained 19.5, 34.5, and 19.2% of the total variation in flowering onset, respectively (Table 2B). Separate bivariate relationships showed that flowering onset was significantly earlier with both increasing Tmin range and MAP range (Figures 2A,B). However, the relationships were driven by C 4 perennials because the models on the other three groups were not significant (Figures 2A,B and Supplementary Table S3A). For the relationship between reproductive period and climatic niche, the best-fitted model (M4) included both mean and range values of climatic niche variables (Table 1C). Specifically, the significant factors of the best-fitted model were MAT mean, Pmax range, and PT × AP, which explained 20.1, 33.4, and 22.9% of the total variation of reproductive period, respectively (Table 2B). Bivariate relationships exhibited that reproductive period was significantly longer with increasing MAT mean and Pmax range (Figures 2C,D). The relationship based on MAT mean was driven by C 4 perennials, while the relationship of Pmax range was significant for both C 3 and C 4 perennials (Figures 2C,D and Supplementary Table S3B).

DISCUSSION
In this study, strong phylogenetic signals in flowering onset, but not reproductive period of grasses indicate the importance of phylogeny in comparative studies on phenology. Significant interactive effects between photosynthetic pathway and life history on reproductive phenology reflect different adaptive strategies. Climatic niche affects reproductive phenology, with species with wider climatic niches flowering earlier and having a longer reproductive period than those with narrower climatic niches.

Phylogenetic Signals in Reproductive Phenology Traits
The high phylogenetic dependence of flowering onset of grasses (λ = 0.50) was consistent with previous studies that strong phylogenetic signals existed in flowering time across 19,631 species in China (Du et al., 2015) and in flowering onset of 62 meadow species (Li et al., 2016). Generally, it was challenging to detect phylogenetic signals within one family due to stronger PNC across families than within families (Losos, 2008), but our data illustrated clear phylogenetic patterns in reproductive phenology within the Poaceae family. Our results also suggested that during the adaptation and evolution processes of grasses, flowering onset was a key trait under stabilizing selection and was therefore phylogenetically conserved (Donoghue, 2008;Anderson et al., 2012), compared with reproductive period. A potential reason was that flowering induction and onset in response to environmental stimuli, such as vernalization, is an elaborate process that involves many molecular regulations and resource investments (Lv et al., 2014;Xu et al., 2019), which tended to be fixed once evolved as an economical strategy. Reproductive period, on the other hand, can be more labile and vary with environmental changes (de Xavier et al., 2019). Environmental stresses that may occur during the reproductive period means that growth and resource allocation to reproduction needs to be more flexible than just initiating flowering in order to ensure survival.
Furthermore, the phylogenetic signals in flowering onset indicated that non-phylogenetic methods were not appropriate due to statistical independence (Revell, 2010), as all the effects

Residuals 178
The phylogenetic generalized least square models on (A) the interactive effects between photosynthetic pathway and life history and (B) key climatic niches affecting the reproductive phenology traits. The degree of freedom (df), F, and P values of each factor are listed. λ Significant factors in each of the best models are shown in bold. PT, photosynthetic type (C 3 , C 4 ); AP, life history (annual, perennial); MAT, mean annual temperature; Ts, temperature seasonality; Tmin, minimum temperature of coldest month; Tmax, maximum temperature of warmest month; MAP, mean annual precipitation; Ps, precipitation seasonality; Pmin, precipitation of driest quarter; Pmax, precipitation of wettest quarter; sr, solar radiation; mean, mean value of each climatic niche; range, range value of each climatic niche; and ns, not significant. **P < 0.01.
of life history, photosynthetic pathway, and climatic niches on flowering onset were significantly related to phylogeny (Table 1). Therefore, we suggest that phylogenetic relationships should be considered when searching for ecological drivers of reproductive phenology in comparative analyses (Davies et al., 2013;Li et al., 2016).

Reproductive Phenology Traits Among Photosynthetic and Life History Types Reflect Their Adaptive Strategies
Due to the interactive effect between photosynthetic pathway (PT) and life history (AP), only C 3 annual grasses, not C 3 perennials, initiated flowering earlier than C 4 grasses ( Figure 1A). This result agreed with previous findings that C 3 annual grasses started flowering the earliest in western United States (Munson and Long, 2017). C 3 grasses may flower earlier and reproduce longer to compensate for growth disadvantages in competing for space and resources with C 4 grasses (Atkinson et al., 2016). However, in our data, C 3 grasses did flower earlier than C 4 grasses, but only within annual species. Most C 4 grasses, which were late flowering regardless of life history, are in tropical/subtropical areas where the evolution of responses to cold/day length in early spring is not as necessary as for the C 3 annuals inhabiting temperate regions . Our results also agreed with previous findings that annuals flower earlier than congeneric perennials in a long-term observational study of 385 plants (Fitter and Fitter, 2002). Thus, the interactive effect (PT × AP) likely arose from perennials because they had more flexibility in flowering onset and reproductive period during multiple growing seasons (Lundgren and Des Marais, 2020). Our hypothesis that reproductive period would be associated with either PT or AP was not supported due to the interactive FIGURE 2 | Relationships between (A,B) flowering onset and (C,D) reproductive period and key climatic niche variables of the 190 subtropical grasses. Points are colored by photosynthetic pathway and life history (C 3 A, C 3 annuals-blue; C 3 P, C 3 perennials-dark blue; C 4 A, C 4 annuals-orange; and C 4 P, C 4 perennials-red). The models are tested for all data and four groups separately, but only the lines of significant models are plotted. The R 2 , P, and Pagel's λ values for the model of all data are listed, while the coefficients for the other models are presented in Supplementary Table S3. The four key climatic niche variables are significant factors of the best-fitted models in Table 2. effect (PT × AP), with C 3 annuals and C 3 perennials showing the longest and the shortest reproductive period, respectively ( Figure 1B). One reason for the similar reproductive period between C 3 and C 4 grasses might be their relatively narrow local climatic ranges in the study region. The local winter temperature (Dec-Mar) was 5-15 • C, while vernalization requires 0-5 • C for 15-30 days. Therefore, some C 3 grasses originally from temperate regions ( Figure 1E) may have lost their coldness stimulation and initiated flowering just as late as C 4 grasses (Munson and Long, 2017), diminishing the differences between the two. We did not find any literature that directly compared the reproductive period of C 3 and C 4 plants, so it is possible that C 3 and C 4 grasses have similar reproductive lengths. Our observation on C 3 wheat and C 4 corn, which both have a reproductive period of 1 to 2 months, supports this. Another possible reason for a longer reproductive period in annuals than in perennials within C 3 species only is that the vegetative growth and the reproduction of annuals occur simultaneously. This might prolong the reproductive period in annuals, while the two stages are more distinct and flexible in perennials (Albani and Coupland, 2010). A longer reproductive period of annuals was consistent with the findings from 302 C 3 tropical tree species (Bawa et al., 2003).
Early flowering is an important strategy to occupy the cold early spring niches, particularly in temperate species (Munguía-Rosas et al., 2011). With a mean winter temperature of 5-15 • C, the subtropical species studied here have rarely experienced a cold early spring; in contrast, the mean temperature of cold winter was at least below 5 • C in the Western United States (Munson and Long, 2017). Thus, the early flowering of C 3 annuals should be an inherited trait. Indeed six out of nine C 3 annual species that flowered earlier than May were from the subfamily Pooideae. Members of this group are reported to flower earlier than those from other Poaceae subfamilies (Liu et al., 2019a) and are mainly distributed in temperate regions  and acclimated to cold stress (Schubert et al., 2019). Therefore, our results reflected different ecological adaptation strategies of different photosynthetic pathways, life histories, and phylogenetic origins in grass species. Consistent with the patterns that we previously found in morphological and physiological traits (Liu et al., 2019b), life history was more important than photosynthetic pathway in explaining differences in both flowering onset and reproductive period among subtropical grasses.

The Four Climatic Niches Affected Reproductive Phenology After Considering Phylogeny
Our hypothesis that species with wider global climatic niches might flower earlier and have a longer reproductive period, based on geographical range theories Thuiller et al., 2004), was supported by four key climatic niche variables. Specifically, species with wider ranges of Tmin and MAP flowered earlier, and species with higher mean MAT and a wider range of Pmax had longer reproductive periods. These relationships were mainly driven by C 4 perennials (Figure 2). After accounting for phylogenetic relationships, the effect size of models based on climatic niches was small (even for the best multiple-factor models, R 2 ranged from 0.16 to 0.29; Table 1). However, the four significant climatic niches still have important ecological implications.
The minimum temperature of the coldest month (Tmin) was a better predictor of flowering onset than other temperature variables, especially in representing the breadth of temperature niche (i.e., range). Clearly, Tmin could capture the coldness in early spring for flowering induction better than MAT. Previous studies have reported that grass species flowered earlier under lower temperatures (Evans, 1971;Fitter and Fitter, 2002;Sherry et al., 2011); however, we did not find this using our phylogenetic models. A possible reason for the contrasting results is that previous studies have matched phenological records with local temperatures, but our data were sampled across each species' global distribution (climatic niches). This partially explained the small effect sizes of all the PGLS models and revealed that, compared with local temperature as the physiological driver, the range of Tmin was the ecological driver of flowering onset.
Species spanning a wider range of precipitation tended to flower earlier. When species undergo drought, their flowering onset tends to be delayed (Craine et al., 2012a). Species of high-precipitation regions (>2,000 mm, subtropical and tropical areas; Figure 1D) tended to extend their reproductive period by flowering earlier and ending later. Therefore, such a flowering strategy could enhance species fitness and was an evolutionary trend based on historical phenological records (Munguía-Rosas et al., 2011;Anderson et al., 2012).
Mean temperature, but not temperature range, affected the reproductive period of grasses. Mean MAT is associated with the accumulation of temperature over time and is therefore important to plant growth and reproduction (Sherry et al., 2011). For example, compared with temperate grasses with a growth season of 6 months or less, subtropical grasses could grow for a whole year. Mean MAT was also the only factor that differed between C 3 and C 4 species (Munson and Long, 2017; Figure 1E).
Species experiencing a wider range of precipitation tended to have a longer reproductive period. This seemingly contrasts previous findings that reproductive phenology was more sensitive to drought rather than high-level precipitation (Franks, 2011;Craine et al., 2012a), but the drier or the wetter the localities a species occupied, the higher the range of Pmax is. The longer reproductive period is a way to increase the fitness of species with sufficient water resources (Figures 2C,D).
Finally, PNC in flowering onset and its relationships with climatic niche indicated that, when climate changes, a temporal shift in flowering onset was predictable and similar within the same lineage (Staggemeier et al., 2010). Meanwhile, no PNC in the relationships between reproductive period and climatic niches may mean that this aspect of reproductive phenology is less sensitive to climatic change. Although many previous studies had reported PNC in flowering phenology (Lessard-Therrien et al., 2014;Du et al., 2015;Li et al., 2016), few had done phylogenetic analyses on their relationships with climatic variables. For example, Li et al. (2016) tested phylogenetic signals in both phenological and climatic data but did not use phylogeny when examining their relationships. Thus, phylogeny was needed to interpret how reproductive phenology was related with climatic variables.
We acknowledged certain caveats about our analyses. Although we found distinct patterns among groups, more accurate flowering onset and reproductive period data (accurate to date and long-term records) would be favored for future studies. The relatively coarse classification of life history (only the two extreme life histories were considered) and the climatic niches calculated for very widely distributed species might also impede the detection of the key factors in affecting phenology. Sampling species in one family and one region might limit the generality of our results. However, phylogenetic tests within one family could be more reliable than across more distantly related families (Staggemeier et al., 2010). Given that both genetic differences and phenotypic plasticity account for plant responses to climate change (Münzbergová et al., 2017;Malanson et al., 2019), plasticity in reproductive phenology might obscure our results. However, we focused on global climatic niches, which were more stable than local climates (Petitpierre et al., 2012). Furthermore, phenological data from one region could avoid phenotypic plasticity in that measured across larger scales (although it was nearly impossible to measure the phenology of many species in all their global localities), where great plasticity in plant traits occurs (Hoffman et al., 2020). Thus, our trait variation reflected intrinsic phylogenetic differences here , but phenotypic plasticity should also be considered in future studies trying to assess the effects of climate change (De Kort et al., 2020). Overall, we propose that climatic niche and phylogeny constraints should be considered simultaneously in studying the reproductive phenology of subtropical grasses.

CONCLUSION
In this study, we integrated the reproductive phenology, globalscale climatic niches, and phylogeny of 190 subtropical Poaceae species. We found strong phylogenetic signals in flowering onset but not in reproductive period, which emphasized the importance of phylogeny in comparative studies of reproductive phenology. There were significant interactive effects between photosynthetic pathway and life history on both flowering onset and reproductive period, such that C 3 annual grasses flowered the earliest and had the longest reproductive period. We found four key climatic niche factors (ranges of Tmin, MAP, and Pmax and mean of MAT) associated with reproductive phenology through phylogenetic multiple stepwise regressions and demonstrated that species with wider climatic niches would flower earlier and have longer reproductive periods. Our results suggested that future studies on the role of climate on grass reproductive phenology should consider their evolutionary relationships.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
HL and YC designed the study. HL collected the data. HL, KL, JW, LQ, YM, and RZ performed the analyses and drafted the first manuscript. All the authors contributed substantially to the revisions.