Potential of Establishing the Universal Critical Nitrogen Dilution Curve for Japonica Rice

Establishing the universal critical nitrogen (NC) dilution curve can assist in crop N diagnosis at the regional scale. This study conducted 10-year N fertilizer experiments in Yangtze River Reaches to establish universal NC dilution curves for Japonica rice based on simple data-mixing (SDM), random forest algorithm (RFA), and Bayesian hierarchical model (BHM), respectively. Results showed that parameters a and b were affected by the genetic and environmental conditions. Based on RFA, highly related factors of a (plant height, specific leaf area at tillering end, and maximum dry matter weight during vegetative growth period) and b (accumulated growing degree days at tillering end, stem–leaf ratio at tillering end, and maximum leaf area index during vegetative growth period) were successfully applied to establish the universal curve. In addition, representative values (most probable number [MPN]) were selected from posterior distributions obtained by the BHM approach to explore universal parameters a and b. The universal curves established by SDM, RFA, and BHM-MPN were verified to have a strong N diagnostic capacity (N nutrition index validation R2 ≥ 0.81). In summary, compared with the SDM approach, RFA and BHM-MPN can greatly simplify the modeling process (e.g., defining N-limiting or non-N-limiting groups) while maintaining a good accuracy, which are more conducive to the application and promotion at the regional scale.


Introduction
Rice (Oryza sativa L.) is the staple crop of and is widely cultivated in China [1]. Nitrogen (N) fertilizer can promote crop yield and affect the formation of final grain quality [2]. However, the N fertilizer overuse has a negative impact on crop production and environmental security [3]. Precise management of N fertilizer helps to prevent waste or shortage and effectively promotes crop growth and development [4,5].
Morphological comparison, chemical analysis, and instrumental detection are often used for determining the crop N nutrition [6][7][8]. The critical N (N C , %) dilution curve is considered as a reliable indicator in crop N status diagnosis [9,10]. Lemaire et al. [11] firstly proposed the concept of N C dilution curve based on the relationship between plant N concentration and dry matter weight (DW, t ha −1 ). Plant N C is defined as the minimum N concentration when DW reaches the maximum growing rate [12]. Under a certain biomass level, it is optimal for plant N concentration to reach the N C value. Plant N concentration higher than the N C value means that crop growth is not limited by N nutrition, while plant N concentration lower than the N C value represents N-limiting status. In addition, the N nutrition index (NNI) was derived from the N C dilution curve, which can effectively diagnose the N status of crops [13].
It is generally believed that an NNI value of 1 represents the appropriate state of crop N nutrition.
At present, some studies have established independent N C dilution curves under specific cultivars, experiments, or other limitations. Parameters a and b varied between different N C dilution curves due to cultivar and environmental condition differences, which increased the uncertainty [14][15][16]. In order to explore the variation sources of uncertainties, some scholars have focused on the uncertainty analysis of parameters a and b from crops' N C dilution curves, and results showed large differences between various species, cultivars, and field managements [17,18]. However, the traditional method of establishing N C dilution curves by sampling and weighing measurements is time-consuming and laborious, and a new approach for rapid and accurate estimation of parameters a and b is urgently needed. Commonly used linear regression methods are difficult to deal with this kind of multi-factor problem. Meanwhile, machine learning methods are widely used in these issues because they can effectively weaken the collinearity influence between predictors [19,20].
In addition to parameter estimation of specific N C curves, determining a universal N C dilution curve is more practical for regional scale. Some works have been carried out to establish universal N C dilution curves and realize the application for specific areas. Greenwood et al. [21] attempted to construct an N C dilution curve for C3 plants, while limited modeling data could not accurately represent all the C3 crops. As for rice, relationships between plant N concentration and aboveground biomass varies under different cultivars, especially Japonica and Indica [22]. Therefore, it is worth exploring to establish the universal N C dilution curve by setting partial restriction conditions (e.g., similar climate, cultivar, and cultivation region). For the regional scale, the universal N C dilution curve is usually represented by simple data-mixing (SDM) of fitting all the N C points determined by defining N-limiting and non-N-limiting groups and the linear-plateau model. Although this traditional approach has a good applicability, the process of determining all the N C points is time-consuming. Moreover, many researchers adopted different growth parameters (e.g., leaf area index [LAI] and leaf dry matter) instead of the traditional plant DW, and the method of averaging the parameters from specific N C dilution curves was verified in application and mathematical statistics [23][24][25]. These provide the possibility for calculating universal curve parameters according to estimated specific curve parameters. In addition to traditional simple linear regressions (e.g., multiple linear regression [MLR]), random forest algorithm (RFA) should also be considered for its high accuracy and wide application in crop parameter estimation. As a learning method integrating multiple decision trees, RFA can process large-scale information and reduce noise [26]. As another way, Makowski et al. [17] also proposed to adopt a Bayesian hierarchical model (BHM) approach to skip the step of distinguishing N-limiting or non-N-limiting groups. Combined with the prior knowledge, the posterior distribution of N C dilution curve parameters can be obtained through the Markov Chain Monte Carlo (MCMC) method. As the main rice-producing area in China, Yangtze River Reaches is planted with plenty of Japonica rice cultivars [27]. Establishing the universal N C dilution curve for Japonica rice will be of great significance to diagnose crops' N status at the regional scale.
In this study, specific N C dilution curves were established based on plant DW and N concentration in the vegetative growth period from numerous N fertilizer experimental data of Japonica rice in Yangtze River Reaches. Overall, the objectives were (a) to analyze the sources and uncertainties leading to differences in each specific N C dilution curve; (b) to explore the potential of estimating parameters from specific N C dilution curves using MLR and a machine learning method (RFA); and (3) to establish a universal N C dilution curve for Japonica rice in Yangtze River Reaches based on SDM, RFA, and BHM approaches. These results will provide theoretical and technical support for field N precise management of Japonica rice.  7), Wuyunjing 24 (WYJ 24), Nanjing 9108 (NJ 9108), Yongyou 2640 (YY 2640), and Wuyunjing 32 (WYJ 32). In addition, experiments were also divided into 5 different crop management groups according to N application rates, fertilization ratio arrangement, and sampling frequency (A group: Exp. 1 and 2; B group: Exp. 3, 6, and 7; C group: Exp. 4 and 5; D group: Exp. 8 and 9; and E group: Exp. 10, 11, 12, and 13). Other specific information is presented in Table 1.

