Changes in soybean cultivars released over the past 50 years in southern Brazil

On-farm soybean yield has increased considerably in the last 50 years in southern Brazil, but there is still little information about how selection and breeding for yield increase has changed the agronomic attributes of cultivars. The objectives of this study were to evaluate the changes in soybean yield, seed oil and protein concentration, and changes in plant attributes that might be associated with yield improvement of 26 soybean cultivars released over the past 50 years in southern Brazil, sown simultaneously in a common field environment for two growing seasons. The average rate of yield gain was 45.9 kg ha−1 yr−1 (2.1% ha−1 yr−1), mainly due increased seed number per area and harvest index. Over year of cultivar release, cultivars became less susceptible to lodging, as well as plant mortality reduced. Meanwhile, the seed oil concentration increased, and seed protein concentration decreased, which could have negative consequences for soybeans use and requires further attention for breeding of future cultivars. Breeders have successfully contributed to the annual rate of soybean yield increase in southern Brazil. By our results, as well as the official on-farm production data, there is no evidence of soybean yield reaching a plateau in the near future in southern Brazil.

the increase of number of seeds per area increased with cultivar year of release while 100-seed weight showed no significant trend 4,13 . These studies show evidence for regional variation in attributes that contributed to yield gain.
The genetic gain analysis enables comparison of the gains obtained in different environments or with the use of different breeding strategies 18 . Evaluating key attributes that increase soybean yield over time is strategic for planning yield increases in future cultivars. From the 2015/2016 until 2017/2018 growing season, the average soybean yields obtained for Brazil, southern Brazil and the Paraná state were similar: around 3.3, 3.3 and 3.4 Mg ha −1 , respectively 2 . These averages are close to those reached in the United States, 3.4 Mg ha −1 and slightly higher than those reached in Argentina, 2.8 Mg ha −1 , respectively 1 . In the last growing seasons, the national soybean yield contest reached yields of 7.1 to 8.9 Mg ha −1 with national champions commonly from the southern part of the country, especially from Paraná state 19 . Such results obtained in soybean yield contest show that there is a room to double the present Brazilian soybean yield by adopting the technology already available to the farmers 20 .
The objectives of this study were to evaluate the changes in soybean seed yield, and seed oil and protein concentration from 26 Brazilian soybean cultivars released over the past 50 years of cultivation. Changes in plant attributes associated with yield and seed constituents were examined since the cultivars were now grown simultaneously in a common field environment for two growing seasons in southern Brazil.

