The study area is located at Guizhou University in Huaxi District, Guiyang, Guizhou Province (106° 39 ′ E, 26° 27 ′ N), and has an altitude of 1,100 m. It has a subtropical monsoon humid climate, with obvious plateau climate characteristics. The annual average temperature is 15.3°C, the annual rainfall is 1,129.5 mm, the frost-free period is 270 d, while the annual sunshine duration is 1,148.3 h. The soil type is yellow soil. The physical and chemical properties of soil are provided below ( Table S1 ). Cultivated grasslands of Medicago sativa (cv. 601) and Dactylis glomerata (cv. Amba) were established in October 2018. The effects of different fertilization gradients on the production characteristics and nutritional quality were explored from May, 2019 to October, 2022.

Medicago sativa (cv. 601) and Dactylis glomerata (cv. Amba), provided by Guizhou Zhongzhiheng Ecological Technology Co., Ltd., were used as the experimental materials. Cultivated grasslands of the two varieties were established in October 2018, with a sowing rate of 22.5 kg/hm². The plot size was 2.0 m × 2.0 m, and a drill sowing method was employed. Each plot consisted of six rows, with a row spacing of 30 cm. At the same time, five fertilization gradients (0, 10, 20, 30, and 40 t/hm²) were set up in a randomized block design. Organic fertilizer was evenly applied on the soil surface of each plot, and the application was repeated three times. There were 15 plots in each of the cultivated grasslands of Medicago sativa and Dactylis glomerata. Before sowing, the land was ploughed and prepared, and weeds were removed. After sowing, 2–3 cm of covering soil was provided. Timely field management was ensured throughout the growth period. A sheep manure-fermented organic fertilizer (pH value: 7.4; dry basis nutrient contents: N=2.23%, P₂O_{5 =} 1.775%, K₂O=3.91%) was used. Fertilization was done during the early stages of the growing season each year (March 20). Meanwhile, production performance indicators were surveyed during the initial flowering stage of Medicago sativa and the heading stage of Dactylis glomerata from 2019 to 2022. A 1 m² quadrat was selected in each plot. The survey quadrat was controlled to be 50 cm away from the plot edge to avoid edge effects, and the height (cm), coverage (%), and biomass (g/m²) of each species in the quadrat were measured. At maturity, the aboveground parts were cut, brought to the laboratory in bags, dried to a constant weight at 80°C and weighed (g/m²). The survey was conducted 5–6 times a year based on growth status. In June, 2019, 2020, 2021 and 2022, the aboveground parts of Medicago sativa and Dactylis glomerata were harvested, packed in envelopes, brought back to the laboratory, dried at 105°C for 30 min, adjusted to 65°C and dried to a constant weight, and then crushed in a pulverizer and screened with a 40-mesh sieve to determine forage nutrition.

Eight main indicators were measured, including plant height, coverage, dry/fresh weight ratio, dry matter yield, crude protein (CP), ether extract (EE), neutral detergent fiber (NDF), and acidic detergent fiber (ADF). The measurement methods were as follows:

Plant height: Randomly select ten plants from each quadrat during sampling, measure the distance from the ground to the highest part of each plant with a tape measure, and calculate the average value as plant height.

Coverage: Record the coverage of forage within each quadrant through visual observation.

Dry/fresh weight ratio and dry matter yield: Immediately weigh fresh forage yield after each cutting, and bring fresh forage to the laboratory in bags. Dry it to a constant weight at 80°C, weigh the dry weight after cooling, and calculate dry/fresh weight ratio.

CP: Kjeldahl nitrogen determination (Jee, 1995).

EE: Soxhlet extraction method (Han et al., 2022), ANKOMAXT15i automatic fat analyzer.

NDF and ADF: measured by the filter bag technology (Mcroberts and Cherney, 2014), using an ANKOM A200i semi-automatic fiber analyzer.

Preliminary data analyses were performed using Office Excel 2019 and one-way ANOVA, while multiple comparisons were conducted using the SPSS 26.0 software. Two-way ANOVA was used to test the significance of the effects of year and fertilization on the production indicators of the two types of cultivated grasslands. Origin 2022 was adopted for drawing. The subordinate function comprehensive evaluation method was employed to comprehensively evaluate the plant heights, coverages, dry/fresh weight ratios, and the grassland yields.

The formula for the subordinate function (Liu et al., 2020) is:

