Field experiments were conducted in 2018–2019 and in 2019–2020 at the Xinxiang Comprehensive Experimental Station of the Chinese Academy of Agricultural Sciences, Henan Province, China (35°18′N, 113°51′E, 78 m a.s.l.). The region has a warm, temperate, semi-humid climate. The brown soil of the region had pH and bulk density of 7.11 and 1.38 g cm^-3, respectively. The water, organic carbon, and total nitrogen levels were 36.1%, 1.7 g kg^-1, and 1.11 g kg^-1, respectively. In 2018–2019 and 2019–2020, the accumulated temperatures were 2559.3 °C and 2437.1 °C and the precipitation levels were 105.5 mm and 110.6 mm, respectively (Ma et al., 2021).

Two hundred and nine varieties were selected among natural wheat populations and Huaimai 40 (stable) and Zhengmai 7698 (sensitive) served as the control varieties. They were sown on October 8, 2018 and October 23, 2019, respectively, and harvested on June 6, 2019 and June 2, 2020, respectively. The experiment was repeated twice using a randomized complete block design. The sowing density, row spacing, sowing depth, and plot size were 270 plants m^-2, 10 cm, 5 cm, and 4.2 m² (3 m × 1.4 m), respectively. Details of water and fertilizer management during the two growing seasons are listed in Table 1 . The SN, KN, GY, and TKW were measured during two growing seasons. Two random samples were taken, after the grain filling stage, from each sample frame (0.5 × 0.5 m) within the test plots (excluding edge areas). The SN in the sample frames was recorded, and the average value was obtained. Before the wheat maturity, 20 representative ears of wheat were randomly selected from each plot and threshed to obtain the average KN. After the wheat matured, it was threshed, dried, and weighed to obtain the GY per unit area (ha). Representative kernels (1000) were then selected and weighed to obtain the TKW. Pest and weed control were implemented in accordance with local productive field standards.

The present study included 209 wheat varieties sampled from different habits (winter, semi-winter, weak-spring, or spring), gluten content types (strong, medium-strong, medium, or weak), breeding years (1980s–2010s), and habitats (China or other countries). The materials were abundant and widely distributed. The medium gluten type had the highest proportion (59.14%) followed by the strong gluten (29.03%) and the weak and medium-strong gluten types (7.53% and 4.30%, respectively). Additionally, the proportion of the semi-winter type was the highest (65.24%) followed by the weak-spring type (17.68%) and the winter and spring types (12.80% and 4.27%, respectively). The crops cultivated in the 2000s had the highest proportion (50.82%) followed by those grown in the 1990s (24.04%) and those raised in the 1980s and the 2010s (14.21% and 10.93%, respectively). The Chinese varieties accounted for the highest proportion (89.00%) and were distributed across nine provinces. The varieties were distributed across six countries and accounted for only 11.00% of the total ( Figure 1 ).

DNA was extracted from the leaves of the wheat seedlings using the high-salt/low-pH method (Guillemaut and Maréchal-Drouard, 1992) and amplified by PCR. A kompetitive allele-specific polymerase chain reaction (KASP) fluorescence detector was used to determine the genotypes for seven functional genes that regulate KW (TaCKX-D1, TaGASR7-A1, TaSus1-7A, TaSus1-7B, TaGS5-A1, TaGW2-6A, and TaKG2-6B), considerably affect its characteristics, and are associated with KASP markers (Rasheed et al., 2019). The primer sequences and amplification conditions for each gene are described in Supplementary Table 1 . KASP detection revealed 209 varieties with 41 KW allelic combinations ( Supplementary Table 2 ). Their wide variety ensured genotypic diversity in this study.

