2 | Materials and Methods

2.1 | Study site and experimental design

Our study site is located at the field research station of Beijing Research & Development Center for Grasses and Environment (40°10’ 45” N, 116°26’13” E) on Xiaotang Mountain in Beijing, China. The mean annual temperature is 11.8 oC (2000-2018). The mean annual precipitation is 526 mm (2000-2018), with >80% of the precipitation falling in the growing season (June-September). Soil organic carbon of top soil (0-10 cm) is around 1%, and the average pH value is 7.47. The study area is dominated by perennial plants, including E. nutans and T. mongolicum .
Our experimental plots were established within an area of 8×12 m in 2019, using a randomized block design with nitrogen (N) addition and plant diversity gradient as main treatment factors. Each block contained twelve treatments, crossing two levels of N addition (no nitrogen addition and add 60 kg N ha-1 yr-1as urea) and six levels of plant diversity (0, 1, 2, 4, 6, and 8). Each diversity level had four plant assemblage types, except for 0 and 8 (Supplementary Table 1). Each treatment had four replicates, resulting in 144 polyvinylchlorid bottom sealed pipes (Supplementary Fig. 1b). Each pipe is 30 cm in diameter and 50 cm in height, three holes were drilled into the bottom to provide drainage, filled with uniformly mixed soil and sand (soil-sand ratio is 3:1), and then buried into soil (Supplementary Fig. 2). Soil inside and outside of pipe was 5 cm lower than its upper edge. According to the germination rate of each species (Supplementary Table 1), seeds was sown at the depth of 1 cm in each plot after fully watered in March, and maintained 120 seedlings in each plot. The eight species we selected were Poa annua , Carex breviculmis , Medicago sativa , Astragalus adsurgens Pall ,Dianthus barbatus , Penstemon campanulatus ,Chrysanthemum maximum , and Allium schoenoprasum , the assemblage types at each diversity level were showed in the Supplementary Table 1. The amount of N addition is twice the background deposition in Beijing (Yu et al., 2019), and urea was applied by spraying on 1 May.

2.2 | Phenology monitoring

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

2.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 (C) and nitrogen (N) content, leaf C/N ratio, biomass and abundance, relative biomass and abundance) are closely related to plant phenology (Lavorel & Grigulis, 2012; Grigulis et al., 2013). Consequently, we determined these traits to explore the mechanisms underlying regulating the response of flowering phenology to experimental N addition and plant diversity gradients. Before the measurements, we investigated the abundance and height of each species in the plot. M. sativa is the predominant species (relative abundance >40% in each pine) and 6 healthy mature individuals were selected to measure the species-level traits in each plot. The functional traits were quantified using standard methods proposed by Pérez-Harguindeguy et al. (2013). The specific leaf area was calculated as the ratio of leaf area to its dry weight. To measure leaf area, length, width, and maximum width, spread leaves were scanned and analyzed Li-Cor 310 (Li-Cor Inc. USA), and then leaves were oven-dried to a constant weight. Finally, the oven-dried leaf samples were ground to determine leaf carbon and nitrogen with an elemental analyzer (PE 2400 II, USA). To measure the biomass of community and individual species, the aboveground part of each plot was clipped in early September (the peak of growing season). Plants clipped from each plot were pooled together, sorted to species, and then oven-dried to a constant weight.
Soil temperature and moisture at the depth of 10 cm were measured every week from April to October with W. E. T sensor kit (Delta-T Devices Ltd, UK). Three soil cores were collected in each plot in early September at the depth of 10 cm, and then mixed together into one sample. Available soil N (Ammonium (NH4+) and nitrate (NO3-)) contents in the extracts were determined colorimetrically by automated segmented flow analysis (Bran + Luebbe AAIII, Germany)

2.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 the abundance distributions (Equation 1, Gross et al., 2009):
\(\text{Mean}_{j}=\sum_{i}^{n}{p_{i}T_{i}}\) (1)
where \(p_{i}\) and \(T_{i}\) are the relative abundance and the plant height of the species j , respectively, and n is the number of species.
Second, we examined how N addition and plant diversity affected environmental factors, the functional traits and flowering events ofM. sativa . We applied linear mixed effects models using “lme4” function (package “NLME”, Pinheiro et al ., 2007) to test the effects of N addition and plant diversity loss on soil temperature and moisture separately in 2019. We set treatments as fixed effects, block and time as a random effect in each model to account for variation among repeated measurements of temperature or moisture. Linear mixed effects models were also used to examine the effect of N addition and plant diversity loss on flowering phenology (first and last flowering day, flowering duration, and flower numbers) and functional traits (leaf and community traits). Treatment was treated as fixed effects, and block was treated as a random effect to account for variation within block.
Third, partial correlation was conducted to evaluate the relationships between the flowering events and the various factors (Chen et al., 2019). For example, after controlling N addition and plant diversity levels, we examined the relationships of the flowering events with light acquisition traits, nutrient acquisition traits, and abiotic factors. Variation partitioning analysis that partitioned the variance shared by all factors was then used to quantify the unique contribution of each group of factors (Chen et al., 2019). Structural equation modelling analysis was employed to evaluate the hypothesized underlying factors that influence flowering phenology (Wang & Tang, 2019b) using the package ‘piecewise-SEM ’ in R software (Shipley, 2000). The model was assessed by Fish C statistic, Akaike information criterion (AIC), and P values.
All statistical analyses and graphs were prepared in R 3.2.2 (R Core Team, 2018). Differences were considered to be statistically significant at P ≤ 0.05.