Data acquisition
At least 15 fresh rice plants were selected at critical growth stages from each plot in 13 experiments and separated by organs (leaf, stem, and spike). A portable li-3000c leaf area meter (Li-Cor., Lincoln, USA) was used for obtaining leaf area, and the LAI was calculated in Eq. 1 [28]. All samples were first heated for 30 min at 105 °C in a forced-draft oven, and then dried at 80 °C to a constant weight. Finally, dry matter weight of samples was weighed by a micrometer balance.
Plant DW includes DW of leaves, stems, and spikes, respectively. The specific formula of organ DW is as follows: where r (cm 2 ) represents the readings of the li-3000c leaf area meter, w (g) stands for the dry matter weight of samples (leaf, stem, or spike), s represents tiller number in 1 m 2 from each plot, and k denotes the tiller number of sampled rice plants.
Stem-leaf ratio (SLR) represents the distribution of plant dry matter; it can be calculated as follows: in which SDW (g) and LDW (g) are the dry matter weight of stems and leaves, respectively.
The dry matter of each organ was ground into powder by a grinder and passed through a 1-mm sieve. Each sample was weighed to 0.15 g and N concentration (N%DW) was determined in a continuous-flow auto-analyzer (BRAN + LUEBEE, Hamburg, Germany) after Kjeldahl digestion and an isochoric process.
Accumulated growing degree days (AGDD, °C) is the cumulative value of growing degree days (GDD, °C). GDD refers to the cumulative effective accumulated temperature experienced by completing a certain growth stage under actual environmental conditions. They are calculated in Eqs. 4 and 5: in which T MAX (°C) and T MIN (°C) are the maximum and minimum temperature in a day, respectively, T B (°C) is the base point temperature of crop development (12.5 °C for Japonica rice) [29], and k is the number of days after rice transplanting.
Rice plant DW reaches 1 t ha −1 at about the end of tillering, and the plant N concentration at this time represents the value of parameter a. Therefore, the end of tillering was closely related to parameters a and b and was set as a time point to construct the indicators. A total of 14 factors were selected for parameter estimation, where thousand grain weight (TGW, g), amylose content (AC, %), plant height (PH, cm), specific leaf area at tillering end (SLA-T), and SLA at vegetative growth end (SLA-V) stood for the cultivar characteristics of Japonica rice; AGDD over the whole growth period (AGDD), AGDD at tillering end (AGDD-T), and AGDD during the vegetative growth period (AGDD-V) reflected the differences of growth days and temperature caused by years and regions; and maximum DW at the tillering stage (DW-T), maximum LAI during the tillering stage (LAI-T), SLR at tillering end (SLR-T), maximum DW during the vegetative growth period (DW-V), maximum LAI during the vegetative growth period (LAI-V), and SLR at vegetative growth end (SLR-V) indicated the impact of various influencing factors (cultivar, region, year, and crop management) on rice growth. TGW, AC, and PH were reference values  of rice cultivar parameters obtained from the China Rice Data Center (https://www.ricedata.cn). They were used to represent the inherent potential of the rice cultivar, and do not mean that the measured value will not change due to cultivation environment. However, there were still great uncertainties in the application of these empirical indicators. SLA is the ratio of one-side leaf area to the corresponding DW [30], which is also an indicator of cultivar parameters in many crop growth models [31,32]. AGDD was calculated from the meteorological information of the local station. In addition, DW, LAI, and SLR were typical growth indicators measured through experimental tests [4,33].

Approaches for establishing the universal N C dilution curve
The following 3 approaches were selected to establish the universal N C dilution curve in this study: a. The SDM approach (i) defined non-N-limiting and N-limiting groups through variance analysis, (ii) determined the N C point based on the linear-plateau model, and (iii) fitted all the N C points to a power function curve.
b. The RFA approach (i) predicted parameters a and b of specific N C dilution curves, and (ii) calculated parameters a and b of the universal N C dilution curve by averaging the estimated parameters of specific curves. RFA takes the decision tree as the basic model and generates a series of decision tree models by constructing different training datasets and feature spaces. It is suitable for solving complex nonlinear problems [34]. NTREE (number of decision trees) and MTRY (number of observations per tree leaf) are critical tuning parameters in RFA. NTREE affected the operation time and was set to 500 as commonly used in previous studies [5]. The MTRY value had a great influence on the modeling accuracy, which was adjusted to achieve the smallest relative root mean square error (RRMSE) according to the input variables (Table 3). In addition, the RFA approach was compared with the MLR method, which is often used for dealing with multi-factor problems [35]. The 2 methods were both used for parameter estimation of specific N C dilution curves, and validated by the leave-one-out cross-validation method. This validation method leaves only one sample as the test set and other samples as the training set. If there are k samples, k times training and testing are needed. It is suitable for a small number of samples and has high sample utilization [36].
c. The BHM approach (i) obtained the posterior distributions of parameters a and b based on the BHM and MCMC, and (ii) selected representative values (Mean, Median, and most probable number [MPN]) from the posterior distributions as universal curve parameters. The BHM was constructed based on the linear-plus-plateau of biomass to N concentration, which also considered the factors of measuring the date's variation and prior knowledge, and the detailed process can be found in Makowski et al. [17]. MCMC, a parameter uncertainty analysis method based on Bayesian statistical theory, was adopted for estimating posterior distributions of universal N C dilution curve parameters. The principle is to construct multiple stable Markov chains according to the parameter prior knowledge and model observation information, and explore the probability distribution of the parameters [24]. This method can be implemented by calling the rjags package in R Programming Language 3.6.0 software (R Foundation for Statistical Computing, Vienna, Austria). In this study, prior 1 (weakly informative priors) was refined without strong limitations (a ranged from 0 to 6, and b ranged from 0 to 1) for better estimation. The algorithm converged after about 50,000 iterations, and then additional 50,000 iterations were performed. Finally, posterior distributions were described by a frequency distribution histogram and kernel density function [37].

Data processing and analysis
N C dilution curve indicates that plant N concentration decreases with the increase of DW, which conforms to the form of power function: in which a stands for plant N concentration when the biomass is 1 t ha −1 , and b is the parameter affecting the slope of the curve. N concentration above the N C dilution curve indicates excess N, N concentration on the curve means suitable N, and N concentration below the curve represents limited N.
Standard deviation (SD) and coefficient of variation (CV) were adopted to describe the basic attributes of the dataset, and the larger values of the two, the more scattered the data. They are calculated in Eqs. 7 and 8: where m is the average value of total samples, and n stands for the number of samples.
Fluctuation range was introduced to represent the variation degree of parameters under different environmental conditions; the specific formula is as follows: in which fluctuation range was calculated as the ratio of distribution range to mean value. R 2 and RRMSE were used to represent the accuracy and stability of models, respectively. The closer R 2 was to 1 and the smaller RRMSE was, the better the model was. They are described in Eqs. 10 and 11: in which m and n stand for the estimated values and measured values, respectively; m and n stand for the average estimated values and measured values, respectively; and k stands for the number of samples.
NNI is the ratio of measured N concentration to N C , which was used in this study to verify the N diagnostic capacity of the established universal N C dilution curves. It can be calculated as: where NC T (%) represents the real-time measured N concentration and N C (%) is calculated by the N C dilution curve. An NNI value greater than 1 indicates N surplus, an NNI value equal to 1 indicates N suitability, and an NNI value less than 1 indicates N deficiency. Statistical analysis was performed using IBM SPSS 25 software (IBM Corporation, Armonk, USA). Significance analysis was completed in R Programming Language 3.6.1 (R Foundation for Statistical Computing, Vienna, Austria) by comparing the distributions of parameters.

Technical flowchart
In this study, we established specific N C dilution curves based on 10 years' experiments in Yangtze River Reaches. Uncertainty factors affecting parameters a and b were analyzed from cultivar, year, region, and crop management. SDM, RFA, and BHM approaches were adopted to explore the establishment of a universal N C dilution curve for Japonica rice. The whole technical flowchart is shown in Fig. 1.

Establishment of specific N C dilution curves
As shown in Fig. 2, a total of 12 specific N C dilution curves were established from field experiments. In general, parameter a ranged from 3.3958 to 3.6275, and parameter b ranged from 0.2390 to 0.4580. There are large differences in the confidence interval widths of curves established by different experiments. The average confidence interval widths of Exp. 1 ( Fig.  2A), Exp. 2 (Fig. 2B), and Exp. 6 ( Fig. 2E) are 0.56, 0.62, and 0.65, respectively, which are much higher than other experiments (0.25 to 0.48). This is closely related to the sampling times, the fewer sampling times, the greater systematic errors, and the larger confidence width. The confidence width of all curves decreased with the increase in biomass and tended to be stable after the biomass reached 2.5 t ha −1 . In addition, when plant DW of rice reached 1 t ha −1 , Japonica rice plants showed similar N concentration due to weak N dilution function at the early growth stage (parameter a ranged from 3.40 to 3.63). Furthermore, confidence interval widths under this aboveground biomass are in the range of 0.45 to 1.37, which indicated a high reliability.

Sources of variation in N C dilution curves
As shown in Fig. 3, a total of 29 specific N C dilution curves were established based on the BHM and MCMC methods. Parameters a and b varied with different genetic and environmental conditions (corresponding to cultivar, year, region, and crop management in actual cultivation), resulting in differences between N C dilution curves. Different cultivars (e.g., WXJ14YZ2010A and LXY18YZ2010A, WXJ14YZ2011A and LXY18YZ2011A, and WYJ19WJ2013B and YY8WJ2013B), years (e.g., WXJ14YZ2010A and WXJ14YZ2011A, LXY18YZ2010A and LXY18YZ2011A, and NJ9108XH2017E and NJ9108XH2018E), and region × crop managements (e.g., WYJ24SH2015D and WYJ24HA2015B, NJ4SH2015D and NJ4HA2015B, and LJ7SH2015D and LJ7HA2015B) all had important influences on parameters a and b. The median values of the posterior distribution were adopted to represent the trend of the specific N C dilution curves, and the descriptive statistics of parameters a and b are shown in Table 2. The values of parameter a ranged from 3.26 to 3.98, while those of parameter b ranged from 0.20 to 0.50. The CV for parameter b reached 23.82%, indicating that values fluctuate greatly due to various environmental conditions.
Considering the actual cultivation environment, the impact of genetic and environmental conditions on parameters a and b was analyzed from cultivar, region, year, and crop management, respectively (Fig. 4). There were significant differences (p < 0.05) in parameters a and b among different cultivars (except parameter a in NJ 9108 and WXJ 14, Fig. 4A-1 and A-2), regions (Fig. 4B-1 and B-2), years (Fig. 4C-1 and C-2), and crop managements (Fig. 4D-1 and D-2). On the whole, the fluctuation range of parameter a was 53.78% smaller than that of parameter b on average. Parameter b was more strongly influenced by the comprehensive effects of genetic and environmental conditions. Parameter a changed synchronously with parameter b, and the large fluctuation range of parameter a was accompanied by the large fluctuation range of parameter b (e.g., cultivar-WYJ 19, region-ZJG, year-2017, and crop management-C). For parameter a, the distribution range of median values calculated from the posterior distribution between different cultivars, regions, years, and crop managements was 3.06 to 3.87, 3.40 to 3.67, 3.42 to 3.74, and 3.36 to 3.64, respectively. The performance of parameter a was relatively stable among different regions, years, and crop managements, while the cultivar had a great influence on parameter a. In addition, the median values of parameter b were at a range of 0.30 to 0.43, 0.24 to 0.46, 0.26 to 0.45, and 0.26 to 0.47 on different cultivars, regions, years, and crop managements, respectively. Although the median fluctuation range of parameter b was relatively smaller among different cultivars, it still fluctuated sharply (Fig. 4A-2).

Parameter estimation of specific N C dilution curves using RFA and BHM
The N C dilution curve is influenced by genetic and environmental conditions, which include cultivar, region, year, and crop management in actual cultivation. Therefore, influencing factors from cultivar (TGW, AC, PH, SLA-T, and SLA-V), year, region, and crop management (AGDD, AGDD-T, AGDD-V, DW-T, LAI-T, SLR-T, DW-V, LAI-V, and SLR-V) were selected respectively for importance analysis (Fig. 5). Coefficient of determination (R 2 ) varied with different factors and parameters, where the 3 factors most related to parameter a were PH, SLA-T, and DW-V, and the 3 factors of largest correlation with parameter b were AGDD-T, SLR-T, and LAI-V. The correlation between TGW and parameters a and b was poor, while SLA-T, SLR-T, and LAI-V achieved good performance. In addition, PH, SLA-T, and DW-V predicted parameter a better than b, while AGDD-T, LAI-V, and SLA-V were more suitable for paramter b estimation.
The 3 highly related factors of parameters a (PH, SLA-T, and DW-V) and b (AGDD-T, SLR-T, and LAI-V) were selected as predictors for estimation, respectively, and leave-one-out cross-validation was adopted for subsequent NNI verification. The RFA method had more advantages than the MLR method in predicting parameters a and b (Table 3). For parameter a, the RFA method was higher than the MLR method 6.58% and 3.64% in modeling R 2 and validation R 2 , and lower than 3.46% in RRMSE value. In addition, the RFA method was also superior to the MLR method in modeling R 2 , verification R 2 , and RRMSE of predicting parameter b by 25.78%, 11.30%, and 4.28%, respectively. Figure 6 shows parameter posterior distributions through the BHM and MCMC methods, and the results were presented   Table 4).