Results
During crop growth, the mean temperature was similar in the two growing seasons (20.6 and 20.3 °C) (Fig. 1A,B). Accumulated solar radiation was also similar in the two growing seasons (2584 and 2537 MJ m −2 ; Fig. 1C,D). For southern Brazil, rainfall accumulation around 800 mm is enough to maximize soybean yield 21 . Cumulative Relationship between plant attributes and year of cultivar release. In both seasons the increasing yield trend with YOR was associated with more seeds per area as consequence of more pods per area and more seeds per pod, besides a slightly greater 100-seed weight. The 100-seed weight was positively correlated with the soybean cultivar YOR in both growing seasons, increasing by 70.4 mg yr −1 (Fig. 2B). The number of seeds and pods per area were also positively correlated with soybean cultivar YOR in both growing seasons, increasing by 20.7 seeds per m 2 and 6.5 pods per m 2 (Fig. 2C,E). The number of seeds per pod also increased with YOR by 0.0075 seeds pod −1 yr −1 (Fig. 2D). Harvest index increased with YOR in both growing seasons by 0.0025 units per year (Fig. 2F).
The aboveground biomass was positively correlated with YOR in both growing seasons, increasing by 3.32 g m −2 yr −1 (Fig. 3A). The remaining plant density at full maturity (R 8 ) had a positive correlation with YOR, in both growing seasons (Fig. 3B), while plant mortality from emergence to harvest decreased over YOR, i.e., more plants that emerged reached full maturity (Fig. 3D). The number of empty pods was negatively correlated with YOR in both growing seasons, reaching values close to 0 in the most modern cultivars (Fig. 3C). Lodging scores also decreased with YOR by 0.07 units per year (Fig. 3E).
The number of lateral branches per area at harvest was positively correlated with YOR increasing by 0.56 branches per m 2 yr 1 (Fig. 4A). Plant height was not significantly correlated with YOR (Fig. 4B). The node number on the main stem per area presented a positive correlation with YOR (Fig. 4C). The node number on lateral branches per area decreased with YOR by 2.80 nodes m −2 yr −1 (Fig. 4D). The number of nodes per area and the height of the lowest pod showed no consistent trend over the YOR (Fig. 4E,F).
The seed protein concentration showed a negative correlation with soybean cultivar YOR, decreasing by 0.84 mg g −1 yr −1 (Fig. 5A). The seed oil concentration was positively correlated with soybean cultivar YOR, increasing by 0.53 mg g −1 yr −1 (Fig. 5B).
Relationships between plant attributes and yield. Seed protein concentration was negatively correlated with yield ( Fig. 5C), while seed oil concentration was positively correlated with yield (Fig. 5D).
For both growing seasons grain yield was positively correlated with seeds per area, pods per area, harvest index, 100-seed weight, aboveground biomass, seed oil concentration, nodes in main stem, plant density at harvest, seeds per pod and number of lateral branches. Yield was negatively correlated with lodging score, growth-cycle length, plant mortality, seed protein concentration, and empty pods per area (Fig. 6). Nodes in lateral branches, total nodes per area, height of lowest pod, and plant height had no consistent correlation with grain yield over the growing seasons (Fig. 6).