Meteorological, soil, field management, and observation data were collected to calibrate and evaluate the parameters of the APSIM-Wheat model. The meteorological data included daily maximum and minimum temperatures, daily sunshine hours, and daily precipitation at the Xinxiang Agrometeorological Station between 1961 and 2020. The phenological stage dataset included emergence, three-leaf stage, tillering, elongation, booting, heading, anthesis, medium milk, and maturity at the Xinxiang station between 2001 and 2013. The soil data used to run the model and calibrate the soil parameters were obtained from Zain et al. (2021) and included soil water and nitrogen distribution during the winter wheat growth periods between 2017 and 2019. The field observation data used to calibrate and evaluate TKW, KN, and GY were derived from the phenotypic observation data (Normalized Difference Vegetation Index (NDVI), TKW, GY, KN, and SN) of the 209 varieties under different managements at the Xinxiang station between 2018 and 2020. They also included data for the shared socioeconomic pathways SSP2-4.5 and SSP5-8.5 of the five global climate models in the Coupled Model Intercomparison Project Phase 6, namely, BCC-CSM2-MR (China), CanESM5 (Canada), EC-Earth3-Veg (Europe), MIROC-ES2L (Japan), and UKESM1-0-LL (UK). Specific on-site farm managements are listed in Supplementary Table 3 . The physicochemical properties of the various soil layers are listed in Supplementary Table 4 .

