Fifteen site-year winter wheat field experiments were conducted in agro-ecological zone 1 (AZ1, 37°43'N and 117°13'E, in Laoling County of Shandong Province) and agro-ecological zone 2 (AZ2, 33°05'N and 119°53'E, in Xinghua City of Jiangsu Province) (Figure 1A), with different patterns of precipitation distribution and monthly average temperature (Figure 1B). The soil types of AZ1 and AZ2 were sandy loam soil and loam soil, respectively.

The design includes seven plot experiments (Exp. 1–7) and eight farmers' field experiments (Exp. 8–15) involving different N rates, sowing dates (SD), and seeding rates (SR). The N fertilizers of Exp. 1–3 were applied with 40% as basal N before sowing, except for farmer's management (FM), which applied 50% as basal N. The remaining N was applied at the stem elongation (SE) stage. The N fertilizers of Exp. 4–5 were applied at the rate of 50% of total N rates before sowing and 50% at the SE stage, respectively. The differences between Exp. 6–7 and Exp. 4–5 are listed in Table 1. With a randomized complete block design, the treatments of plot experiments were carried out with three repetitions. In addition to plot experiments, eight farmers' field experiments (Exp. 8-15) were conducted across the Nanxia village in Shandong Province to compare different N management strategies. In each farmers' field experiment, there were three treatments: FM, regional optimum N management (RONM), and PNM. RONM, integrating high-yielding and high-NUE, was managed with 81 kg N ha⁻¹ as basal N and 138 kg N ha⁻¹ as topdressing N (Zhou et al., 2017). For PNM, the basal N rate was the same as RONM, and topdressing N was applied according to an active canopy sensors-based algorithm (Cao et al., 2017b). Other managements were the same as RONM. Irrigation was done two times at the sowing and SE stages in AZ1 and once after sowing in AZ2. All N source was urea. For all plots, wheat was not limited by phosphate and potash fertilizers and was kept free of pests and diseases during the growing seasons.

The data were analyzed using the following steps: (1) establishing the SLR models only using NDRE; (2) identifying optimal variables based on variable selection strategies; (3) training and validating the RFR models using selected variables by updating hyperparameters and cross-validation; (4) comparing and evaluating the robustness of different models; (5) determining the variable contributions in different AZs, and (6) exploring the effective combinations of multi-source data for N status estimation (Figure 2).

A total of 396 observations were obtained, with 277 for training and 119 for the test (Table 2). PNU showed the highest variation on training dataset (CV = 60.53, 51.04 and 52.29% in AZ1, AZ2 and AZ1+2, respectively), followed by AGB (CV = 50.2, 38.53, and 41.56% in AZ1, AZ2 and AZ1+2, respectively) and NNI (33.67, 34.75, and 32.35% in AZ1, AZ2, and AZ1+2, respectively). Similar results were presented in the test dataset. In general, experimental data is suitable to evaluate the predictive performance for N status estimation.

The performance of SLR models for estimating N status indicators varied slightly across the different AZs (Figure 3). Comparatively, NDRE explained 59% of AGB variation, 62% of PNU, and 57% of NNI in AZ1, which performed better than AZ2 with 36, 52, and 47%, respectively. The performance of models in AZ1+2 was similar to that in AZ2 (44% for AGB, 54% for PNU, and 46% for NNI). The validation results confirmed that models for AZ1 (Figures 4A,D,G) had lower or similar RMSE and RE than models for AZ2 (Figures 4B,E,H) or AZ1+2 (Figures 4C,F,I), respectively.

For the PCC, the correlation and corresponding p-values were analyzed between candidate independent variables and dependent variables. Notably, there were both positive and negative correlations in each group (Figure 5). For the temperature-related variables, GDD was selected first due to its highest correlation with AGB in AZ1, and RAD was selected for its high correlation with AGB but lower correlation with GDD (<0.5). Similarly, PPT and SDI were selected from the water-related variables, while BN and SR from the management variables (Figure 5A). The above steps were repeated for AGB, PNU, and NNI in AZ2 (Figure 5B) and AZ1+2 (Figure 5C), respectively. Overall, six variables were selected in both AZ1 and AZ2 and five variables in AZ1+2.

