This study was conducted in Hebei Province at the Mazhuang Experimental Station (MZ, 37°47′52″N, 115°18′28″E) during the 2021/2022 wheat growing season and at the Malan Experimental Station (ML, 37°58′28″N, 115°12′2″E) during the 2021/2022 and 2022/2023 wheat growing seasons ( Figure 1 ). The MZ and ML experimental stations were located in the Piedmont Plain of Mt. Taihang and the Heilonggang Plain, respectively, at elevations of 53 and 37 m in the Haihe Plain (Hou and Hao, 2010). Both the Heilonggang Plain and the Piedmont Plain of Mt. Taihang have a warm temperate semi-humid continental monsoon climate, whereas the Piedmont Plain of Mt. Taihang has relatively good soil conditions (Asadi Zarch et al., 2015). The annual precipitation in 2021/2022 and 2022/2023 was 804.7 mm and 731.0 mm, respectively, with the wheat growing season receiving 174.1 mm and 181.9 mm of precipitation, respectively. The annual average temperatures for the 2021/2022 and 2022/2023 seasons were 17.2°C and 15.14°C, respectively, with the wheat growing season averaging 8.7°C and 9.4°C, respectively. During the 2022/2023 season, the lowest temperature reached −9.3°C ( Supplementary Figure S1 ). The organic matter, total nitrogen, alkaline-dissolved nitrogen, available phosphorus, and potassium in the study area during the experiments were 22.4 g kg^–1, 1.26 g kg^–1, 131.2 mg kg^–1, 25.8 mg kg^–1, and 122.2 mg kg^–1, respectively, in MZ, and 19.7 g kg^–1, 1.09 g kg^–1, 125.0 mg kg^–1, 15.32 mg kg^–1, and 101.3 mg kg^–1, respectively, in ML.

The field arrangement followed a single-factor randomized block design with three replicates. A total of 29 winter wheat cultivars, the main cultivars in the Haihe Plain in certain years, were selected as treatments. The experimental materials were provided by the Hebei Provincial Agricultural Biological Resources Preservation Center ( Table 1 ). The plot area was 30 m² (6 m × 5 m; row spacing: 15 cm). Basal fertilizer consisting of 120 kg ha^–1 N, 120 kg ha^–1 P₂O₅, and 120 kg ha^–1 K₂O was applied in ML and MZ at Zadoks Stage 0, respectively. Wheat seeds were sown in ML and MZ at Zadoks Stage 0 (Zadoks et al., 1974) with a seeding rate of 375 seed m⁻² and harvested in ML and MZ on 6 June and 10 June 2021 and in ML on 8 June 2022. Using the local irrigation system, irrigation was applied at the jointing at Zadoks Stage 31 and anthesis at Zadoks Stage 65 (Zadoks et al., 1974) with 90 mm of groundwater each time (Fang et al., 2017). Additionally, 120 kg ha^–1 of N was applied during the jointing stage with irrigation, and weed and pest control were performed regularly during the growth period.

The logarithmic mean Divisia index (LMDI) method, which is a well-known decomposition model, was used to decompose the contribution of grain yield variation in different years to spike number, kernel number per spike, and grain weight (Xie et al., 2019). The grain yield (ΔY) changes from years t−1 to years t was calculated using Equation 16, as follows:

where, ΔS, ΔG, and ΔK, respectively, represent the influence of the spike number, kernel number per spike, and kernel weight on the changes of grain yield, whose relevant decomposition factors between t and t−1 was calculated using Equations 17–20, as follows:

where

The genetic gain (GN) in grain yield and its components based on cultivar release date was calculated using Equations 21, 22 (OrtizMonasterio et al., 1997):