where U(X_ij) denotes the subordinate value of the j^th indicator under the i^th fertilization gradient; X_ij represents the measured value of the j^th indicator under the i^th fertilization gradient; X_jmax and X_jmin denotes the maximum and minimum values of the j^th indicator under all fertilization gradients, respectively;U(X_ij) denotes the average subordinate value of all indicators under the i^th fertilization gradient, that is, the larger the average value, the better the production performance.

The grasslands of M. sativa and D. glomerata were surveyed 5–6 times per year from 2019 to 2022. The dry matter yield of M. sativa presented a trend of increasing first and decreasing afterward with increasing fertilization gradient, except in November 2019, September 2020, July 2021, and September 2022; and peaked under the fertilization gradient of 20 t/hm², except in June 2019 ( Figure 1 ). The seasonal changes in the maximum dry matter yields of M. sativa under fertilization also presented a trend of increasing first and decreasing afterward, except in 2020 ( Figure 1 ). The maximum values occurred in July or August each year, and were 2.70, 3.21, and 3.08 t/hm², respectively.

The dry matter yield of D. glomerata showed a trend of gradual increase with increasing fertilization gradient, except in 2019 and 2020; and peaked under the fertilization gradient of 40 t/hm² ( Figure 1 ). In 2019, the maximum dry matter yields of D. glomerata under fertilization showed a trend of gradual decrease over time. Conversely, in 2020, the maximum dry matter yields of D. glomerata showed a gradual increasing trend. In 2021 and 2022, the maximum dry matter yields of Dactylis glomerata presented a trend of increasing first and decreasing afterwards. The maximum values occurred in August, and were 1.17 and 0.61 t/hm², respectively ( Figure 1 ).The annual dry matter yield of M. sativa presented a trend of increasing first and decreasing afterward with increasing fertilization gradient in 2019–2022, and peaked under the fertilization gradient of 20 t/hm² ( Table 1 ). The maximum values were 8.14, 13.87, 11.85, and 13.08 t/hm², respectively, and the maximum values were 2.17–12.72 times the minimum values (under the fertilization gradient of 0 t/hm²). From 2019 to 2022, the annual dry matter yield of M. sativa presented a trend of increasing first and decreasing afterward under the same fertilization gradient. Under the fertilization gradients of 0 and 10 t/hm², the maximum annual dry matter yields of Medicago sativa all occurred in 2021, and were 5.47 and 6.94 t/hm², respectively. Under the fertilization gradients of 20, 30, and 40 t/hm², the annual dry matter yield of Medicago sativa showed a trend of gradual increase under the same fertilization gradient, and the maximum values all occurred in 2022 and were 13.08, 9.90, and 8.16 t/hm², respectively ( Table 1 ).

The annual dry matter yield of D. glomerata showed a trend of gradual increase with increasing fertilization gradient in 2019–2022, and peaked under the fertilization gradient of 40 t/hm². The maximum values were 3.29, 5.10, 3.53, and 2.29 t/hm², respectively, and the maximum values were 1.64–3.19 times of the minimum values (under the fertilization gradients of 0 and 10 t/hm²) ( Table 1 ).

Figure 2 shows that the annual average plant height of M. sativa showed a trend of gradual increase, except in 2019. Annual average plant height and coverage both presented a trend of increasing first and decreasing afterwards with increasing fertilization gradient within the four years, and both peaked under the fertilization gradient of 20 t/hm². The maximum annual plant heights were 1.38–2.16 times the minimum annual plant heights. The maximum annual coverages were 1.22–2.24 times the minimum annual coverages.

From 2019 to 2022, the annual average plant height of D. glomerata showed a gradual increase. The maximum annual plant heights (19.46, 20.16, 22.75, 24.85, and 25.01 cm) were 1.22–1.45 times the minimum annual plant heights. The annual average coverage of D. glomerata showed a trend of gradual increase with increasing fertilization gradient. The maximum annual coverages (34.47%, 40.53%, 49.4%, 55.00%, and 64.00%, respectively) of D. glomerata were 1.44–2.18 times the minimum annual coverages ( Figure 2 ).

Table 2 shows that the CP and EE contents of M. sativa both presented a trend of increasing first and decreasing afterward with increasing fertilization gradient, and peaked under the fertilization gradients of 30 and 10 t/hm², respectively. The maximum values (21.33% and 11.34%) were 1.42 and 1.43 times of the minimum values. The NDF and ADF contents of M. sativa presented “N-shaped” and “M-shaped” change trends with increasing fertilization gradient, and reached their minimum values under the fertilization gradients of 30 and 20 t/hm², respectively. Their minimum values (22.62% and 30.95%) were 72.80% and 94.13% of their maximum values, respectively.