For the VIRRFR, cross-validation accuracy improved significantly within the top two or three variables and then fluctuated slightly (Figure 6). Variables were determined by taking the highest R², lowest RMSE, and RE into account. For example, when the AGB model achieved the largest R² (0.71), the lowest RMSE (0.71 t ha⁻¹) and RE (25.3%) in AZ1, the 10 variables in the models were selected as the optimal variables (Figure 6A). In the same way, the number of selected variables for AGB (13 in AZ2 and 11 in AZ1+2, Figures 6B,C), PNU (5 in AZ1, AZ2, and 12 in AZ1+2, Figures 6D–F), and NNI (9 in AZ1, 11 in AZ2, and 12 in AZ1+2, Figures 6G–I) were determined.

The importance of variables selected by the PCC for estimating N status indicators was explored by the RFR models (Figure 9). The results showed that NDRE was consistently viewed as the most important variable except for the estimation of NNI in AZ1, and BN was considered important for all N status indicators. It should be noted that the important variables differed with the AZs and N status indicators. In AZ1, SDI and RAD were ranked as the top two variables for AGB, PPT and BN for PNU, and PPT and BN for NNI (except for NDRE). The water-related variables were consistently identified to be the vital variables (Figures 9A–C). Compared to AZ1, temperature-related variables were also of great importance to AZ2 (Figure 9D). Generally, water and temperature-related variables were listed as the most critical variables across AZs (Figures 9G–I). With respect to N status indicators, temperature-related variables were more important for AGB (Figures 9D,G), while water-related variables were more important for PNU and NNI (Figures 9B,C,E,F).

Different combinations of variables were further investigated and evaluated for the estimation of N status indicators based on dRMSE. The results indicated that the more the sources, the better the model's performance (Figure 10). Specifically, the dRMSE of two-source data (NDRE+C or NDRE+M) was more significant than that of single-source data (NDRE), and full sources of data achieved the best results, with 33.72–40.95% for AGB, 38.56–46.99% for PNU, and 48.54–52.50% for NNI. The combination of NDRE+C produced better results than those of NDRE+M (23.55, 27.21, and 42.83% vs. 13.62, 3.42, and −2.1% for AGB, PNU, and NNI, respectively), indicating that climatic data make more important contributions to the estimation of N status indicators in AZ1. On the contrary, the combination of NDRE+M performed better than NDRE+C in AZ2. Adding the climatic variables produced the greater dRMSE for AGB in AZ2 than that in AZ1, while the opposite was true for NNI. For PNU, the performances of adding the climatic variables were similar in AZ1 and AZ2. In AZ1+2, the contributions of the climatic data to NNI were greater than PNU.

The SLR models are commonly used to estimate plant N status and guide variable N applications (Li et al., 2019; Wang et al., 2019). However, models using NDRE could only explain 48–68% of N status indicators variabilities for AZ1, 36–38% for AZ2, and 52–59% for AZ1+2 based on the test dataset (Figure 4). Those results agreed with the findings of wheat N status estimation only using VIs across farmer's fields in a village (Chen et al., 2019). The performance difference of these models between AZ1 and AZ2 may be attributed to the different experimental treatments (Table 1) and climatic and management factors (Figure 1B). Consequently, it is necessary to combine climatic and management factors with NDRE in RFR models to improve the estimation accuracy. In agreement with Cummings et al. (2021) and Wang et al. (2021), the multi-variable RFR models performed better than SLR models (Figure 3), because the RFR models included multi-source data and could analyze nonlinear and complex relationships (Wang et al., 2021).

Generally, variable selection prior to model construction could reduce the data dimension and measurement while differentiating meaningful variables from the noise and improving calibration efficiency (Heremans et al., 2015; Feng et al., 2020). Two variable selection strategies, the PCC and VIRRFR, were adopted in this study. In Figure 7, the performance of full variables models was not better than those of using fewer variables with variable selection, demonstrating the effectiveness of variable selection. Similar results were observed in the previous studies (Cai et al., 2019; Zhang et al., 2019; Wang et al., 2020). The PCC showed more potential than VIRRFR to estimate N status, except for the estimation of PNU in AZ1 and AZ2. The variables were selected based on the PCC strategy from a statistical analysis point of view, resulting in more informative and diverse variables. Many variables were correlated with each other at the p < 0.001 significance level but still selected in the VIRRFR strategy (data not shown), such as T_max and T_mean of AZ1 (Figure 5).