where y_i is the grain yield and its components of cultivar i, ln(y_i ) is the natural log of y_i , x_i is the year in which cultivar i was released, a is the estimated intercept of each equation, and μ is the residual error. The linear function (Equation 21) provides an estimate of the grain yield and its components increase based on cultivar release date in absolute terms (i.e., b measures grain yield and its components gain), while the logarithmic function (Equation 22) gives the relative grain yield and its components increase (dy_i /dx_i ) such that 100b estimates the percentage grain yield and its components gain per year. All parameters were estimated using a restricted maximum likelihood linear mixed effects model (RMLM). The RMLM was also used to evaluate the effects of cultivar replacement and planting location on the grain yield, yield components, and grain-filling parameters of winter wheat. Winter wheat cultivar nested in the founder parent as a random effect for a lack of independence within cultivar and founder parent between year of cultivar release, experimental year, and site. The RMLM was conducted in the “lme4” package of R (version 4.2.1) (Bates et al., 2015). The variances of each fixed variable, random variable, and unexplained factor were calculated using variance decomposition, which was performed in the “partR2” package of R (Stoffel et al., 2020). Additionally, redundancy analysis (RDA) was conducted using the “vegan” package of R (Oksanen et al., 2023), with the grain-filling parameters as independent variables, and cultivar replacement and experimental location as dependent variables. All figures were generated using ArcGIS 10.2, OriginPro 2021 and the “ggplot2” package in R (Wickham, 2016).

The results of the official cultivar approval showed that the yield of wheat increased, the spike number per unit area decreased, and the kernel number per spike and kernel weight increased over the analyzed timespan ( Figures 3A–D ). The LMDI results indicated that the spike number had a negative contribution yield increase, except in the 2020s. Kernel number had a greater positive contribution to yield before the 2000s, while the kernel weight contributed positively to yield after the 2000s. Notably, the contribution of kernel weight in the 2020s increased by 43.95% compared to the 1980s ( Figure 3I ). Similarly, field experiment results showed that the yield of main winter wheat cultivars in the Haihe Plain has increased from the 1950s to the present. This increase in yield was accompanied by a decrease in spike number, and an increase in kernel number per spike and kernel weight ( Figures 3E–H ). According to the LMDI results, except for the 2000s, the spike number consistently had a negative contribution to yield increase. The LMDI indexes for the spike number were −3,222.1, −486.7, −767.7, −1,642.5, 1,137.1, −1,183.44, and −2,858.01 from the 1960s to the 2020s, respectively. Kernel number had a greater positive contribution before the 2000s, and the LMDI indexes for the kernel number were 2,126.3, 289.9, 72.7, −333.7, 605.3, −105.7, and −491.9 from the 1960s to the 2020s, respectively. After the 2000s, kernel weight contributed positively to yield, and the LMDI indexes for the kernel weight were −282.8, 152.5, 1,384.4, −627.5, 310.2, 1,087.3, and 1,350.9 from the 1960s to the 2020s, respectively ( Figure 3J ).

In this study, the average grain yields from the 1950s to the 2020s were 6,311.1, 6,029.8, 6,299.1, 7,742.8, 7,662.2, 8,146.3, 8,536.5, and 9,014.2 kg ha^–1, respectively. The average spike number was 638.7, 501.7, 527.4, 509.5, 432.9, 461.1, 453.1, and 373.7 ×10⁴ ha^–1, respectively, from the 1950s to the 2020s. The average kernel number per spike was 32.8, 36.4, 35.9, 38.3, 39.9, 37.9, 39.1, and 39.4, respectively, from the 1950s to the 2020s. The average thousand kernel weight was 41.9, 39.6, 42.4, 43.5, 44.8, 44.9, 48.4, and 52.7 g, respectively, from the 1950s to the 2020s. The results of RMLM showed that grain yield significantly increased by 39.12 kg ha^–1 and kernel weight significantly increased by 0.15 g year^–1, respectively, spike number significantly declined by 24,350 hm^–2 year^–1 (p < 0.001), but the kernel number per spike did not significantly increase (p > 0.05) ( Figures 4A–D ). Genetic gain refers to the improvement in average genetic value in a population or the improvement in average phenotypic value due to selection within a population over cycles of breeding. In this study, genetic gain specifically refers to the annual percentage increase or decrease in yield and yield components. GN of grain yield, spikes per unit area, kernel number per spike, and kernel weight were 0.629%, –0.574%, 0.158%, and 0.332%, respectively ( Figures 4A–D ). Furthermore, the results of variance decomposition showed that the year of release accounted for 25.6% of the variation in grain yield, whereas planted location, unexplained factors, and experimental year accounted for 11.8%, 50.9%, and 7.38% of the variation, respectively ( Figure 4E ). This finding indicates that in addition to varietal replacement, cultivation conditions, experimental year, and other factors have a large impact on grain yield. Kernel weight was affected more by the year of release, accounting for 37.5% of the variation, and spike number was affected more by the experimental year, accounting for 53.2% of the variation. Kernel number was mainly affected by the experimental location, accounting for 15.7% of the variation ( Supplementary Table S1 ; Figure 4 ).