Similar to M. sativa, the CP and EE contents of D. glomerata also presented a trend of increasing first and decreasing afterward with increasing fertilization gradient, and peaked under the fertilization gradients of 20 and 30 t/hm², respectively ( Table 2 ). The maximum values (18.04% and 3.76%) were 1.14 and 1.67 times of the minimum values. The NDF and ADF contents of D. glomerata presented “V-shaped” and “M-shaped” change trends with increasing fertilization gradient, and peaked under the fertilization gradients of 30 and 20 t/hm², respectively. Their minimum values (33.37% and 30.46%) were 92.34% and 95.55% of their maximum values, respectively.

According to the two-way ANOVA ( Table 3 ) of forage production performance indicators using fertilization and year, year, fertilization, and Y × F all significantly affected the plant height, coverage, dry/fresh ratio, and yield of M. sativa (P<0.05), except that Y × F had a non-significant effect on coverage. Year, fertilization, and year × fertilization (Y×F) all significantly affected the plant height, coverage, dry/fresh ratio, and yield of D. glomerata (P<0.05), except that Y × F had a non-significant effect on plant height (P>0.05).

According to the correlation analysis between the forage indicators of the cultivated grasslands of M. sativa and D. glomerata ( Figure 3 ), the plant height of M. sativa was significantly positively correlated with coverage, yield, and CP content and significantly negatively correlated with ADF content (P<0.05), but had no significant correlation with dry/fresh ratio, EE content, or NDF content (P>0.05). The coverage of M. sativa was significantly negatively correlated with dry/fresh ratio and ADF content and significantly positively correlated with yield and CP content (P<0.05). The EE content of M. sativa was significantly positively correlated with NDF content (P<0.05), but had no significant correlation with any other indicator (P>0.05) ( Figure 3A ).

The coverage, dry matter yield, CP content, and EE content of D. glomerata were significantly positively correlated (P<0.05). The plant height of D. glomerata was significantly negatively correlated with dry/fresh ratio and NDF content (P<0.05). The dry/fresh ratio of D. glomerata was significantly positively correlated with NDF content, and significantly negatively correlated with plant height, dry matter yield, CP, and EE (P<0.05). The ADF content of D. glomerata had no significant correlation with any other indicator (P>0.05) ( Figure 3B ).

Comprehensive evaluation of the subordinate function based on height, coverage, dry/fresh ratio, and yield ( Table 4 ) indicated that the average subordinate function values of M. sativa under the fertilization gradient of 20 t/hm² were all greater than those under other fertilization gradients. In 2019–2022, under the fertilization gradients of 0, 10, and 40 t/hm², the average subordinate function values of M. sativa all presented a trend of increasing first and decreasing afterward over time, and peaked in 2020, 2021, and 2021, with maximum values of 0.41, 0.55, and 0.70, respectively. Under the fertilization gradient of 20 t/hm², the average subordinate function values of M. sativa presented a trend of increasing first, decreasing afterward, and then increasing again over time, and peaked in 2022, with a maximum value of 0.83. Under the fertilization gradient of 30 t/hm², the average subordinate function values of M. sativa kept increasing over time, and peaked in 2022, with a maximum value of 0.76.

The average subordinate function values of D. glomerata under the 40 t/hm² fertilization gradient were all greater than those under other fertilization gradients. In 2019–2022, the average subordinate function values of D. glomerata all presented a trend of increasing first and decreasing afterward over time. The average subordinate function values under the five fertilization gradients all peaked in 2020, with maximum values of 0.33, 0.40, 0.61, 0.72, and 0.74, respectively.

The establishment and management modes of cultivated grasslands have greatly changed the biological and abiotic factors of grassland ecosystems, and modified their energy and material circulation patterns. Fertilization is one of the most direct and effective management modes of cultivated grasslands. For cultivated grasslands, forage height and coverage significantly correlate with forage yield and species. In this study, the height and yield of the grassland of M. sativa were both significantly higher than those of the grassland of D. glomerata. The main reason is that M. sativa is an erect-type top-grass, while D. glomerata is a rhizomatous top-grass, leading to significant differences in forage height and yield between different types of cultivated grasslands. Reasonable fertilization can increase the forage height and coverage of grasslands. In this study, the plant height and coverage of M. sativa both peaked under the fertilization gradient of 20 t/hm², and were 30.14%–57.16% and 172.04%–304.41% higher than those of untreated M. sativa, respectively. The plant height and coverage of D. glomerata gradually increased with increasing fertilization gradient, both peaked under the fertilization gradient of 40 t/hm², and were 19.64%–28.55% and 85.69%–211.21% higher than those of untreated D. glomerata, respectively. These findings are consistent with a previous study (Zhou et al., 2005), which also found that the forage height and coverage of cultivated grasslands increased by 15.40%–83.10% and 10.60%–174.0% after fertilization, respectively. Based on the results from comprehensive evaluation of subordinate function ( Table 4 ), the contribution of coverage to the subordinate function of M. sativa was the greatest under the fertilization gradient of 20 t/hm², suggesting that 20 t/hm² is the optimal organic fertilizer application rate for M. sativa. The subordinate function values of the height and coverage of D. glomerata showed a trend of gradual increase with increasing fertilization gradient. However, under the fertilization gradient of 20–40 t/hm², the difference in the subordinate function value of height was only 0.09%–0.18%, and that in the subordinate function value of coverage was only 19.27%–23.51%. In terms of fertilization costs, 20 t/hm² is likely to be the optimal organic fertilizer application rate for D. glomerata as well.

