## Classification of production systems and construction of database

We categorized China’s wheat systems into seven ecological zones (i.e., NE, R-NW, I-NW, I-NCP, R-NCP, YR, and SW) covering 2096 counties. This classification was determined by considering factors such as soil properties, weather patterns/characteristics, agricultural practices, and geographical conditions (Figure S2). The information regarding sample number, N application rate, basic soil properties, and weather patterns of seven wheat production zones is given in Table S2.

We have searched in China Knowledge Resource Integrated Database (http://www.cnki.net/) and the Web of Science (https://www.webofscience.com/) or literature published between 1980 to 2020 to collect data on the effects of N fertilizer on wheat yield in China. The search terms employed “wheat”, “nitrogen”, and “yield”. Studies meeting the following criteria were selected: 1) clear documentation of the experimental time and locations; 2) involving wheat field trials conducted in China (excluding potted plants, greenhouses, etc.); 3) documenting yearly yield data (instead of multi-year average yields); 4) having three or more N rate treatments; 5) demonstrating a yield response reaching a yield plateau; 6) utilizing urea or ammonium as N application sources, excluding enhanced-efficiency fertilizers; and 7) excluding trials involving irrigation treatments. A total of 369 papers, comprising 3726 observations, were selected for analysis. Data presented in the form of figures was extracted using GetData Graph Digitizer (GetData Pty Ltd, Kogarah NSW 2210, Australia). The data obtained from the selected literature includes the title, the author, trial time and locations, wheat yield, N rates, and soil properties (i.e., pH, SOM, TN). Figure S11 illustrates the geographical distribution of the experimental sites. Various biophysical and socioeconomic factors, including pH, SOM, TN, GAT, GR, GS, and FI, were collected. More detailed information can be accessed in the Supplementary methods.

## Regional optimal nitrogen rates

The division of government departments can hinder the progress of the SDGs as they tend to prioritize their mandates. For instance, the Ministry of Agriculture and Rural Affairs focuses on matters such as fertilizer input, farmers’ income, and crop yield, while environmental issues like Nr losses are the domain of the Ministry of Ecology and Environment21. Considering this, three scenarios with specific targets in mind were conceptualized to assess how coordinated efforts by different government departments contribute to sustainability. To achieve this, the study developed three scenarios for optimizing N usage.

Scenario I: S1, focused on optimizing a single sustainable objective (maximizing either crop yield or economic benefits, or minimizing environmental pollution). YON was determined by identifying the N application rate that either crop maximum relative yield, allowing yield data to be compared with each other across growing seasons and locations. ECON was established to determine the N application rates that led to maximum normalized economic benefits using the range normalization method to calculate net income as shown in Eq. (1). ENON was defined as the N application rate that aimed to classify the N application rates that achieved a balance (apparent N balance) between N input and output. Apparent N balance46, as an environmental index, was used to evaluate the difference between fertilizer input and crop removal. This approach is simpler and more cost-effective46 as systematic nutrient balance that considers all input-output pathways. The calculation process followed Eq. (2):

$$\text{Economic benefits}\text{ }\text{=}\text{ }\left({\text{W}}_{\text{y}}\text{ }\text{×}{\text{ }\text{W}}_{\text{p}}\right) \text{-}\text{ }\left({\text{N}}_{\text{r}}\text{ }\text{×}\text{ }{\text{N}}_{\text{p}}\right)$$

1

where Wy, Wp, Nr, and Np represent wheat yield, wheat price, N application rate, and N fertilizer price, respectively. The price data for different years comes from the China Statistical Yearbook (http://www.stats.gov.cn).

$$\text{N balance}\text{ }\text{= }{\text{N}}_{\text{input}} \text{-}{\text{ }\text{N}}_{\text{har}}$$

2

where N balance, Ninput, and Nhar represent apparent N balance, N fertilizer input, and crop N harvest. Nhar was calculated from the product of N requirements per megagram grain and grain yield.

Most commonly used parametric models, such as linear, quadratic, linear plateau, and quadratic plus plateau functions, were employed to estimate crop yield responses. However, the inherent uncertainty in the model structure can introduce significant noise into the results. To address this issue, a combination of parametric and nonparametric models without specific assumptions about the data distribution characteristics was used to estimate the optimal N rates in this study. The choice of the models was based on their definitions, as outlined in Table S3. The methodology of weighing different models based on the inverse of the root mean square error (RMSE) was in line with approaches employed in previous studies47. This approach assisted in determining the best YON, ECON, and ENON values, while also enhancing the interpretability of the model.

Scenario II: S2, involved optimizing multiple sustainable objectives by assessing regional ON thresholds based on the interval ranges of YON, ECON, and ENON (Figure S12).

Scenario III: S3, the focus was on simultaneously optimizing multiple sustainable objectives, particularly coordinating yield, economy, and environment. SN was calculated as the N rates at the maximum NEEB (Eq. (3)). Notably, the SN is needed to fall within the established optimal N thresholds. If exceeded this boundary, it was constrained to the closest value (Figure S12).

$$\text{N}\text{E}\text{E}\text{B}=\left({\text{W}}_{\text{y}}\text{×}{\text{W}}_{\text{p}}\right)\text{-}\left({\text{N}}_{\text{r}}\text{×}{\text{N}}_{\text{p}}\right) \text{-}\text{ }\text{E}\text{cost}$$

3

where Ecost represented the environmental costs including the damage costs of water eutrophication, GHG emissions, soil acidification, and human health cost. The detailed information is listed in Table S4.

A bottom-up approach was employed to estimate the total wheat production, considering different strategies. The estimation was conducted based on the wheat planting area in 2018 (http://www.cnki.net/).

## Calculation of sustainability indicators

The sustainability of a multi-objective N optimization was evaluated in this study, considering environmental, food, and economic indicators. The results were compared with FN obtained from Li et al.48 to evaluate their effectiveness.

$$\text{N}\text{r} \text{l}\text{o}\text{s}\text{s}\text{e}\text{s}\text{=}{\text{ Nr-NH}\text{3}\text{+Nr-N}\text{2}\text{O+Nr-leaching and runoff}}_{ }$$

5

$$\text{N} \text{f}\text{o}\text{o}\text{t}\text{p}\text{r}\text{i}\text{n}\text{t}\text{= Nr losses / Yield}{ }_{ }$$

6

$$\text{NUE}\text{= }\text{N}\text{har}\text{ /(N}\text{fer}\text{ + N}\text{dep}\text{ +N}\text{fix}\text{ ) × 100\%}{\text{ }}_{\text{ }}$$

7

$$\text{S}\text{N}\text{M}\text{I}\text{=}\text{ }\sqrt{{(1-\text{N}\text{Y}\text{i}\text{e}\text{l}\text{d}\text{*})}^{2}+{(1-\text{N}\text{U}\text{E}\text{*})}^{2}}{ }_{ }$$

8

$$\text{N}\text{Y}\text{i}\text{e}\text{l}\text{d}\text{*}=\left\{\begin{array}{c} NYield/NYieldref (NYield\le {\text{N}\text{Y}\text{i}\text{e}\text{l}\text{d}}_{\text{r}\text{e}\text{f}})\\ 1 (NYield>{\text{N}\text{Y}\text{i}\text{e}\text{l}\text{d}}_{\text{r}\text{e}\text{f}})\end{array}\right.$$

9

$$\text{N}\text{U}\text{E}\text{*}=\left\{\begin{array}{c} NUE/NUEref (NUE\le {\text{N}\text{UE}}_{\text{r}\text{e}\text{f}})\\ 1 (NUE>{\text{N}\text{UE}}_{\text{r}\text{e}\text{f}})\end{array}\right.$$

10

where Nfer, Nfix, and Ndep are N fertilizer, biological N fixation, and atmospheric deposition, respectively. NYield*, NUE*, NYieldref, and NUEref represent the normalized yield, normalized NUE, reference Yield48, and reference NUE27 (NUEref = 0.8).

## Meta-analysis

To investigate the impacts of weather patterns, soil properties, and farmers’ income on optimal N application, a random-effects model was employed to analyze environmental variables. The effect size of N treatment was quantified as the natural logarithm of the response ratio (RR) following the approach outlined by Hedges et al.49:

$$\text{ln}\text{RR= }\text{ln}\left(\frac{{\text{X}}_{\text{T}}}{{\text{X}}_{\text{C}}}\right)$$

11

where XT and XC represent the means of N optimization treatments and the non-fertilized treatments for weather, soil, and farmer income indicators. The variance (V) of lnRR (V(lnRR)) and 95% confidence interval was calculated as:

$$\text{V}\text{(lnRR) }\text{= }\frac{{\text{SD}}_{\text{T}}^{\text{2}}}{{\text{N}}_{\text{T}}{\text{X}}_{\text{T}}^{\text{2}}}\text{+}\frac{{\text{SD}}_{\text{C}}^{\text{2}}}{{\text{N}}_{\text{C}}{\text{X}}_{\text{C}}^{\text{2}}}$$

12

$$\text{95% confidence interval =}\text{ln}\text{RR}\text{±1.96}\sqrt{\text{V}(\text{ln}\text{R}\text{R})}$$

13

where \({\text{S}\text{D}}_{\text{T}}\) and \({\text{S}\text{D}}_{\text{C}}\) represent the standard deviations of optimal N application and non-fertilized treatments, respectively. \({\text{N}}_{\text{T}}\) and \({\text{N}}_{\text{C}}\) represent the sample sizes for optimal N application and non-fertilized treatments. The weighted effect sizes (\({\text{l}\text{n}\text{R}\text{R}}_{++}\)) were calculated as:

$$\text{ln}{\text{RR}}_{\text{++}}\text{=}\frac{\sum {\text{w}}_{\text{i}}\text{×}\text{ln}\text{R}{\text{R}}_{\text{i}}}{\sum {\text{w}}_{\text{i}}}$$

14

$${\text{w}}_{\text{i}}\text{ }\text{=}\frac{\text{1}}{{\text{V}}_{\text{i}}}$$

15

where \({\text{w}}_{\text{i}}\) and \({\text{V}}_{\text{i}}\) are the weighted factors of effect size and variance in the ith study, respectively.

Differences between optimal N application and non-fertilized treatments were determined by assessing whether the 95% confidence interval overlapped with 0 (with a significance level of P < 0.05)49. The percent change was calculated to quickly capture the response of the optimal N application as:

$$\text{The change of yield = }\left(\text{R}\text{R}\text{-1}\right)\text{×}\text{ }\text{100}$$

16

## Nitrogen management practices

Nitrogen use efficiency, a sustainability indicator, is influenced by soil quality, N fertilizer input, and weather conditions. To enhance NUE, this study paid special attention to soil quality and N fertilizer inputs as modifiable factors. The trials conducted under the N optimization strategy were classified into different management levels: low management level (yield below the yield corresponding to the average optimizing N application) and high management level (yield above the yield corresponding to the average optimizing N application). A crop-based estimation method was used to assess soil quality, taking into account the mean yield, variation, and frequency distribution of the no N-rate control. Soil quality was classified as high, moderate, or low based on this method (Table S1 and Figure S13).

The EU28 evaluation system provides criteria for determining best management practices. These criteria include N output ≥ 80 kg N ha− 1, NUE falling within the range of 50–90%, and Nsur ≤ 80 kg N ha− 1. There is a significant opportunity for improvement if these criteria are not met. In cases where NUE ≥ 50% but Nsur > 80 kg N ha− 1, priority is given to the use of slow-release fertilizers to minimize residual N released in the environment. Figure S9 presented other recommended methods for improving N management practices.