The logistic regression model accurately captured the trends in the wheat grain-filling process (R ² > 0.9) ( Supplementary Table S2 ). Kernel weight growth showed an “S”-shaped increase in kernel weight during the grain-filling process when visualized. With the replacement of cultivars, the kernel weight of modern cultivars was higher than that of earlier cultivars ( Supplementary Figure S2A ). The grain-filling rate increased first and then decreased. The grain-filling rate indicated that the grain weight initially increased until it reached T_max , the time to reach the maximum grain-filling rate, after which it began to decrease. The T_max advanced with the replacement of cultivars, which were 20.3, 19.1, 20.1, 19.7, 18.9, 18.9, 18.3 and 17.1 days for the 1950s–2020s, respectively. R_max also increased gradually, except in the 1950s, with values of 2.7, 2.2, 2.2, 2.1, 2.2, 2.2, 2.3, and 2.4 mg grain^–1 d^–1 in the 1950s–2020s, respectively ( Supplementary Figure S2B ; Supplementary Table S3 ).

RDA was used to determine the response of the grain-filling characteristics to the replacement of cultivars, planting locations, and experimental year. These three constraining variables explained 61.0% of the variation in grain-filling characteristics ( Supplementary Table S4 ). The results showed that RDA1 accounted for 51.4% and RDA2 for 7.3% of the total variation in the grain-filling characteristics ( Figure 5A ). From left to right (MZ to ML, 2021/2022 to 2022/2023), v ₁ , T ₃ , p, T ₂ , T _0.99, and T_max decreased as the cultivation ecological environment deteriorated, whereas v_mean, v ₂ , R_max, v ₃, and T ₁ increased ( Figure 5A ). Moreover, T _0.99, T ₃, p, v_mean , T ₂, v ₁, v ₂, and R_max values were positively correlated with the year of cultivar release, whereas T_max , T ₁, and v ₃ were negatively correlated with the year the cultivars were released ( Figure 5A ). RDA was also used to determine the response of grain yield and kernel weight to grain-filling characteristics ( Figure 5B ). Grain yield and kernel weight explained 72.3% of the variation in grain-filling characteristics ( Supplementary Table S4 ). The results showed that RDA1 accounted for 67.5% and RDA2 accounted for 4.8% of the total variation in grain yield and kernel weight, respectively, but RDA2 was not significant (p > 0.05). From right to left, T _0.99 , T ₂, T ₃, T_max , and v ₁ increased as grain yield and kernel weight increased, whereas T ₁ , v ₂ , v ₃ , R_max , and v_mean decreased.

To further demonstrate the relationship between wheat-filling parameters and year of release, RMLM was performed ( Figure 6 ). During the year of release, p, T ₂, T ₃, and v ₁ increased significantly, whereas T_max and T ₁ decreased significantly (p < 0.05). However, there were no significant relationships between T _0.99, R_max , v_mean , v ₂, and v ₃ and the year of release (p > 0.05) ( Supplementary Table S2 ). The results showed that the replacement of wheat cultivars in the Haihe Plain had an impact on grain-filling characteristics. Specifically, it has extended the active grain-filling period, particularly during the rapid filling stage, with durations of 5.7, 6.6, 7.1, 7.6, 7.7, 7.8, 8.2, and 8.0 days for the 1950s–2020s, respectively, and the late grain-filling stage, with durations of 8.9, 13.3, 15.8, 16.9, 17.7, 18.1, 19.0, and 20.0 days for the 1950s–2020s, respectively. Additionally, it has enhanced the average grain-filling rate, with values of 1.1, 1.0, 1.0, 1.1, 1.0, 1.0, 1.1, and 1.3 mg grain⁻¹ d⁻¹ for the 1950s–2020s, respectively. However, it has also shortened the duration of the initial grain-filling stage, with durations of 15.5, 13.6, 13.1, 12.3, 11.6, 11.4, 10.5, and 9.5 days for the 1950s–2020s, respectively.