The NNI validations of the established universal N C dilution curves
The NNI value calculated by each specific N C dilution curve based on the traditional method was taken as the reference value (measured value). As shown in Fig. 7, the NNI verified R 2 of MPN reached 0.82, while Mean and Median had a poor ability of N diagnosis. There were 2 other approaches (RFA and SDM) for establishing the universal N C dilution curve (Table 4), and the corresponding validation results of NNI are shown in Fig. 7. The RFA approach aimed to determine parameters a and b by averaging the estimated parameter values from the specific curves above (Table 4, RFA, a = 3.56, b = 0.34) with an NNI validation R 2 of 0.81 (Fig. 7C). In addition, the average value of specific N C dilution curves also showed a good validation result of NNI (Table 4, Average, a = 3.56, b = 0.34, Fig. 7B, The SDM approach aimed to fit all the N C points and build the mixed curve (Table 4, SDM, a = 3.36, b = 0.30), which obtained a high NNI validation R 2 of 0.8138 (Fig. 7A). In general, there was no significant statistical difference of NNI validations between the universal N C dilution curves based on SDM, Average, RFA, and MPN. However, Mean and Median were far less capable of N diagnosis than MPN. Results indicated that universal N C dilution curves established by SDM, Average, RFA, and MPN had strong applicability and robustness.

Uncertainty analysis of parameters a and b
Fewer data points and more outliers resulted in higher confidence interval and uncertainty of specific N C dilution curve in each experiment, which is due to the lower sampling frequency (e.g., Exp. 1, Exp. 2, and Exp. 6). This can also well explain the N C dilution curve in Exp. 6, which had a large confidence width because of no sampling points before plant DW reached 3 t ha −1 .
A smaller sample number also resulted in a larger posterior distribution of parameters a and b based on the BHM approach (e.g., cultivar-WYJ 19, region-ZJG, and year-2017). At about the tillering end of rice, the plant DW reaches 1 t ha −1 , and the N concentration at this time represents the value of parameter a. This is an important stage closely related to the value of parameter a, which is also the reason for using this stage to construct parameter predictors (e.g., AGDD-T, DW-T, LAI-T, SLR-T, and SLA-T). Parameter a can reflect the N absorption characteristics of rice plants, which is mainly determined by cultivars. This study aimed at different Japonica rice cultivars; thus, the fluctuation range of parameter a was relatively stable even under different environmental conditions. Compared with region, year, and crop management, parameters a and b both had a smaller fluctuation range among different cultivars. There were many factors affecting parameter b a lot, which made it difficult to predict it accurately. AGDD-T, SLR-T, and LAI-V were the factors highly related to parameter b. Plant N concentration dilution is caused by the physiological process of crop growth and development. Plants grow higher to obtain more sunlight, leading to the change of dry matter distribution, which is greatly associated with parameter b [38] and reflects on SLR. Moreover, the shading phenomenon between leaves becomes more obvious with the progress of crop growth.
In order to make effective use of N nutrition in crop canopy,  N in shaded leaves was transported to the top leaves, resulting in the decrease of leaf N concentration. This process is closely related to LAI, which is an important indicator of parameter b. AGDD-T indicated that parameter b was affected by growth days and temperature. As another characterization of crop growth status other than LAI, DW (DW-V) also had a great impact on parameter a. Nevertheless, rice cultivar characteristics affected parameter a a lot, and PH and SLA-T are both such indicators. Moreover, DW-V represents the change in crop biomass influenced by various factors. There was no temperaturerelated factor that showed great impact on parameter a, indicating that the environmental conditions' experimental sites were similar and relatively stable at the early growth stage of Japonica rice.

Parameter estimation of specific N C curves
This study provided a new approach for predicting parameters a and b from specific N C dilution curves, while the unity of predictors still needs to be verified in the future due to the limited experimental data. In particular, PH showed a very high contribution to parameter a estimation, which may be caused by a limited dataset. However, the reference values (TGW, AC, and PH) from the China Rice Data Center were empirical with strong uncertainty. Although only PH was used in the prediction of parameter a in this study, more reliable cultivar indicators need to be explored. For the selection of indices, there was still a lot of room to explore (e.g., soil indicators and management measurements). Dry matter accumulation and N concentration changes of japonica rice involve complex physiological processes, and simple linear regression models are difficult to deal with such nonlinear problems. Nevertheless, this issue can be better handled by the machine learning method (RFA), which can be compatible with input samples from various dimensions, and also achieve accurate prediction in the case of lacking value [39]. RFA has been widely used in crop yield prediction and other aspects [40,41]. In order to enhance the applicability of parameter estimation in specific N C dilution curves, only highly related factors were selected as predictors. In fact, factors affecting the value of parameter a and b are different, and relatively stable factors with higher importance are expected to be found in the future. In addition, the parameter tuning process can be carried out on the RFA method in this study, which may further improve the prediction accuracy [42]. Furthermore, more machine learning methods (e.g., artificial neural network and support vector machine) should also be taken into consideration to compare and explore better solutions [43,44].

The establishment of universal N C dilution curves
It has always been a hot issue to establish a reliable universal N C dilution curve with strong applicability [45]. This study explored 3 feasible approaches for Japonica rice in Yangtze River Reaches. SDM, as a traditional approach, achieves a high verification accuracy of NNI through fitting all the N C points. However, the determination of N C point and subsequent data processing are undoubtedly time-consuming [46], which is not conducive to the actual production application. Therefore, simplifying the complex steps while ensuring N diagnostic capacity is of great significance. In this study, results showed that universal N C dilution curves determined by averaging parameters a and b from specific curves also had good N diagnostic ability, which is consistent with Yao et al. [24]. Moreover, on the basis of predicting specific N C dilution curve parameters with highly related indicators, it can be extended to the rapid estimation of universal curves by averaging parameters a and b. This approach can greatly improve modeling efficiency, but more practice is needed in the selection of predictors. In addition, the BHM for N C dilution curves proposed by Makowski et al. [17] can also skip the step of distinguishing N-limiting and non-N-limiting groups. Parameter posterior distributions can be easily obtained from basic data and prior knowledge through BHM and MCMC methods. Previous research always adopted the MCMC method for uncertainty analysis of model parameters, while few focused on parameter estimation and verification [47,48]. According to the posterior distribution of parameters a and b, 3 representative values (Mean, Median, and MPN) were selected to explore the potential of establishing universal N C dilution curves. Although MPN achieved excellent results of NNI validation, there were still some uncertainties. On one hand, prior knowledge almost directly determined the upper and lower limits of parameter posterior distributions. On the other hand, selected representative values may not fully represent the actual values, even if they have been applied with high precision. The applicability of MPN was proved in this study, and its stability is expected to be verified in the future. In this study, the sources of the differences among specific N C dilution curves were from cultivars, years, regions, and crop managements. For Japonica rice, parameter a is mainly determined by the N absorption characteristics of cultivars, while parameter b is related to cultivars and environmental conditions (e.g., AGDD and rain-heat resources). Highly related indicators (PH, SLA-T, DW-V, AGDD-T, SLR-T, and LAI-V) were selected from 14 factors to explore the potential of predicting parameters a and b. In the parameter estimation of specific N C dilution curves, the RFA method was, on average, 14.69% and 6.98% higher than the MLR method in the accuracy of modeling and verification, respectively. In addition, SDM, Average, RFA, Mean, Median, and MPN determined parameters of the universal N C dilution curves for Japonica Rice in Yangtze River Reaches, among which, SDM, Average, RFA, and MPN obtained a good performance with NNI validation R 2 above 0.80. The good validation accuracy of Average provided theoretical support for the implementation of the RFA approach. Compared with the SDM approach, RFA and MPN approaches simplified complex steps, improved operational efficiency, and showed strong applicability.
Funding acquisition, Resources, Project administration, and Supervision. Competing interests: The authors declare that they have no competing interests.

Data Availability
The data used to support the findings of this study are included within the article.