The study was conducted in Xilinhot, Inner Mongolia, China. The climate is a semiarid temperate continental climate, and nearly 70% of precipitation occurs from May to September. The vegetation type is typical steppe, and the soil is Kastanozem according to the WRB (Dong et al., 2021). The growing season starts at the end of April and ends in late September.

In this study, we analyzed data collected from two different grassland habitats in Xilinhot ( Figure 1 ). The first site is located at the Inner Mongolia Grassland Ecosystem Research Station (IMGERS, 43°33′N, 116°40′E), approximately 70 km southeast of Xilinhot (site ERS hereafter). The altitude is approximately 1,190 m above sea level. The mean annual temperature and precipitation are 1.1°C and 331.4 mm, respectively (1982–2020). In 1979, a fence (24 ha, 600 m × 400 m) was constructed, and since then, all grazing has been excluded. Before fence construction, the grassland was severely degraded due to heavy grazing. After decades of restoration, the vegetation has been significantly restored. The dominant species are Leymus chinensis and Stipa grandis.

The second site is located at the Xilinhot National Climate Observatory (44°08′N, 116°20′), approximately 15 km east of Xilinhot (site NCO hereafter). The altitude is 1,129 m above sea level. The mean annual temperature and precipitation are 3.2°C and 273.7 mm, respectively (1982–2020). The climate at NCO is drier and the species richness is lower than that at ERS. NCO has free grazing with low intensity before the experiment began, and the dominant species is S. grandis. In 2012, three 1.44-ha (120 m × 120 m) fences were constructed to prevent grazing (more details are provided by Dong et al., 2021).

Sampling at ERS was conducted once every half month in the growing season. Sampling began on May 15 and ended on September 15, so nine samplings were performed each year. For each sampling, 20 adjacent 1 m × 1 m quadrats were set and surveyed. Sampling at NCO began on May 20, and samples were taken on the 10th, 20th, and 30th of each month until September 20. That is, a total of 13 samplings were performed each year. For each sampling, 15 1 m × 1 m quadrats (3 fences × 5 randomly established quadrats) were set and surveyed.

In each quadrat, all the species were recorded, and for each species, living aboveground tissue was clipped from individual plants and stored in a paper bag. Plant materials were dried in an oven (at 75°C, approximately 20 h), after which the weights were determined.

In this research, we used two databases. The first database is a short-term database of quadrat data collected from ERS and NCO for seven consecutive years (from 2014 to 2020). It is worth noting that at NCO, data are not available for 30 August 2018, 20 May 2020, and 20 June 2020. The second database is a long-term database for site ERS from 1984 to 2020 (except in 1995 and 1996, i.e., 35 years). The long-term database contains only data from two periods in the peak growing season in each year (i.e., sampling on August 15 and 30).