Domestic and foreign studies have shown that grain yield gradually increases with the replacement of the wheat crop in general; however, conclusions regarding yield components have been inconsistent (Zhou et al., 2007b; Feng et al., 2018; Qin et al., 2018, 2019; Sakuma and Schnurbusch, 2019). For foreign wheat cultivars, Sakuma and Schnurbusch (2019) showed that wheat yield was strongly influenced by changes in kernel number per spike; however, the effects of spike number and kernel weight were not significant. However, other studies have shown that a steady increase in the 1000-kernel weight and HI leads to a sustained increase in wheat yield and higher yield potential (Lo Valvo et al., 2018). Haihe Plain is located in northern China. In Northern China, Feng et al. (2018) reported that a significant increase in wheat yield was caused by an increase in kernel number per spike and kernel weight, although there was little change in spike number per unit area. Zhou et al. (2007a) also showed that there was a gradual increase in kernel weight during the replacement of wheat cultivars. Qin et al. (2015) showed that the replacement trends of wheat cultivars differed under different cultivation environments. Therefore, although previous studies have agreed that the replacement of wheat cultivars resulted in increased grain yield, differences in the timescales of studies and cultivation environments prevented consistent conclusions on the contribution of yield components to yield increases, but the increase in wheat yield in the Haihe Plain is closely related to kernel weight. This study demonstrated that both yield and grain weight increased significantly with the replacement of wheat cultivars. However, there are inconsistent conclusions regarding the reasons for the increase in yield during the process of cultivar replacement (Zhou et al., 2007b; Feng et al., 2018). In our study, the increase in kernel number was not significant, whereas kernel weight increased significantly during the cultivar replacement process. In contrast, Philipp et al. (2018) found that kernel number was the main reason for the increase in wheat yield. This discrepancy may be due to the different statistical methods used for measuring kernel number: the study used kernel number per square meter (Gonzalez-Navarro et al., 2016), while we measured kernel number per spike. Additionally, wheat cultivars have strong regional characteristics, which may also contribute to different reasons for yield increases during cultivar replacement in different areas. Given the significant increase in kernel weight observed in this study, it is necessary to further investigate the factors contributing to the increase in kernel weight during cultivar replacement.

Phenotype characteristics of wheat are influenced by both genetics and the environment; therefore, while wheat yield increased significantly with the year of release, there was also some influence of environmental factors, such as temperature, precipitation, and solar radiation, and soil factors such as topsoil thickness, soil water content, organic matter, nutrients, and trace elements in different ecological regions, differentially affecting the yield of the same wheat cultivar planted (Motzo and Giunta, 2004; Lobell et al., 2009; Subira et al., 2014). In the Haihe Plain, ML was located in the Heilonggang Plain zone, and MZ was located in the Piedmont plain of Mt. Taihang zone; the ecological conditions of the Piedmont plain for wheat planting were better than those of the Heilonggang Plain (Wu et al., 2023), which was further verified by comparison of soil characteristics in this study ( Supplementary Table S5 ). Combining the two wheat growing seasons, the variance decomposition results showed that the spike number was more affected by the experimental year, which was due to the environmental factors of frost damage in the second year, which reduced the spike number compared to the first year, whereas the kernel weight was more affected by the year of release.

