*4.1. Study Site and Experimental Design*

Our study site (40◦100 45" N, 116◦260 13" E, and 50 m above sea level) is located at the Station of the Institute of Grassland, Flowers and Ecology on Xiaotang Mountain in Beijing, China. The mean annual precipitation is 526 mm, the mean annual temperature is 11.8 ◦C (2000–2018) in our study site. The common garden experiment was established in 2019, using a split-plot experiment design with N addition (0 and 6 g N m−<sup>2</sup> year−<sup>1</sup> ) as the main treatment factors and plant species richness (species richness is 1, 2, 4, 6, and 8) as the subfactors (Figure S2 and Table S1). The experiment comprised the following: 2 N addition levels × 4 replicate blocks × 17 combinations of species per N treatment/replicate (four combinations per level of species richness with 1, 2, 4, and 6 species and one with 8 species) = 136 pipes. Each pipe (30 cm in diameter and 50 cm in height polyvinylchloride sealed pipes) was filled with uniformly mixed soil (natural soil with stones and roots sieved out) with sand (in a 3:1 soil–sand ratio) and then buried into the ground without any shelters. Eight plant species were selected in our study to build assembled grassland communities [48]. Different species of seeds were mixed up well and distributed randomly in the pipes, maintained at around 60 individuals in each pipe at the beginning of the experiment, and did not re-sow later. Nitrogen was added as urea (12.86 g m−<sup>2</sup> year−<sup>1</sup> ), dissolved in N-free water, and then applied by spraying on 1 May of each year; the control treatment received equal N-free water. Additional information regarding the experimental design was provided by Wang et al. [48].

#### *4.2. Phenology Monitoring*

To track the flowering phenology of *M. sativa*, phenology was monitored every 3–4 days during the growing season from May to September 2019. Three individuals in each pipe were randomly selected, marked, and monitored across the growing season. The first and last date of a flower observed for each marked individual was recorded as the first and last flowering day, and the periods between the first and last flowering day were recorded as flowering duration. Flowering numbers were counted for each marked individual. Flowering phenological events and flowering numbers were averaged for three individuals in each pipe.

#### *4.3. Functional Traits and Abiotic Factors Measurements*

Light acquisition traits (plant height and relative height, leaf mass and area, leaf length and width, and specific leaf area) and nutrient acquisition traits (leaf carbon content, leaf nitrogen content, leaf C/N ratio, abundance and relative abundance of plant species) are closely related to plant phenology [49,50]. Consequently, we determined these traits to explore the mechanism underlying regulating the response of flowering phenology to experimental N addition and plant species richness in the peak growing season in 2019. Before the measurements, we investigated the abundance and height of each plant species in the pipes. *M. sativa* is the predominant species (relative abundance >40% in each pipe) (Figure S3), three healthy individuals, which we measure flowering phenology were selected to measure the species-level traits in each pipe, and six leaves on each individual were selected to measure leaf traits once at full flowering, and the leaves in each pipe were collected on the same day. The traits were quantified using standard methods proposed by Pérez-Harguindeguy et al. [51]. Specific leaf area was calculated as the ratio of leaf area to its dry weight. Leaf area, length, width, and maximum width, spread leaves were scanned and analyzed by Li-Cor 310 (Li-Cor Inc., Lincoln, USA), and then leaves were oven-dried to a constant weight. The oven-dried leaf samples were ground to determine leaf C and N concentration with an elemental analyzer (PE 2400 II, PerkinElmer Ltd., CT, USA) and then to calculate the C/N ratio. To measure the biomass of *M. sativa*, the aboveground part of each pipe was clipped in early September (the peak of the growing season) in 2019. *M. sativa* clipped from each pipe were pooled together and then oven-dried at 65 ◦C to a constant weight.

Soil temperature and moisture at a depth of 10 cm were measured every week from April to October with a W. E. T sensor kit (Delta-T Devices Ltd, Cambridge, UK). Three soil cores were collected in each pipe in early September at a depth of 10 cm and then mixed into one sample. Available-soil N (Ammonium (NH<sup>4</sup> + ) and nitrate (NO<sup>3</sup> −)) concentrations in the extracts were determined calorimetrically by automated segmented flow analysis (Bran + Luebbe AAIII, Bran + Luebbe Ltd, Hamburg, Germany).