Discussion
In the past 50 years of soybean cultivation in Brazil, the average on-farm yield has increased in a linear fashion 1 . Our evaluation of old and modern soybean cultivars in the same environment showed that yields increased by 45.9 kg ha −1 yr −1 ( Fig. 2A) in two harvest seasons, which was similar to the average on-farm yield gain obtained in Brazil, southern Brazil and Paraná state of 44, 43 and 39 kg ha −1 yr −1 , respectively 2 . This shows that the release of improved cultivars has been a substantial driver for on-farm yield increases.
In this study, old and modern soybean cultivars were evaluated in the context of current agronomic practices, since in order to determine the genetic basis for increased yield, cultivars must be evaluated in a common environment 7,22 . Modern cultivars tend to have a shorter vegetative growth period (V E -R 1 ) compared to older cultivars (Fig. 7).
As management practices evolved over time along with the release of newer cultivars, the rate of on-farm yield increase is a result of the interaction between the cultivar and the environment. More studies are necessary to separate the effects of management from the effects of cultivar on increased on-farm yield; and to understand how older cultivars from Brazil perform under less favorable conditions that occurred more often in the past, e.g., low fertility, acid soils, or high aluminum concentrations in the soil. Modern cultivars from United States produce more than the older ones even in low yield environments, but especially in high yield environments 7 . For cultivars and environmental conditions from Brazil this is unknown and further studies with a historical set of cultivars under less favorable conditions would be needed to test the interaction.
The rate of yield gain in cultivars from southern Brazil is similar to those obtained for cultivars from Argentina (43 kg ha yr −1 ) 4 . Previous studies carried out in Paraná, southern Brazil, with soybean lines released from 1981 to 1986 indicated yield increases of 45 and 37 kg ha −1 yr −1 for early and semi-early maturity cultivars, respectively 14 . More recently, genetic gains of 39.4 and 40.7 kg ha −1 yr −1 were reported 17 . Several studies from major soybean producing countries have also reported yield gains, with different sets of historical cultivars and maturity groups ( Table 1).
The genetic improvement in yield of soybean cultivars from Brazil and Argentina were greater than those reported in the United States, China, Canada, and India. As commercial soybean production is more recent in South America, greater rates of yield increase show that breeders were quick to develop adapted and productive For cultivar breeding, high-yielding environments contribute to maximize the expression of genetic yield potential, even if the yield potential is not reached under farmers' field conditions 7,22 . In addition, the selection of cultivars in unfavorable environments, with lower yield potential, such as low fertility, high weed incidence and so on stressful conditions, helps to identify more resilient cultivars 31  www.nature.com/scientificreports/ affected this performance in modern cultivars compared to the older ones, especially the response to high temperatures and drought.
For the cultivars studied, the number of seeds per area, number of seeds per pod, number of pods per area, and 100-seed weight were all correlated to YOR and yield gain (Fig. 6), similar to results from cultivars from China 10,11,29 . However, for soybean cultivars from the United States and Canada, no consistent relationship between soybean cultivar YOR and 100-seed weight were found 7,12,22 . The major driver for soybean yield increases in our study was the number of seeds per area, being consistent with other studies 9,12 . Lodging was reduced over the YOR (Fig. 3E) and plant mortality was also reduced (Fig. 3D), made more plants survive from emergence until harvest (V E -R 8 ) and with upright canopy (Fig. 3B) which allowed greater formation of pods and seeds per area (Fig. 2C). www.nature.com/scientificreports/ The different plant mortalities between cultivars (Fig. 3D) were probably due to differences in canopy architecture. Modern cultivars had more upright canopies as evidenced by the lodging score (Fig. 3E), as well as in the past 50 years of cultivar release the increase in the harvest index was higher (1.52-fold, Fig. 2F) than the increase in aboveground biomass (1.19-fold, Fig. 3A). As diseases and pests were adequately controlled during the study, possibly the reasons for the different lodging score and plant mortality are due to competition within the soybean canopy caused by attributes such as stem strength, root support, and/or canopy weight (excessive vegetative growth), but this should be further evaluated in future studies.
The increase in yield potential in soybean cultivars from United States is associated with more efficient light interception and energy conversion efficiency 6,22 . In modern cultivars from United States the yield is strongly associated with the number of seeds per area and positively correlated with crop growth rate 32 .  www.nature.com/scientificreports/ The ability to increase the aboveground biomass up to the end of seed filling stage (R 5 ) is a critical component for increasing the number of seeds per area, since the seed filling rate as well as the 100-seed weight show little or no variation 32 . Experiments performed with elevated [CO 2 ] concentration suggest that the capacity of the source (photoassimilates) is greater than that of the sink (pods and seeds), as soybean cultivars showed increases in leaf photosynthesis of 24%, while yield increased only by 15% and the harvest index decreased 33,34 . In addition, there is significant genetic variation in soybean response to elevated [CO 2 ] concentration 35 , so it seems that in the future more productive cultivars need a greater reproductive sink, that is, more nodes, pods, and seeds.
In this study, the height of the lowest pod from both old and modern cultivars was higher than 10 cm and was suitable for mechanical harvest (Fig. 4F). For some cultivars, the height of the lowest pod was around 20 to 30 cm, having infertile nodes in these lower canopy layers. The production of pods and seeds in these nodes could be necessary to improve cultivar yield potential, being a desirable attribute. Soybean cultivars with plant architecture that allow greater light distribution within the canopy have higher yield, mainly due to high pod fixation and seed filling rate 36 .
Aboveground biomass presented a slightly increased trend with YOR in both growing seasons (Fig. 3A). However, some old cultivars, produced as much aboveground biomass than some of the modern cultivars, despite having lower seed yield, which occurred in part due the increase trend in harvest index over YOR (Fig. 2F), and by the longer vegetative growth period in some old cultivars (Fig. 7).
The aboveground biomass positive correlation with YOR (Fig. 3A) and negative correlation with lodging (r: − 0.39, significant at p ≤ 0.001, data not shown) suggests that the net canopy photosynthesis of modern cultivars increased, but it should be further investigated. In our study plant lodging started mainly after the flowering (R 1 ) stage, when some plants overtopped others. Similar results were found for soybean cultivars from United States, where greater photosynthetic efficiency was achieved by the more upright plants with less susceptibility to lodging 6 .
The number of empty pods decreased with YOR in both growing seasons (Fig. 3C). This shows that modern cultivars have greater capacity to fill seeds and/or that modern cultivars are more efficient at aborting pods that they are not able to fill. In soybeans from the United States, the modern cultivars under low plant densities www.nature.com/scientificreports/ produce more compensatory yield in lateral branches than older cultivars 24 . In this study, there was a slight trend to reduce the number of lateral branches with YOR, but this needs to be further investigated in lower plant densities.
No-consistent trend was observed between plant height and YOR in our study, although there was a great difference among cultivars (Fig. 4B). Some cultivars from China also did not show any trend in plant height with YOR 11 . However, a negative relationship between plant height and soybean cultivar YOR was observed in studies conducted with cultivars from the United States 5,7,8 , Canada 12 and China 9 . In the search for cultivars more resistant to lodging, breeders opted for shorter plants in the United States which increase seed production per area and harvest index 8 .
The node number on the main stem had a slight upward trend with YOR (Fig. 4C). This attribute is strongly related to temperature, photoperiod, and the cultivar growth habit 37 . The node number on lateral branches had no consistent trend with YOR over growing seasons (Fig. 4D). The non-significant trend occurred because it is a complex trait and because in our study, we used cultivars with both determinate and indeterminate growth habit to represent the cultivars that predominated over time, as the use of determinate growth habit predominated in the past and indeterminate growth habit increased and prevails in the last 20 years 3 .
The lodging resistance is an attribute that is positively related to YOR and yield in Brazil (Fig. 3E). It corroborates with results obtained with cultivars from Argentina 4 , United States 5-8 , Canada 12 , and China 9,10 . Upright plants are more able to intercept photosynthetically active radiation, especially during the seed filling period 6 . The harvest index consistently increased with soybean cultivar YOR in this study (Fig. 2F), which also occurred in cultivars from the United States 6 and China 9 .
From this set of cultivars, the seed protein concentration had a negative relationship with seed yield in the average of growing seasons (r: − 0.70, significant at p ≤ 0.01) in the same way as it occurred with cultivars from the United States 7 . These information shows that for future cultivars, breeding must address the seed composition desired by the consumer market, or by industry, as it is a value aggregator to soybean production.
Selection for greater yields has come at the cost of seed protein concentration in cultivars from different countries. In our study seed protein concentration reduced at a rate of 0.84 mg g −1 yr −1 (Fig. 5A) and showed a significant negative correlation with yield (Fig. 5C). A reduction in seed protein concentration was also reported in the United States between rates of 0.16 to 0.22 mg g −1 yr −17 and of 0.18 to 0.35 mg g −1 yr −18 . In Canada, a reduction in seed protein concentration was also reported at a rate of 0.54 mg g −1 yr −112 . However, in Chinese cultivars, the seed protein concentration was not consistently modified in the modern cultivars in relation to the old ones, although the rate of yield improvement has been smaller in China compared to South and North America 9,11 .
Seed oil concentration increased significantly with soybean cultivar YOR at a rate of 0.53 mg g −1 yr −1 in this study (Fig. 5D), it also increased in cultivars from China between rates of 0.17 to 0.60 mg g −1 yr −111 , and in cultivars from United States between rates of 0.05 to 0.14 mg g −1 yr −17 and of 0.09 to 0.27 mg g −1 yr −18 . In Canada, an increase in seed oil concentration was also reported at a rate of 0.45 mg g −1 yr −112 . In another study with cultivars Table 1. Summary of the soybean yield gain results, authors, country, environments (E) and/or seasons (S), latitude, years of breeding period (BP), maturity group (MG), number of cultivars (n), and yield gain (YG). www.nature.com/scientificreports/ from Northeast China, no increase in seed oil concentration was observed with the YOR, but the cultivars used in the past had already higher seed oil concentration, around 200-220 mg g −19 .
In the soybean breeding program, there is no unique path to increase yield through trait improvements as there is a complex interaction among traits 29 . Although seed concentration of oil and protein, plant height, and 100-seed weight are important attributes, the primary considerations in the decision to release a new cultivar are yield, days to maturity and lodging 7 .
For the subtropical conditions such of this study, maintaining the plant density that emerged until the harvest period (V E -R 8 ) proved to be an important attribute. The sample of cultivars used in this study also evidences that breeding improved the canopy architecture of the most modern cultivars, making them less susceptible to lodging.