Meteorological data from ERS and NCO were obtained from the weather stations of IMGERS (http://www.cnern.org.cn) and the China Meteorological Data Service Centre (http://data.cma.cn), respectively. We considered 2017 to be a drought year because it was the second driest year in the past 15 years. Precipitation at ERS and NCO was equivalent to only 84.05% (278.2 mm) and 61.31% (168.6 mm) of the annual average precipitation, respectively ( Figure S1 ). We chose eight meteorological indices to analyze the impact of meteorological conditions on community structure. These eight indicators were classified into two categories. The first category included four past meteorological indicators (i.e., those before the growing season), namely the average temperature of the last coldest month (i.e., January, T _Jan), extremely low temperature in January (T _min), precipitation of the past years (P _past), and precipitation of the years with less precipitation in the last 2 years (P _less2). The second category included four current meteorological indicators (i.e., those of the current growing season), namely the average temperature of the last 15 days before sampling (T ₁₅), current accumulated temperature >0°C (from 1 January until the sampling day, CAT0), accumulated precipitation of the last 15 days before sampling (P ₁₅), and current total accumulated precipitation (from 1 January until the sampling day, CAP).

We built a phylogenetic tree using all vascular plants that were recorded at both sites. The phylogenetic tree was constructed in R package V.PhyloMaker (Jin and Qian, 2019). V.PhyloMaker uses GBOTB (i.e., GenBank taxa with a backbone provided by Open Tree of Life) as a backbone to generate a phylogeny (Smith and Brown, 2018). GBOTB includes all accepted families of vascular plants (according to APG IV) in the world. The species that were absent from the mega-tree were bound to the halfway point of the family (or genus) branch (by setting “Scenario = S3” in function phylo.maker) (Jin and Qian, 2019; Qian et al., 2019).

We used the net relatedness index (NRI), weighted by relative dry aboveground biomass of individual species (i.e., set “abundance.weighted = TRUE”), as an index of community phylogenetic structure. Relative dry aboveground biomass (RDAB) of species i was calculated as follows:

The NRI was calculated as follows:

where NRI is the net relatedness index, MPD_obs is the observed mean phylogenetic distance (MPD) in the quadrat, MPD_null are 999 random MPD values under the null model, and sd (MPD_null) is the standard deviation of MPD_null (Webb et al., 2002). The NRI was calculated in package “picante” (NRI is equivalent to −1 times the output of “ses.mpd”; see Kembel et al., 2010). A positive NRI indicates phylogenetic clustering (i.e., taxa are more phylogenetically similar than expected by chance), while a negative NRI indicates phylogenetic overdispersion (i.e., taxa are more phylogenetically distant than expected by chance) (Swenson et al., 2006).

All statistical analyses were performed using R 3.5.3 (R CoreTeam, 2019). We used one-way ANOVA followed by Duncan’s multiple range test (significance level p < 0.05) (De Mendiburu, 2014) to determine the significance of intra- and interannual differences in the NRI. If the NRI in the peak growing season was significantly lower than at the beginning and end of the growing season, habitat filtering would play a greater role in community assembly at the beginning and end of the growing season than in the peak growing season and, thus, lead to the intra-annual dynamic of the phylogenetic structure. We used multiple linear regression and structural equation models (SEMs) to quantify the relative contributions of meteorological factors before and during the current growing season to short-term community phylogenetic structure at both sites. The best multiple linear regression model was identified as that with the minimum Akaike information criterion (AIC). The SEM was fitted by maximum likelihood estimation using the “lavaan” package (Rosseel, 2012). An acceptable SEM fit was defined as a non-significant χ ² test result (p > 0.05), high goodness-of-fit index (GFI > 0.95), low standardized root mean square residual (SRMSR < 0.08), and low root mean square error of approximation (RMSEA < 0.05) (Schermelleh-Engel et al., 2003).

We also used yearly NRI anomalies (i.e., NRI of the last year – NRI of the current year, NRI_an) to reflect the influence of meteorological factors on community phylogenetic structure under long time series. A negative NRI_an indicates a declining NRI when compared to that of the last year, and a positive NRI_an indicates an increasing NRI when compared to that of the last year. We evaluated the effects of meteorological factors on NRI_an with classification and regression trees (CARTs) using the “rpart” package (Therneau et al., 2015). To avoid overfitting, the complexity parameter (cp) was set according to the smallest cross-validation error (xerror). In addition, we used linear mixed-effects models (LMMs) to evaluate the relationships between meteorological factors and NRI_an, incorporating potential meteorological factors as fixed effects and year as a random intercept. The LMMs were fitted using the “lme4” package (Bates et al., 2014), and the variance explained by the LMMs was calculated using the “MuMIn” package (Barton, 2020).

Eighty-eight species belonging to 29 families and 64 genera were recorded in this research. The species richness at ERS was significantly higher than that at NCO (75 species at ERS and 39 species at NCO, respectively, Figure 2 and Figure S2 ). In addition, the NRI was significantly greater than 0 at both sites, indicating phylogenetic clustering. The average NRI at ERS was significantly lower than that at NCO. We also calculated the average NRI and species richness in the peak growing season (i.e., in late August). We found that at site ERS, in the peak growing season, the NRI and species richness were significantly lower and higher than the annual mean values, respectively. In contrast, the difference in the NRI and species richness between the peak growing season and annual mean at the NCO site was not significant ( Figure 3 ).

At ERS, the intra-annual NRI exhibited an inverse unimodal curve (except in 2016 and 2017); that is, the NRIs were the highest at the beginning of the growing season and tended to decrease as the season progressed but then increased again at the end of the growing season (only significant in 2014 and 2018, Figure S3 ). In 2016, the NRIs tended to increase from the beginning of the growing season until the middle of the growing season and became inversely unimodal ( Figure 4 ). At NCO, the NRIs did not show a clear intra-annual trend. Only in 2014, 2016, and 2017, the NRIs were higher in the late-middle stages of the growing season than in the early-middle stages of the growing season ( Figure 5 and Figure S4 ).

The community phylogenetic structure at both sites showed significant interannual variation based on one-way ANOVA. At ERS, the annual mean NRIs showed an inverse unimodal curve from 2014 to 2020. The annual mean NRIs in 2018 (i.e., the year following the drought year) were the lowest, followed by those in 2019 ( Figure 6 ). The annual mean NRIs were the highest in 2014. At NCO, the lowest annual mean NRI values were also observed in 2018 and 2019, and the highest annual mean NRI values were observed in 2015 and then in 2016 ( Figure 7 ).

Multiple regression showed that meteorological factors before the current growing season accounted for 71.12% of the variation in the NRI at ERS. Among these factors, the average temperature in January (T _Jan) and lower precipitation in the last 2 years (P _less2) accounted for 15.53% and 52.20% of the variation, respectively. In addition, the interaction between T _min and P _less2 accounted for 1.98% of the variation in the NRI. Meanwhile, the meteorological factors of the current growing season (i.e., CAP, CAT0, and T ₁₅) accounted for only 7.48% of the variation in the NRI ( Table 1 ).

Multiple regression showed similar results at NCO. Meteorological factors before the current growing season (i.e., T _Jan and P _less2) accounted for 61.20% of the variation in the NRI. Current accumulated precipitation (CAP) and current accumulated temperature >0°C (CAT0) accounted for 1.23% and 1.09% of the variation in the NRI, respectively ( Table 2 ).

Our SEMs fit the meteorological data well (χ ² = 2.987, p = 0.394, GFI = 0.985, SRMR = 0.020, RMSEA < 0.001 and χ ² = 0.852, p = 0.837, GFI = 0.997, SRMR = 0.017, RMSEA < 0.001 at ERS and NCO, respectively). In general, meteorological factors had only a slight direct impact on the phylogenetic structure at both sites. However, meteorological factors showed a strong indirect effect on the NRI via the proportion of perennial grass biomass (P.Grass).

At ERS, P _less2 and T _min had positive and negative direct effects on the NRI value, respectively (standardized path coefficient, β = 0.10, p = 0.03 and β = −0.05, p = 0.24, respectively). Notably, P _less2 and T _min had significant positive direct effects on the proportion of perennial grass biomass in the community (i.e., P.Grass, β = 0.68 and 0.46, respectively). Meteorological conditions of the current growing season (Current.T, including CAT0 and T ₁₅) had non-significant direct negative effects on the NRI value (β = −0.09, p = 0.09). In addition, the indirect effects of P _lesst2, T _Jan, and Current.T via P.Grass on the NRI were 0.60 (p < 0.01), 0.40 (p < 0.01), and −0.26 (p = 0.02), respectively ( Figure 8A ).

At NCO, P _less2 and T _min had non-significant positive and negative direct effects on the NRI value (β = 0.10, p = 0.74 and β = −0.08, p = 0.74, respectively). Meanwhile, Current.T had a non-significant direct effect on the NRI value (β = 0.23, p = 0.72). The indirect effect of P _less2 on the NRI was 0.57 (p = 0.04). However, the indirect effects of T _min and Current.T on the NRI were not significant (β = −0.19, p = 0.55 and β = 0.36, p = 0.67, respectively) ( Figure 8B ).

The CART results showed that only meteorological factors before the growing season (i.e., T _min, P _less2, and P _past) were used in tree construction, and NRI_an was mainly related to T _min and P _less2 ( Figure 9 ). Specifically, NRI_an was divided into two branches by a T _min of −38°C; that is, T _min <−38°C resulted in a negative NRI_an (node 2, mean NRI_an = −0.89), and T _min ≥−38°C resulted in a positive NRI_an (node 3, mean NRI_an = 0.44). In addition, the tree was further divided into five branches according to P _less2 and P _past. Notably, the interaction between T _min <−38°C and P _less2 <289 mm resulted in the lowest NRI_an (mean = −2.45, n = 8), followed by the interaction between T _min ≥−38°C and P _past <280 mm (mean = −0.25, n = 8).

The LMM results showed that only T _min and the random factor (i.e., year) were significantly correlated with NRI_an (p = 0.01, Figure S6C ). LMMs also showed that the amount of variation in NRI_an explained by the random factor (i.e., year) was greater than that explained by fixed factors (e.g., T _min) ( Table 3 ).

At both sites, the NRI was significantly greater than 0, indicating phylogenetic clustering (Webb et al., 2002). There are two processes that may result in phylogenetic clustering (under the assumption of niche conservatism) (Webb et al., 2002; Donoghue, 2008; Mayfield and Levine, 2010), habitat filtering (i.e., abiotic filtering), and competitive exclusion between distantly related phylogenetic clades because of competitive ability differences (i.e., biotic filtering) (Mayfield and Levine, 2010; Goberna et al., 2015), although it is difficult to detect abiotic and biotic filtering by phylogenetic analyses alone (Swenson and Enquist, 2009). Our results showed that the NRI was highly positively correlated with the proportion of perennial grass ( Figures S5 , S6A ). This suggests that the relative competitive advantage of grasses over forbs and shrubs is an important driver of community phylogenetic clustering. Meanwhile, the community phylogenetic structures of both sites were related to meteorological factors, especially considering that the degree of clustering of community phylogenetic structure in early spring and autumn was higher ( Figures 4 , 5 ; see below), suggesting that abiotic filtering factors may also have an impact on community phylogenetic structure. However, given that life strategies and competitive abilities are strictly related to the environment, abiotic filtering factors may still be the major driver of community phylogenetic clustering (Saito et al., 2016). Our results were partly consistent with those of Yang et al. (2018), who found that in an undisturbed (i.e., control treatment) Inner Mongolian grassland, the community showed a positive NRI (although not significantly different from 0), and nitrogen addition promoted the colonization of species distantly related to the residents. Qian et al. (2019) also found that, in general, the angiosperm community showed phylogenetic clustering in northwestern China.

In general, the intra-annual NRI showed an inverse unimodal curve at ERS ( Figure 4 ). This is consistent with our hypothesis 1 that habitat filtering would play a greater role in community assembly at the beginning of the growing season than in the peak growing season. The reason may be related to the harsh environment at the beginning and end of the growing season. To avoid facing low temperature and precipitation, some species (e.g., Heteropappus altaicus) delay revival at the beginning of the growing season. For the same reason, some species (especially annual plants, e.g., Dontostemon dentatus) complete their life cycle before the autumn arrives. That is, competition between species is most intense during the peak growing season (i.e., July to August) (Doležal et al., 2019) when the environmental filtering effect is the smallest. As a result, the community showed less phylogenetic clustering during the peak growing season.

Notably, both the intra- and interannual variability (indicated by coefficients of variation, CVs) at ERS were greater than those at NCO ( Figures 4 – 7 ). This result indicates that the community phylogenetic structure at NCO is more stable than that at ERS, which has more species richness than the former ( Figure 3B ). The relationship between biodiversity and ecosystem stability remains a controversial topic (Fonseca and Ganade, 2001; Mazzochini et al., 2019). Traditionally, more ecologists believe that biodiversity is positively related to ecosystem stability (i.e., resistance and resilience) because larger species pools can contain more species and functional groups that adapt to different disturbances (McNaughton, 1978; Yachi and Loreau, 1999). In recent years, increasing evidence has shown that biodiversity may also be negatively correlated with ecosystem stability. For instance, Pfisterer and Schmid (2002) found that species-poor systems were more resistant to disturbance than species-rich systems. Van Ruijven and Berendse (2010) found that biodiversity enhances grassland community recovery, but its relationship with resistance is highly dependent on biomass. However, previous studies mainly used the dynamics of biomass as an indicator of community resistance and resilience to perturbation (Tilman and Downing, 1994; Tilman et al., 2006; Isbell et al., 2015). Here, we provide a new perspective on how biodiversity affects ecosystem stability. That is, in view of community phylogenetic structure, we provide evidence that communities with low biodiversity may have higher stability than communities with high biodiversity.

Our results showed that meteorological factors before the growing season explained more of the variation in the NRI and NRI_an than meteorological factors during the current growing season on both short-term and long-term scales ( Tables 1 – 3 , Figures 8 , 9 ). Vegetation and community production show a time-lag response to climate that has been widely reported (Bradley et al., 1999; Piao et al., 2003; Bai et al., 2004; Vicente-Serrano et al., 2013; Wu et al., 2015); however, the lag response of community phylogenetic structure still receives little attention. Here, we provide evidence that the dynamics of grassland community phylogenetic structure are deeply influenced by meteorological conditions before the growing season rather than during the current growing season, which suggests the time lag of the community phylogenetic structure responses to meteorological factors (i.e., hypothesis 2).

Our results highlight the fact that extreme climatic events, such as drought and extreme cold, may impact grassland community structure and reduce the degree of phylogenetic clustering. Additionally, the SEMs showed that meteorological factors mainly indirectly changed community phylogenetic structure by influencing the proportion of perennial grass biomass ( Figure 8 ). As we mentioned above, the competitive advantage of grasses is one of the key reasons why the community remains phylogenetically clustered. In other words, drought and extreme cold inhibit the dominance of grass ( Figure 8 ). We found that after drought years (i.e., in 2018 and 2019), the lower NRIs at ERS were mainly related to an abnormal decrease in perennial grasses and an increase in forbs (e.g., Potentilla bifurca and Klasea centauroides) in the community. In contrast, due to the lack of perennial forbs, niches were occupied by annual forbs (e.g., Salsola collina and Dysphania aristata) at NCO. This was consistent with the finding of VanderWeide and Hartnett (2015) that drought decreased grass but not sedge and forb bud and shoot densities. Bai et al. (2004) found a negative correlation between dominant grass species and non-dominant species in a typical steppe, which is a key mechanism maintaining the stability of the ecosystem. Drought or extreme cold may decrease the proportion of perennial grass because grasses and forbs have different reproduction strategies. Specifically, compared to grasses, forbs rely more on sexual reproduction (Song et al., 2001; Stampfli and Zeiter, 2004), and seeds are more likely to survive low temperatures and in other unfavorable environments than buds (Carter et al., 2012). Stampfli and Zeiter (2004) suggested that recolonization from seeds was the key to recovery after interference. However, the sexual reproduction ability of dominant grasses in the Mongolian Plateau (i.e., S. grandis and L. chinensis) is usually lacking (Teng et al., 2006). Therefore, extreme weather events are more likely to lead to the death of perennial grasses. In addition, plants tend to allocate more resources to sexual reproduction under resource shortages or in unfavorable habitats (Olejniczak, 2001; Dong et al., 2015). As a result, extreme meteorological events before the growing season inhibited the reproduction of grasses and further reduced the degree of phylogenetic clustering.

Our results were consistent with those of Evans et al. (2011), who found that drought may result in a reduction in dominant grass and correspond with an increase in forbs. The authors attributed this to the higher colonization ability of forbs under drought conditions (Evans et al., 2011). Other researchers have suggested that grass generally has shallow and fibrous roots and thus is more sensitive to drought than forbs with deep roots and taproots (Buckland et al., 1997; Dunnett et al., 1998). In addition, in Mediterranean grasslands, researchers found that both dominant perennial grasses and annuals increased substantially after extreme events (such as burning), but the latter increased even faster for the reason of surviving as seeds (Vidaller et al., 2019). As a result, when faced with an unfavorable environment, perennial grasses in grasslands may be more affected, while forbs and annuals occupy the niche vacancies of perennial grasses, thus changing community species composition and affecting community phylogenetic structure.

Notably, at both sites, the NRIs in 2018 and 2019 were significantly lower than those in other years ( Figures 6 , 7 ). Multiple regression showed that precipitation of the years with less precipitation in the past 2 years (i.e., P _less2) explained more variation than precipitation of the last year ( Tables S1 , S2 ). This implies that the effects of drought last more than 1 year. Our results were consistent with those of VanderWeide et al. (2014), who found that the impact of drought on grasslands may last more than 1 year. Oesterheld et al. (2001) also found that interannual variation in primary production of a semiarid grassland was related to previous-year production. The authors suggested that seed, plant, and tiller density might be reduced after a drought year with low precipitation and, thus, might either constrain the ability of ecosystems to respond to a subsequent increase in water availability (Oesterheld et al., 2001).

We found that extremely low temperature (i.e., T _Jan and T _min) contributed to the reduction in phylogenetic clustering, especially on a long time scale ( Figure 9 and Table 3 ). However, extremely low temperature had a significant influence on NRI_an rather than the NRI, suggesting that community phylogenetic structure depended on not only meteorological conditions but also the community structure of the last year. LMM analyses additionally showed that the year explained more of the variation in NRI_an than meteorological factors (the random factor explained 0.77 and 0.82 of the variation in NRI_an at ERS and NCO, respectively). That is, the interannual dynamics of community phylogenetic structure are obvious regardless of changes in climatic conditions. Therefore, grassland community phylogenetic structure should not be treated as quiescent, and assessment of specific grassland communities’ phylogenetic structure needs to include a multiyear average rather than one sampling; in particular, the impact of extreme meteorological events before sampling needs to be fully considered.

Our results were partially consistent with those of Zheng et al. (2019), who found that winter temperature was a major factor leading to phylogenetic overdispersion of the shrub community on the Mongolian Plateau. Our results suggested that extremely low temperatures play an important role in steppe biodiversity by inhibiting the competitive ability of grass. Therefore, with global climate warming, the reduction in extreme low-temperature events (i.e., frost days; see Sillmann et al., 2013) may adversely affect grassland biodiversity. The impact of extreme climate conditions on grassland community structure still needs more exploration.