#### *4.4. Statistical Analyses*

We analyzed experimental data with the following three steps. First, we scaled the species-level height to the community level by calculating the mean of abundance distributions (Equation (1) [52]):

$$Mean\_{\mathbb{C}} = \sum\_{i}^{n} p\_{i} T\_{i} \tag{1}$$

where *p<sup>i</sup>* and *T<sup>i</sup>* are the relative abundances and the plant height of the species *i*, respectively, and *n* is the number of plant species. Hence, the average height of *M. sativa* divided by *Mean<sup>c</sup>* is the relative height.

Second, we applied linear mixed-effects models using the "*lme*" function (package "*nlme*" [53]) to test the effects of N addition and plant species richness on soil temperature and moisture. We set N addition and species richness levels as fixed effects; the date, block, and plant combination were set as random effects in each model to account for variation among repeated measurements. In addition, linear mixed-effects models were also used to examine the effect of N addition and plant species richness on flowering phenology and functional traits. Nitrogen addition and species richness levels were treated as fixed effects, and the block and plant combinations were treated as random effects.

Third, variation partitioning analysis that partitioned the variance shared by all factors was then used to quantify the unique contribution of biotic and abiotic factors. Structural equation modeling was employed to evaluate which are the major factors that influence flowering phenology [54] by the package '*piecewise-SEM*' in R software [55]. The SEM requires establishing an a priori framework based on the hypothesized causal relationships among these variables. Second, the relationships between these variables were examined by bivariate correlations. Finally, models with lower *Fisher's C* and Akaike information criterion (*AIC*) and higher *p* values (*p* ≥ 0.05) were selected in our analysis (Figure 5). All statistical analyses and graphs were prepared in R 3.2.2 [56]. Differences were considered to be statistically significant at *p* ≤ 0.05.

#### **5. Conclusions**

The study highlighted the influence of functional traits on flowering phenology following nitrogen addition levels and plant species richness gradients in an assemblage grassland through a common garden experiment. It was observed that the first flowering day was delayed 0.31 days, the last flowering day advanced 0.64 days, and the flowering duration was shortened by 0.95 days with per-plant species increase, but the effects of plant species richness on flowering phenology did not interact with nitrogen addition, which indicates that nitrogen addition could change plant flowering phenology by changing biodiversity, but the effects would be independent with the effects of biodiversity. Moreover, flowering phenology changed following nitrogen addition levels, and plant species richness gradients were mainly driven by the intraspecific variation in functional traits, which suggests that variation in functional traits among communities may be a good predictor for the dynamic of plant phenology under global changes.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants12101994/s1. Figure S1: Design and photo of our experiments; Figure S2: Photos of experimental arrangement; Figure S3: The average biomass, relative biomass, relative abundance, and relative height of *M. Sativa* along the plant species richness gradients under different nitrogen (N) addition levels; Figure S4: Pearson's correlation between the biotic and abiotic factors; Figure S5: Relative contributions of light acquisition traits, nutrient acquisition traits, and abiotic factors to flowering events; Figure S6: Original structural equation modellings of N addition and plant species richness on the first flowering day, last flowering day, flowering duration, and flower number. Table S1: Assemblage types at different species richness levels.

**Author Contributions:** Conceptualization, C.W. and X.Z.; methodology, C.W., Y.M., X.L. and Y.H.; writing—original draft preparation, Y.M. and C.W.; writing—review and editing, Y.M., X.Z. and C.W.; funding acquisition, X.Z., C.W. and Y.H.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the Beijing Natural Science Foundation (Grant No. 5232006), National Natural Science Foundation of China (Grant No. 31901173), China Key Technologies Research and Development Program (Grant No. 2020YFD1000201), and the Excellent Talents Innovation Project of Shanxi, China (Grant No. 201805D211018).

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** We thank Wenjun Teng, Chao Han, and Ruibin Xue for their help with the fieldwork.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