Conclusions
Over the last 50 yr. breeders have successfully contributed to the annual rate of soybean yield increase in southern Brazil. The average yield gain rate evaluated in this study, 45.9 kg ha −1 yr −1 , is close to the on-farm rate of yield increase, showing a similar linear trend, which illustrates the important role of genetic improvement in Brazilian soybean production. Increased seed number per area and harvest index were the main contributors to increased yield for the evaluated cultivars. Plants more resistant to lodging provided a more upright canopy, which increased the number of seeds per area. Also, the reduced lodging over the year of cultivar release while plant mortality also reduced, made more plants survives from emergence until harvest and with upright canopy, which allowed greater formation of pods and seeds per area. The seed oil concentration increased, and seed protein decreased with breeding, which could have negative consequences for the use of soybeans and requires further attention in the development of future cultivars.

Methods
Site and experimental design. Two field experiments were carried out during the 2016/2017 and 2017/2018 growing seasons at the Universidade Estadual do Centro-Oeste research site, in a subtropical environment located in Guarapuava, Paraná State, Brazil (25° 23′ S, 51° 29′ W, altitude 1029 m). By convention, we adopt the year of harvest to refer to each growing season. The soil in the area is classified as very clayey Oxisol (USDA Soil Survey). The climate of the location is classified as Cfb by Köppen's climate classification system 38 .
In both growing seasons, predecessor crops followed a common rotation scheme from commercial fields. The experimental site for the 2016/2017 growing season had been cultivated with potatoes in the previous summer and black oat in the winter, in a conservation tillage system. The experimental site for the 2017/2018 growing season had been cultivated with maize in the previous summer and black oat in the winter, in a no-till cultivation system. Based on soil chemical analysis (Table 2), limestone was applied at a rate of 850 kg ha −1 and 1000 kg ha −1 in the winter season of 2016 and 2017, respectively.
The experimental design consisted of completely randomized blocks with three replications of 26 historical soybean (Glycine max (L.) Merrill) cultivars (Table 3). These cultivars, released from 1965 to 2015, were selected based on a survey with farmers, discussion with researchers and from the scientific literature.
In the past, the relative maturity of cultivars in Brazil were classified as early, mid, and full-season based on location 39 , but it was not successful in describing maturity for different latitudes and environments that occur in the soybean production region of Brazil 39 . This classification was gradually replaced from the 2000s onwards by the maturity group (MG) method, which groups cultivars based on photoperiod responsiveness and adaptation area, later introduced mainly by foreign soybean breeding companies 40 . During the choice of cultivars, we selected early and mid-season growth cycles between those released before the 2000s, and group V and VI for those launched after the 2000s (Table 3).
Fifteen seeds per old cultivar were obtained from germplasm bank of Brazilian Agriculture Research Corporation and were multiplied prior to sowing to have enough seed number, as well as high seed physiological potential. Modern cultivars were obtained from breeder companies or seed suppliers. The plant experiments were performed in accordance with relevant guidelines and regulations.
Cultivar plots contained 4 rows, each 5 m long with 0.45 m row spacing. The experimental area was limited to the two central rows, excluding 0.5 m from the edges. Before sowing, seeds were treated with Pyraclostrobin Seeds were sown on November 4 in both 2016 and 2017 in adjacent field areas. This sowing date was within agroclimatic zoning that benefits potential soybean yield in this subtropical environment 21 . Basic fertilizer application consisted of 80 kg ha −1 of P 2 O 5 , 80 kg ha −1 of Ca, 53 kg ha −1 of S (single superphosphate) and 70 kg ha −1 K 2 O (potassium chloride).
The seeding rate interval recommended by breeders was different between cultivars and we tried to adjust it within this recommended range, aiming for the best condition for the development of the cultivars. Approximately 44 seeds m −2 were sown and were thinned out during the V C -V E growth stages 41 to fit within the recommended plant density (Table 3). Weeds were controlled by herbicide before crop emergence and mechanically removed over the growing season. Pests and diseases were adequately controlled.
Evaluations. Four plants per plot were harvested at the grain filling (R 5 ) growth stage when each cultivar reached this stage. Leaf and stem biomasses were dried in a forced air drier at 60 °C for 48 h.
An area of 3.6 m 2 per plot was harvested at the full maturity (R 8 ) growth stage to determine yield, 100seed weight, number (nº) of seeds per area, nº of seeds per pod, nº of pods per area, height of the lowest pod, www.nature.com/scientificreports/ aboveground biomass, plant density at harvest, nº of empty pods per area, plant mortality, nº of lateral branches per area, plant height, node nº on the main stem per area, node nº on lateral branches per area, lodging score, harvest index, seed protein concentration, and seed oil concentration. The 100-seed weight was evaluated from a subsample of 600 seeds. The height of the lowest pod was measured from the first pod insertion point to soil surface. The aboveground biomass used was the sum of dry leaf biomass at R 5 , plus dry biomass of stem, seed, and pod shell at R 8 . Seed moisture was adjusted to 130 g kg −1 . Plant mortality was assessed by the relationship between plant density at harvest and the initial plant density. Death plants were considered the plants that did not reach full maturity, i.e., produced no seeds or pods at all. The Lodging score was evaluated through visual qualitative evaluation with grades from 0 to 5, where 5 indicates 100% lodging and 0 no lodging at all. Apparent harvest index was determined by dividing seed mass by aboveground biomass at R 8 (stem, seed, and pod). Seed samples were finely ground using a plant mill prior to seed oil and protein evaluation. Seed nitrogen concentration was determined by indophenol blue spectrophotometric method 42 after sulfuric acid digestion (digestion block). Seed protein concentration was calculated as seed nitrogen concentration × 6.25 and expressed on a dry weight basis. Seed oil concentration was determined using the Soxhlet extraction technique with the solvent petroleum ether pro analysis (Method 945.16 from AOAC 43 ). Seed protein and oil concentrations are expressed on a 130 mg g −1 moisture basis.  www.nature.com/scientificreports/ Meteorological data. Daily meteorological data including temperature (Fig. 1A,B), solar radiation ( Fig. 1C,D), and rainfall (Fig. 1E,F) were obtained from a meteorological station (SIMEPAR/Brazil) located around 100 m far from the experiments. A sequential water balance 44 was calculated to identify phases with water deficit during the crop growing season (Fig. 1E,F).

Data analysis.
Analysis of variance (ANOVA) was determined (Table 4) using a mixed model with cultivar effect considered fixed, with a random intercept for growing season, block, and block within growing season, using the PROC MIXED procedure of SAS On Demand for Academics (SAS Institute, Inc., Cary, NC). The total variation (σ 2 p) for each treatment was partitioned into variance components -cultivar (σ 2 c), growing season (σ 2 gs), and cultivar × growing season interaction (σ 2 c × gs) variance-using the VARCOMP procedure of SAS. As σ 2 c were higher than σ 2 c × gs for the evaluated attributes, data were averaged across growing seasons according to Gomez and Gomez criteria 45 . Linear regression analysis was performed taking grain yield and evaluated attributes as dependent variable and year of cultivar release as independent variable using PROC REG procedure of SAS (Table 4). Each evaluated plant attribute was plotted against the year of cultivar release to illustrate their changes over time. Pearson's correlations among evaluated attributes by year of release or yield were calculated to establish relationships and the correlations significance was evaluated by the Student t-test (α = 0.05).