Fertilization has the most significant effect on promoting the growth and development, yield increase, and income increase of forage. It boosts the tillering of D. glomerata and M. sativa, and increases the dry matter weights of aboveground and underground parts. Fertilizer application rate is positively correlated with forage yield and quality within a certain range of fertilization, but yield begins to decrease beyond the critical application rate of the organic fertilizer. In this study, the critical application rate of the organic fertilizer for M. sativa was 20 t/hm². The critical application rate for D. glomerata was greater than 40 t/hm², under which the yield-increasing effect on D. glomerata was already weaker than that under a low fertilization gradient. Thus, the critical application rate for D. glomerata should be around 40 t/hm². In this study, the critical application rate of the organic fertilizer for M. sativa was lower than that for D. glomerata, mainly because legume plants can coexist with rhizobia and provide a large amount of nitrogen, which reduces the demand of M. sativa for the organic fertilizer. Bai and Qiu et al. identified the nitrogen application rate of 105–300 kg/hm² per year as the threshold of mature or degraded grasslands, and that, beyond this application rate, changes in grassland biomass, species diversity, and functional group structure would no longer be obvious (Qiu and Luo, 2004; Bai et al., 2010). In this study, the critical application rates of the organic fertilizer for M. sativa and D. glomerata could be converted to nitrogen at 446 and 892 kg/hm², which were much higher than the threshold of nitrogen addition. This can be explained from two aspects: First, the cultivated grasslands in this study were newly cultivated land with low soil nutrients ( Table S1 ). Second, cultivated grasslands consume a lot of soil nutrients due to continuous cutting and output, thus only a high application rate could offset the absorption of soil nutrients by forage. The yields of M. sativa and D. glomerata fluctuated significantly between years. The yield of M. sativa was low only in 2019, but remained at a high level over the next three years. In 2019, the first year after the establishment of the grassland of M. sativa, fertilization preferentially promoted root growth, and most nutrients were allocated to roots, thereby accelerating the establishment of fine root meristems and roots (Beck, 1996; Takei et al., 2001). As a result, the growth of aboveground parts lagged behind that of roots. Meanwhile, the yield of D. glomerata was low in 2019 and lowest in 2022, which might be related to the climate anomaly of that year. That is, the summer of 2022 had a higher temperature and less rainfall relative to previous years, and, when the temperature was above 28°C, the growth and tillering of D. glomerata slowed down (Zhang et al., 2019). Moreover, D. glomerata is a shallow-rooted forage, the yield of which responds significantly to temperature rise and drought (Thivierge et al., 2016).

Nutrition quality is the basis for evaluating the quality of forage, and nutrient content is one of the most important indicators to consider while measuring the feeding value of forage. CP and EE contents largely determine the nutritional quality of forage (Hu et al., 2014). In this study, the CP and EE contents of M. sativa and D. glomerata presented a trend of increasing first and decreasing afterward with increasing fertilization gradient. The CP and EE contents after fertilization were significantly higher than those without fertilization. Under the fertilization gradient of 10–30 t/hm², the CP and EE contents of M. sativa peaked, while the NDF and ADF contents reached their minimum values. In terms of forage nutrition, the optimal organic fertilizer application rate for M. sativa should be 10–30 t/hm². D. glomerata had the highest CP and EE contents and the lowest NDF and ADF contents under the fertilization gradient of 20–30 t/hm². In terms of forage nutrition, the optimal organic fertilizer application rate for D. glomerata should be 20–30 t/hm². In practical production, fertilizer application rate depends on the goal pursued. For example, a moderate application rate can be considered if economic benefits are taken as the goal, whereas a high application rate can be implemented if yield is the goal.