As the study focused on the mechanisms by which KW responds to sowing date and genotype, 209 wheat varieties were screened here. First, it was assumed that varieties with KN and flowering period NDVI in close proximity had other genetic factors near the KW-regulating gene. Hence, varieties with similar KN (range < 5,000) and flowering period NDVI (range < 0.025) were selected. The data analysis was verified using the Kruskal–Wallis test. Eighty-one varieties met the foregoing criteria and did not significantly differ in terms of KN or NDVI. The Kruskal–Wallis test was run using the “npmc” package in R v. 3.6.2 (http://www.R-project.org/). The KN and NDVI for the 81 selected and 128 unselected varieties during their growth seasons are shown in Figure 2 . The allelic combinations and percentages are listed in Table 2 . The KW distributions for the selected 81 varieties in eight allelic combinations are shown in Figure 3 .

The eFAST method was used in this study to elucidate the relationships among the variety and ecotype parameters and accurately simulate KW, KN, and GY. The eFAST method has been used to assess the sensitivity of various models (Campolongo et al., 2000) and was applied in four steps.

Both the first-order sensitivity and the global sensitivity indices are used to assess the sensitivity of the system to the input parameters. However, the former only considers the effect of a single input parameter on the output parameter, while the latter considers the interaction among multiple input parameters. The closer the value of the sensitivity index is to 1, the greater will be the influence of this input parameter on the output parameter; the closer the value is to 0, the smaller will be the influence on the output parameter.

The phenological, soil, variety, and ecotype model parameters directly and indirectly affecting KW formation were calibrated and evaluated. The wheat phenological parameters were calibrated and evaluated via batch parameter adjustment ( Supplementary Table 6 ). The simulation was conducted under water- and nitrogen stress-free conditions. The physical properties were then calibrated and evaluated ( Supplementary Table7 ). Based on the results of the parameter sensitivity analysis, the variety and ecotype parameters related to KN were then calibrated and evaluated. The range of parameter values was obtained for the optimally simulated KN results for the 209 varieties. The variety and ecotype parameters related to KW and GY were then calibrated and evaluated. The range of parameter values was obtained for the optimally simulated GY and KW results for the 209 varieties ( Supplementary Table 8 ). Root mean square error (RMSE) was used to determine the total difference between the observed and simulated values. Linear regression (R²) between the measured and simulated values was used to evaluate model performance. Consistency (D) between the measured and simulated values was also determined.

Simulations were run to (1) explore the responses of KW to sowing date and genotype under the background of climate warming, and (2) screen adaptive sowing dates and favorable allelic combinations that enhance KW. The APSIM-Wheat model was used to simulate KW under various sowing dates, allelic combinations, and climate scenarios as follows:

R v. 3.6.2 (http://www.R-project.org/) was used to analyze and test the significance of the differences in KW under various sowing dates, allelic combinations, and climate scenarios. The AOV function in the “stats” package of R was used to perform the analysis of variance (ANOVA) and the lsd.test function in the “multcomp” package of R was used to perform the analysis of multiple comparison.

The first-order and global sensitivity indices of the factors affecting wheat KW, GY, and KN were consistent. Of the 41 parameters affecting KW, the variety parameters “max_grain_size” and “grains_per_gram_stem” had first-order sensitivity indices > 0.2 and global sensitivity indices > 0.4. The ecotype parameters “y_swdef_leaf” and “fr_lf_sen_rate” had first-order sensitivity indices > 0.05 and global sensitivity indices > 0.1. The ecological parameter “node_sen_rate” had a first-order sensitivity index > 0.025 and a global sensitivity index > 0.1. The ecotype parameters “y_rue” and “eo_crop_factor_default” had global sensitivity indices > 0.05. Of the 41 parameters affecting GY, the variety parameters “max_grain_size”, and “grains_per_gram_stem” and the ecotype parameters “y_swdef_leaf”, “oxdef_photo”, “node _sen_rate”, and “y_swdef_pheno_flowering” had first-order sensitivity indices > 0.05 and global sensitivity indices > 0.1. Of the 41 parameters affecting KN, the species parameter “grains_per_ gram_stem” and the ecotype parameters “y_swdef_leaf”, “oxdef_photo”, and “y_swdef_ pheno_flowering” had first-order sensitivity indices > 0.05. The species parameter “grains_per_gram_stem” and the ecotype parameters “y_swdef_ leaf”, “initial_root_depth”, “oxdef_photo”, and “y_rue” had global sensitivity indices > 0.1.

Figure 4 shows the validation of the simulated phenological growth stages of wheat including emergence, three leaves, tillering, elongation, booting, heading, anthesis, medium milk, and maturity, at the experimental base in Xinxiang, Henan Province, China between 2007 and 2013. The simulated and measured phenological stages were consistent. We compared the simulated and measured values for the growth stages between 2001 and 2007 to evaluate the parameters independently. The fitted linear equation was y = 1.007x - 0.618 and its slope was ~1. As RMSE = 7.5 d, the deviation between the measured and simulated values for the phenological stage was 8 d. The D values calculated from the measured and simulated values was 0.998. The R² was 0.993 and P < 0.05 ( Figure 5 ).

Figures 6 , 7 show the validation of simulated water and nitrate nitrogen content in the 0–80-cm soil layer at the experimental base in Xinxiang, Henan Province, China between 2017 and 2018. The simulated results for the water and nitrate nitrogen content were good for the surface soil layer but poor for the deep soil layer. This could be owing to a lack of accurately measured soil data. Using data obtained from the literature may not fully reflect the actual local conditions, resulting in inaccurate simulation results for the deep soil layer. We compared the simulated and measured values for the soil water and nitrate nitrogen content in 2018–2019 to evaluate the parameters independently. The fitted linear equations were y = 0.635x + 0.073 and y = 0.82x + 1.421 ( Figure 8 ). The RMSE were 0.06 mm mm^-1 and 2.892 mg kg^-1, respectively, and the D values were 0.669 and 0.844, respectively. The R² were 0.235, 0.548 (P < 0.05), respectively.

Figures 9A–C shows the validation of the simulated KW, GY, and KN and their comparison against the measured values for the 81 wheat varieties in 2018–2019. The simulated and measured KW, GY, and KN values were highly consistent. We also validated the KW, GY, and KN for the 81 wheat varieties in 2019–2020 ( Figures 9D–F ). The fitted equations were y = 0.784x + 7.94, y = 0.985x - 375.08, and y = 0.369x + 13146.67, respectively. The RMSE were 2.602g TK^-1, 689.772 kg ha^-1, and 2,008.485 KN spike^-1, respectively. The D values calculated from the measured and simulated values were 0.819, 0.802, and 0.613, respectively. The R² were 0.5, 0.673, and 0.167 (P < 0.05), respectively.

Based on the current sowing dates, varieties, and field managements at the Xinxiang station, the calibrated APSIM-Wheat model was used to simulate the phenological stages, biomass, and KW of winter wheat under the baseline (1991–2020), SPP2-4.5 (2031–2060), and SSP5-8.5 (2031–2060) climate scenarios. The simulation demonstrated a temporal increase in the average temperature during the winter wheat phenological stages, shortening of the latter, increasing biomass (reproductive stages) ( Figure 10A–D ), and decreasing KW ( Figure 11 ). For the preceding climate scenarios, the average temperature for winter wheat were 9.3°C, 10.4°C, and 10.6°C, respectively; the mean phenological stages 235.21 d, 228.93 d, and 227.20 d, respectively; the mean reproductive stages were 33.83 d, 33.81 d, and 33.73 d, respectively; the mean biomass values for these reproductive stages were 1239.50 g m⁻², 1374.99 g m⁻², and 1386.21 g m⁻², respectively; and the mean KW were 39.90 g TK-1, 39.58 g TK-1, and 39.65 g TK-1, respectively.

The calibrated APSIM-Wheat model was used to simulate winter wheat KW under the aforementioned climate scenarios. ANOVA was then performed based on the simulation results. Allelic combination, climate scenario, and sowing date extremely significantly affected KW (P < 0.001) while the interaction between allelic combination and climate scenario significantly affected it (P < 0.05) ( Table 3 ). Figure 12 shows the non-significant effect of changing the sowing date on TKW and KW under SSP2-4.5 and SSP2-8.5 scenarios. The KW was 39.58 g TK⁻¹ and 39.65 g TK⁻¹ for the normal sowing period on October 11 (DOY 284), respectively. The average KW increased gradually from October 5 to October 23 (DOY 278–DOY 296; an average increase rate of 0.092 g TK⁻¹ and 0.126 g TK⁻¹, respectively). Figure 13 shows that changing the sowing date for SSP2-4.5 and SSP2-8.5 scenarios had a significant effect on the accumulated temperature at the reproductive stages ≥10°C. The accumulated temperature at the reproductive stages (≥10°C) increased with delay in sowing, with an average increase rate of 1.57°C·d and 1.56°C·d, respectively. Figure 14 shows the favorable allelic combination B6 (TaCKX-D1b+ Hap-7A-1+ Hap-T+ Hap-6A-G+ Hap-6B-1+ H1g+ A1b) with a KW of 40.13 g TK⁻¹. For the SSP2-8.5 scenario, the favorable allelic combinations are B5, B6, and B8, with a KW of 40.53 g TK^-1, 40.85 g TK^-1, and 40.19 g TK^-1, respectively, and the KW values of allelic combination B6 (TaCKX-D1b+ Hap-7A-1+ Hap-T+ Hap-6A-G+ Hap-6B-1+ H1g+ A1b), from the first quartile (25%) to the third quartile (75%), from the box position, were higher than those of B5 and B8. Stable, favorable allelic combinations can mitigate the negative impact of climate change on wheat TKW.

Figure 15 shows the values of the species-type parameters “potential_grain_growth_rate,” “potential_grain_filling_rate,” and “max_grain_size” for different allelic combinations in the APSIM-Wheat model. In the B6 high-KW allelic combination, “max_grain_size” had the highest value with a mean value of 0.0503; the “potential_grain_growth_rate” parameter had the highest value (top 25%) with a mean value of 0.00128; the “potential_grain_filling_rate” parameter had the median value (front 50%) with a mean value of 0.00167.

Parameter localization is a prerequisite for effective crop model use. Parameter sensitivity analysis facilitates targeted model parameter calibration (Kirui et al., 2022). Process-based physiological and ecological models simulate crop growth, development, and water and nutrient uptake. They can explore the dynamic relationships among atmosphere, crop, and soil and enable the prediction of crop responses to climate, genotype, farm management, and their interactions (Briak and Kebede, 2021). Prior studies focused on crop variety parameters. However, Zhao et al. (2014) reported that the APSIM-Wheat model includes other variables that may also affect the output. To minimize model uncertainty, researchers performed a sensitivity analysis of eight variety parameters and 33 ecotype parameters. Kernel weight (KW) was sensitive to the variety parameters “max_grain_size”, “potential_grain_growth_rate”, and “potential_grain_filling_rate”. This finding was consistent with the results of Zhao et al. (2014). Nevertheless, KW was also sensitive to the ecotype parameters “y_swdef_leaf”, “y_rue”, “eo_crop_factor_default”, “fr_lf_sen_rate”, and “node_sen_ rate”. Other studies reported contradictory results (He et al., 2015; Casadebaig et al., 2016). Broadening the parameter range reduces the uncertainty of the model and improves its performance at simulating winter wheat growth and development (Silva et al., 2021).

Calibrating and evaluating APSIM-Wheat is the foundation for determining the optimal sowing time and allelic combination high KW in winter wheat under future climate scenarios (Ahmed et al., 2016). The present study validated the variety and ecotype parameters for KW and GY based on the validation of the model parameters phenological period, soil moisture content, nitrate nitrogen, and particle number. This analysis provides theoretical support for subsequent investigations of the source-sink relationships in wheat KW. Figures 5 , 8 show that after calibration and evaluation, the APSIM-Wheat model simulated phenological period, soil moisture content, and nitrate nitrogen reasonably well. The R² were 0.993, 0.124, and 0.286, respectively, and the RMSE were 7.588 d, 0.06 d, and 2.892 d, respectively. Hence, the APSIM-Wheat model lays a theoretical foundation for optimizing the parameters related to KW and GY. Substantial differences were observed in the sowing dates over 2018–2020, the weather conditions during the sowing periods, and the multi-year phenological data from the agrometeorological stations. Nonetheless, the calibrated and evaluated APSIM-Wheat model reflected the responses to different sowing dates. Figure 9D shows that the calibrated and evaluated APSIM-Wheat model accurately simulated KW (R² = 0.575; RMSE = 3.076 g TK^-1). The APSIM-Wheat model had an input of 209 varieties and accurately simulated the phenological period and the soil water-nitrogen balance. Thus, the calibrated and evaluated model reflected the responses to diverse allelic combinations and is an important tool for selecting optimal sowing dates and allelic combinations in future studies.

The APSIM-Wheat simulation-based study revealed that climate warming is an important cause of KW reduction in wheat. KW is influenced by biomass accumulation during vegetative growth in winter wheat. It is also constrained by grain size and filling rate (Almeida and Lidon, 2009; Liu et al., 2020). Juknys et al. (2017) reported that according to long-term phenological data, climate warming can shorten the nutritional growth stages of winter wheat and result in insufficient biomass accumulation during the pre-reproductive stages as well as inadequate KW during the post-reproductive stages. Jenner et al. (1993) demonstrated that elevated temperatures reduce KW by downregulating the enzymes biosynthesizing soluble starch. In this manner, high temperatures shorten grain filling. The APSIM-Wheat model simulation of the present study forecasted that the average temperature will increase and the phenological stages will decrease in future winter wheat growing seasons ( Figures 10A, B ). However, biomass will increase ( Figure 10D ). Hence, the KW source will increase and the sink (grain-filling period from flowering to maturity) influences the decrease in KW. Field and greenhouse experiments conducted by Xie et al. (2015) showed that faster and longer grain filling were associated with higher KW. In the present study, the APSIM-Wheat model parameters “max grain size”, “potential_grain_growth_rate”, and “potential_grain_filling_rate” reflected the post-flowering grain size and the grain filling and growth rates which determine KW. In APSIM-Wheat, variety parameters “max_grain_size”, “potential_grain_growth_rate”, and “potential_grain_filling_rate” reflected maximum grain size, potential growth rate, and filling rate under ideal climatic conditions, respectively, and are affected by climatic conditions (primarily temperature). In this study, the actual grain filling rate in APSIM-Wheat decelerated and the growth rate accelerated under the future climate scenario with an increase in the average daily temperature ( Figure 10A ). Lobell et al. (2012) reported that an increase in the average daily temperature increased the transpiration rate of the grains, leading to water loss and negatively affecting the grain filling rate. In addition, Barlow et al. (2015) reported that an increase in average daily temperature will increase the growth rate of the grain, leading to a reduction in the duration of grain filling. In this study, the reproductive growth phase of winter wheat was shortened under the future climate scenario ( Figure 10C ). Therefore, the deceleration of filling rate and shortening of filling duration caused by climate warming are responsible for the reduction in KW.

Slight changes in the sowing window have only a marginal impact on KW. Sowing date modification has the lowest cost of all farm management measures and is a key factor in KW. Sowing date is mainly under human control and its modulation may subject winter wheat to different weather factors such as temperature and accumulated temperature that affect KW. For example, Shah et al. (2020) stated that delaying the sowing date by four weeks lowered the accumulated temperature from a normal sowing time value of 514°C to a wintertime value of 226°C, negatively impacted early growth, and reduced KW. As the present study used a maize-wheat rotation system, the wheat had to be harvested in a timely manner to allow for maize sowing. Therefore, we used the APSIM-Wheat model to establish seven sowing windows that started on October 5 and were repeated every 3 d until October 23. A delayed sowing window slightly increased KW because the adjustments to the sowing window led to a slight increase in the accumulated temperature ( Figure 13 ). For this reason, limited modifications to the sowing window resulted in slightly increased KW. Zheng et al. (2016) implemented the CSM-CERES-Wheat model in southern Khuzestan, Iran to adjust the sowing window to start on October 25 and repeat every 10 d until January 5. They found that fewer adjustments to the sowing window had a negligible impact on KW.

The variety parameters and their relative importance in the APSIM-Wheat model are closely related to KW. In general, high KW is accompanied by the upregulation of functional alleles associated with favorable KW (Zhang et al., 2014; Chegdali et al., 2022; Tillett et al., 2022). The APSIM-Wheat model parameters “max_grain_size”, “potential_grain_ filling_rate”, and “potential_grain_growth_rate” had the strongest influences on wheat KW (Zhao et al., 2014). Our first-order and global parameter sensitivity analyses showed that “max_grain_size” had the greatest impact on wheat KW followed by “potential_grain_filling_rate” and “potential_grain_growth_rate”. Several studies demonstrated that “max_grain_size” adjusts the grain size, “potential_grain_filling_rate” adjusts the grain filling rate, and “potential_grain_growth_rate” adjusts the grain growth rate (Zhao et al., 2014; Xie et al., 2015; Kobata et al., 2018). The present study showed that the values for the parameters characteristic of winter wheat varieties with high KW were consistent with high KW allelic combinations. In the B6 high-KW allelic combination, “max_grain_size” had the highest value ( Figure 15 ). This observation may be explained by upregulation of the advantageous A1b allele of the KW functional gene TaGS5-A1 which promotes cell division and the formation of large kernels (Wang et al., 2016). In the B7 low-KW allelic combination, “max_grain_size” had one of the lowest values of all parameters (bottom 25%) ( Figure 15 ). This finding may be explained by the upregulation of the advantageous A1a allele of the KW functional gene TaGS5-A1. An earlier study showed that the KW functional genes TaSus1-7A and TaSus1-7B encoding sucrose synthases have the advantageous alleles Hap-7A-1/2 and Hap-T that promote starch sediment biosynthesis and facilitate grain filling (Hou et al., 2014). In the B6 high-KW allelic combination, the “potential_grain_growth_rate” parameter had the highest value (top 25%), and the “potential_grain_filling_rate” parameter had the median value (front 50%) ( Figure 15 ). This discovery may be explained by the upregulation of the advantageous allele Hap-6B-1 of the KW functional gene TaGW2-6B, which regulates the number of endosperm cells, promotes grain development, and is superior to the favorable allele “Hap-6B-2” (Qin et al., 2014).The foregoing reports revealed that the APSIM-Wheat model variety parameter values effectively captured KW expression in different allelic combinations.

Determining the responses of KW functional genes to climate warming facilitates the identification of high-KW allelic combinations. Here, multivariate ANOVA revealed that allelic combinations and climate scenarios extremely significantly influenced TKW (P < 0.001) while their interaction significantly affected it (P < 0.05). Under different climate scenarios, the B6 allelic combination had higher KW and was stable. The B7 allelic combination had lower KW than the other allelic combinations. The preceding results suggest that the values of the variety parameter determining high KW were consistent with the high-KW allelic combinations. The variety parameter values were highest for the B6 allelic combination which included the KW functional genes TaSus1-7A, TaSus1-7B, TaGW2-6B, and TaGS5-A1 and their favorable alleles Hap-7A-1, Hap-T, Hap-6B-1, and A1b, respectively. The B7 allelic combination had low variety parameter values but included the KW functional genes, TaSus1-7A, TaSus1-7B, TaGW2-6A, and TaGW2-6B, and their favorable alleles, Hap-7A-2, Hap-T, Hap-6A-A, and Hap-6B-2, respectively. A recent study revealed that in TaSus1-7A, Hap-7A-1 and Hap7A-2 had strong and weak effects on wheat KW, respectively (Sehgal et al., 2019). The favorable allele Hap-6B-1 of TaGW2-6B has a stronger impact on KW than the favorable allele Hap-6B-2. However, both had a greater effect on KW than TaGW2-6A (Hou et al., 2014). The favorable allele of TaSus1-7A, “Hap-1/2,” was measured in multiple environments with high TKW and stable expression of the trait (Qin et al., 2014).Therefore, the B6 allelic combination can achieve maximum KW. The present and preceding studies provide critical reference for breeding winter wheat and certain other cereal crops with high kernel weight under climate warning.