For the wheat crop in general, Madani et al. (2010) showed that with the year of release, the increase in kernel weight was due to an increase in sink activity in modern wheat cultivars. Álvaro et al. (2008) found that the increase in the size of the sink in modern cultivars led to an increase in the capacity of the wheat stem to transport substances to the grain, thereby limiting the source capacity and lowering the kernel weight of early cultivars. Near the Haihe Plain, Zhou et al. (2007a) reported a significant increase in winter wheat yield caused by an increase in kernel weight. The mechanism was then established by Tian et al. (2011) who showed that the coordination of the source–sink relationship and increased pre- and post-anthesis transport capacity were the main reasons for the steady increase in kernel weight. Therefore, the kernel weight of modern wheat cultivars has increased because of the coordination of the source–sink relationship and the enhancement of sink activity in the Haihe Plain. In this study, the grain-filling rate gradually accelerated with the replacement of wheat cultivars in the Haihe Plain, which was consistent with previous studies. However, the grain-filling rate of wheat in the 1950s differed from that of other cultivars. It had a fast grain-filling rate but a short grain-filling duration, which may be attributed to the limitation of source capacity (Smith et al., 2018). Therefore, this study focused solely on sink activity and did not consider source capacity, which should be addressed in future research.

Furthermore, Khodarahmi et al. (2023) found that the average grain-filling rate increased gradually with the year of release; however, there was no significant association between grain-filling duration and the year of release. Other studies have also shown that the extension of grain-filling days is the main reason for higher kernel weights and yields of modern cultivars (Wiegand and Cuellar, 1981). For example, Wang and Shangguan (2015) showed that the key to increasing kernel weight is extending the duration of rapid and late grain-filling. The grain-filling rate in the rapid period was three times faster than that in the initial and late periods and, thus, had the greatest contribution to kernel weight (Xie et al., 2015). Wang et al. (2021) also showed that an increase in the grain-filling rate in the late period and a reduction in its fluctuation were key to increasing kernel weight, and the wheat grain-filling process is dominated by the regional ecological environment, especially temperature and light. Grain-filling characteristics vary depending on environmental conditions (Mirosavljević et al., 2018).

In the present study, with year of release, the duration of the initial grain-filling period became shorter whereas the filling rate increased, and the time to reach the maximum grain-filling rate advanced. Therefore, although the duration of the initial grain-filling stage of modern cultivars was shortened, the increase in filling rate resulted in an increase in kernel weight during the initial grain-filling stage. In addition, the significantly longer durations of the rapid and late grain-filling stages resulted in a longer active grain-filling period, which increased the kernel weight. This may be related to the fact that the North China Plain is often subjected to hot and dry winds, and the early cultivars are prone to early senescence in the middle and late filling stages (Cai et al., 2022). Furthermore, the Heilonggang Plain exhibited more rapid rates of grain-filling, including maximum and average rates, as well as during the rapid and late stages of grain-filling. In contrast, the Piedmont plain of Mount Taihang boasted longer durations for grain-filling periods, such as the effective grain-filling duration, the active grain-filling period, and the durations of both the rapid and late grain-filling stages. In addition, because of the lowest temperature of −9.3°C in the second year, which was close to the temperature of mild freezing damage (CMA, 2018), wheat suffered freezing damage, resulting in a shorter duration of filling period and active filling period in the second year. However, compared with the first year, the average and maximum filling rates increased. Beche et al. (2018) investigated the duration of developmental phases and eco-physiological traits related to grain yield in Brazilian wheat cultivars released across different decades, findings that suggest that longer durations in the growth and development stages can improve grain-filling characteristics and yield.

Kernel weight is one of the determinants of grain yield (Olmedo Pico et al., 2023). Maeoka et al. (2019) showed that the kernel weight of modern wheat cultivars increased by 27.28% and the spike number decreased by 14.83% compared with those of older cultivars. Similarly, Zhou et al. (2007a) showed that the kernel weight of modern cultivars in China increased by 49.31% and spike number decreased by 27.96% compared to those of earlier cultivars. Zheng et al. (2011) also showed that in the northern wheat region, the genetic gains for spike number, kernel number, and kernel weight were –0.62%, 0.75%, and 0.96%, respectively. In this study, kernel weight of modern cultivars reached 53.9 mg for Malan1, which was 26.7% higher than that in 1950s cultivars, but the spike number decreased gradually with the year of release, while modern cultivars decreased more rapidly by 62.6% compared to the 1950s. Sadras and Slafer (2012) showed that the number of spikes per unit area and kernel number exhibited greater phenotypic plasticity. Therefore, in the future, more attention should be paid to spike traits and canopy architecture.