Consistent with prior studies, NDRE, a useful proxy reflecting plant growth and development, has been deemed the most important variable for crop N status estimation (Osco et al., 2020; Colaço et al., 2021). Nonetheless, it is climatic and management factors that affect winter wheat growth and development. Accordingly, this study focused on the values of climatic and management factors for wheat N status estimation in different AZs. It is not surprising that these factors had different contributions in different AZs and for different N status indicators. Water-related variables played a more important role in AZ1, while both water-related and temperature-related variables were important in AZ2 (Figures 9D–F). Indeed, unlike AZ1, with a dry climate and sufficient sunshine, AZ2 belongs to a subtropical monsoon climate and is characterized by abundant precipitation and relatively less sunshine. In addition, the contributions of temperature-related variables to improve the estimation accuracy of AGB were greater than those of water-related variables, which were more correlated to PNU and NNI. A possible reason is that temperature has a greater effect on crop yield that is determined by the accumulation and redistribution of dry matter (Lee and Tollenaar, 2007) than rainfall (Schlenker and Lobell, 2010), because the dry matter of plants is formed by intercepting solar radiation and physiological and biochemical processes (Xue et al., 2002). Another possible reason is that the crop root system influenced by the water directly or indirectly is critical to absorbing nutrients from the soil, like N (Walsh et al., 2012; Sharma et al., 2018).

Overall, the results substantiated the viewpoint that integrating more sources of data will lead to better estimation performance (Chen et al., 2021; Gua et al., 2021; Liu et al., 2021; Li et al., 2022). Regarding their combinations, the results varied with sites and N status indicators. As for AGB, adding climatic factors improved the model performance more in AZ2 than in AZ1. Metabolic imbalances induced by lower average temperature during the vegetative stages could retard crop germination and plant growth (Ye et al., 2020). The difference of phenological information (SE stage, Table 1) would change AGB accumulation because it could, directly and indirectly, influence photosynthesis and respiration (Gua et al., 2021). On the contrary, adding the climatic factors improved the estimation of NNI in AZ1 more than AZ2, which could be caused by the influence of the length of the growth period or the responses of winter wheat cultivars to climate variability (Tao et al., 2014). In AZ1+2, the contributions of the management factors to PNU were greater than climatic factors because PNU showed a positive correlation to N rates (Egan et al., 2019; Sandaña et al., 2021). The sowing date in AZ2 was later than that in AZ1, leading to lower AGB accumulation, the number of tillers and plant height due to low temperature in the vegetative stages and higher temperature in the reproductive stages (Fazily, 2021). In contrast, a suitable sowing date can make full use of natural resources, such as light, heat, and water, enhancing plant population and facilitating the accumulation and uptake of N before winter. Different variables could interact with each other to influence crop growth and N status. For example, higher seeding rates could make up for stagnant tiller development, which would benefit cultivars with fewer tillers (Staggenborg et al., 2003).

In this study, the estimation models across two AZs showed similar results to that of a single AZ. Remarkably, it demonstrated the robustness and application of the established RFR approach.

While the results of this study were encouraging, there is still room for improvement. A total of 13 variables were considered in this study, which may be insufficient for constructing a reliable model for wider applications due to limitations in data amount and quality, model representativeness, and the number of key variables (Chlingaryan et al., 2018). Although the active canopy sensors are promising for evaluating N status, as indicated in this study, unmanned aerial vehicles and satellite remote sensing are more efficient for larger regions or farms (Huang et al., 2017; Chen et al., 2019). In this regard, follow-up research is needed to expand the scale of the study regions, supplement more climatic and management factors that directly or indirectly influence winter wheat growth, and develop more applicable N estimation models using unmanned aerial vehicles or satellite remote sensing images.