Нужно ли нам собираться для региональной адаптации в Сойбе? - Оценка генотипного взаимодействия и стабильности трансвеститов в Сойбе в Германии.

1. Introduction

Soybean ( Glycine max (L.) Merr.) is a leguminous crop of globally high economic importance. Due to its high protein content, it is mainly used as animal feed. However, soybean products, such as tofu and dairy alternatives, are also very popular for human consumption with a steadily increasing demand in Germany [ 1 , 2 ] and worldwide [ 3 ]. In addition, soybean is also grown for its high oil content. As a crop that was domesticated in temperate China [ 4 ], soybean production requires warmer temperatures [ 5 ]. The largest producers are the USA, Brazil, and Argentina [ 6 ]. As a result of selection and breeding efforts, cultivars that mature earlier and are more tolerant to cooler temperatures are available nowadays [ 7 , 8 ]. This enables the successful cultivation of soybean at higher latitudes, like in Germany, despite the shorter growth periods and the generally lower temperatures, especially in spring. In Germany, the majority of the soybean acreage lies in the southern half and to date, there is only very little soybean cultivation in the more northern regions [ 9 ]. The southern regions are generally better suited for soybean cultivation due to temperature, rainfall patterns, solar radiation, and shorter photoperiod [ 5 ]. Photoperiod is a decisive factor in soybean cultivation since soybean is originally a short-day plant and as such requires short daylengths to induce flowering [ 10 ]. However, due to breeding efforts targeting this response, there are cultivars available that are less sensitive to long-day conditions, enabling soybean production at higher latitudes [ 11 ].

Germany has a high demand for soybean that is mainly covered by imports. Only around 2% of the amount of soybean imported into Germany is produced by German farmers [ 12 ]. In exporting countries, extensive soybean production leads to environmental issues such as the deforestation of rainforests [ 13 ]. In the context of the “Protein Plant Strategy”, with the aim to increase sustainable plant-based protein production and decrease the dependence on imports, the German government has promoted the expansion of the production of soybean and other legumes [ 14 ]. To increase the acreage and yield of soybean grown in Germany and to expand the cultivation across the whole country, especially to regions at higher latitudes, further breeding efforts are necessary. To successfully develop cultivars suitable for the different regions, a thorough understanding of the patterns of genotype-by-environment interaction is crucial.

Genotype-by-environment interaction describes the phenomenon of different genotypes responding differently to different environments, meaning the relative performance of the genotypes changes between environments. More specifically, crossover interaction leads to rank changes of genotypes between different environments. In order to deal with genotype-by-environment interaction, trait stability can aid in assessing a genotype’s suitability across a number of environments. Trait stability can be estimated following the static or dynamic concept [ 15 , 16 ]. Static stability describes the ability of a genotype to show the same absolute level of performance regardless of the environment. A dynamically stable genotype, by contrast, shows the same level of performance across environments but always relative to the environmental mean. In consequence, the absence of genotype-by-environment interaction leads to dynamic stability [ 17 ].

Two approaches to managing genotype-by-environment interaction in a breeding program are to either avoid it by selecting genotypes that show stable and high performance across environments (wide adaptation) [ 17 ] or to exploit it by selecting different genotypes for different mega environments (specific adaptation), which are environmentally similar subsets of environments. Since high stability is generally associated with low levels of genotype-by-environment interaction [ 17 , 18 , 19 ], it can be challenging to breed for wide adaptation. Consequently, when aiming to breed for wide adaptation, low levels of genotype-by-environment interaction or at least non-crossover interaction are required [ 20 ].

Genotype-by-environment interaction, or when analysed within a single year also referred to as genotype-by-location interaction, is a very common phenomenon for quantitative traits such as yield [ 21 ], but little is known for soybean in high latitude regions such as Germany in Central Europe. The objectives of this study were therefore to (i) evaluate the suitability of different regions in Germany for soybean production, (ii) investigate the magnitude and pattern of genotype-by-location interaction for soybean cultivated in Germany, and (iii) draw conclusions for soybean breeding in high latitude regions.

2. Results

This study analyses single-year data from seven locations with a focus on soybean performance in different regions in Germany. The locations in our study ( Figure 1 , Table 1 , Figure S1 ) were separated into two environmental target regions based on their latitude, as this represents an important factor of soybean adaptation. Furthermore, this approach separates locations in the south, where soybean is already established, from the locations at higher latitudes, where there is barely any soybean cultivation. Therefore, the environmental target region in the southern half of Germany (STH) comprised the locations Eckartsweier (EWE), Hohenheim (HOH), Moosburg an der Isar (MOS), Niedertraubling (NIT), and Landshut (LDH). The environmental target region of the northern half of Germany (NTH) was assigned to the locations Biendorf (BID) and Gülzow-Prüzen (GÜL).

Across all seven locations, the mean seed yield was 30.37 dt ha −1 ( Table 2 ) and across locations within the northern and southern half of Germany, the mean seed yield was 26.83 dt ha −1 and 31.82 dt ha −1 , respectively. Within single locations, mean seed yield ranged from 24.26 dt ha −1 for BID to 36.57 dt ha −1 for EWE. Protein content was high with an average of 45.09% across all locations. In single locations mean protein content ranged from 42.87% in MOS to 47.29% in EWE. Note that data on protein content was not available for GÜL. Seed yield was significantly ( p < 0.05) positively correlated between all locations except for GÜL and MOS ( Figure S2 ). For protein content, all locations were highly significantly ( p < 0.001) positively correlated ( Figure S3 ). Most trait correlations were not significant ( Figure S4 ). Exceptions were protein and oil content that were significantly ( p < 0.001) negatively correlated ( R = −0.91). Oil content was significantly ( p < 0.05) positively correlated with seed yield at a moderate level ( R = 0.31) and consequently protein content and seed yield were moderately negatively correlated ( R = −0.20), though this was not significant. In addition, protein content was significantly ( p < 0.05) positively correlated with days to maturity ( R = 0.43), indicating that later maturing genotypes tended to have a higher protein content.

The genotypic variance was significant for all traits across locations, within environmental target regions STH and NTH, as well as within single locations ( Table 2 ). In all cases was the genetic variance larger than the genotype-by-location interaction variance. The location variance was larger than the genotypic variance across all locations for all traits and within NTH for seed yield. Heritability across all locations was high to very high for all traits, ranging from 0.79 for plant height to 0.90 for seed yield and 0.94 for protein content ( Table 2 ). In addition, heritability was very high for seed yield in STH (0.91), but in NTH was substantially lower at 0.55. For seed yield, repeatability in single locations ranged from moderate in GÜL and BID (0.53 and 0.54) to very high in NIT and MOS (0.88 and 0.92). For protein content, repeatability within single locations ranged from moderate in BID (0.55) to very high in EWE and NIT (0.94).

In Figure 2 a,c, the scaled seed yield and protein content for each genotype within each location are shown. For example, genotype G35 showed a high seed yield relative to the other genotypes within all locations and G38 showed very high protein content relative to the other genotypes in all locations. Figure 2 b,d show boxplots of the best linear unbiased predictors (BLUP) within each location. The red lines indicate that for both seed yield and protein content, the genotypes performing best in at least one location had higher values than average in most locations. However, clear crossing-over interaction for seed yield, even for the highest-performing genotypes, could be observed ( Figure 2 b).

The angles between vectors in the Genotype plus Genotype-by-Environment (GGE) biplot for seed yield indicate a positive correlation between HOH, NIT, and LDH and between BID, GÜL, and EWE ( Figure 3 a). The performance vs. stability biplot ( Figure 3 b) indicates that G35 had the highest seed yield across locations while being mostly stable (relatively short red dotted line), whereas G39 was the lowest yielding genotype but had similar stability. Shukla’s stability variance ranged from 2.27 for G24 to 42.48 for G44 for seed yield ( Table S1 ). As seed yield and stability variance for seed yield were significantly ( p < 0.05) negatively correlated ( R = −0.37) ( Figure 4 a), rankings between both measures differed.

For protein content, the vector angles in the GGE biplot indicate positive correlations between all locations ( Figure 3 c). The respective performance vs. stability biplot ( Figure 3 d) shows G38 as the highest-performing but unstable genotype, and G47 as lowest performing while being relatively stable for protein content. The stability variance for protein content ranged from 0.03 for G43 to 6.01 for G50 ( Table S2 ). Protein content and stability variance for protein content were significantly ( p < 0.05) negatively correlated ( R = −0.32) ( Figure 4 b). As observed for the means for seed yield and protein content ( Figure S4 ), also their stability variances were not significantly correlated ( Figure 4 c). In NTH, seed yield ranged from 21.6 dt ha −1 to 30.47 dt ha −1 , whereas in STH seed yield was generally higher, ranging from 21.54 dt ha −1 to 38.77 dt ha −1 . Seed yield between NTH and STH was significantly ( p < 0.05) positively correlated at R = 0.59 ( Figure 5 c). Further, the top performing genotype (G35) was the same in both NTH and STH ( Figure 5 a,b, Table 3 ), but for the other genotypes, some moderate crossing-over can be observed in the respective boxplot ( Figure 5 b). In addition, stability rankings according to Shukla’s stability variance were different between both environmental target regions. Stability variance for NTH ranged from 0 for G3, G8, G15, G24, G25, and G46 as the most stable genotypes to 83.21 for G44 as the least stable genotype. For STH, stability variance ranged from 0.40 for G40 as the most stable and 38.57 for G44 as the least stable genotype ( Table 3 ). However, stability analysis for five genotypes (G5, G7, G16, G21, G42, and G50) was not possible due to missing observations in GÜL. GGE biplots, ranking genotypes based on both their performance and stability within both environmental target regions compared to a computed fictional ideal genotype, show differences in the most suitable genotypes for NTH and STH, except for G35, which emerged as most suitable for both ( Figure 5 d,e).

3. Discussion

3.1. Phenotypic Performance and Variance Components

Very high heritabilities across locations for both seed yield (0.90) and protein content (0.94) in conjunction with a lower, however highly significant ( p < 0.001), genotype-by-location interaction variance component ( σ G L 2 ), illustrate a moderate level of genotype-by-location interaction ( Table 2 ). For both traits, genotype-by-location interaction of the crossing-over type was observed, indicating the benefit of forming environmental target regions as subgroups of locations ( Figure 2 b,d). Importantly, to evaluate the repeatability of the observed interaction patterns, further evaluation across several years is required. The heritability for seed yield in STH was very high, but for NTH was substantially lower at 0.55. Accordingly, the two northern locations GÜL and BID forming the environmental target region NTH both showed lower repeatabilities (0.54 and 0.53), indicating a lower data quality in comparison to the other locations. This could be due to more inhomogeneous field conditions that could not sufficiently be accounted for by the field design or adaptation issues of at least some of the genotypes.

3.2. Performance and Trait Stability of Genotypes

In practical breeding, combining high performance and trait stability is a key objective to ensure reliable yields across different locations and years. In the case of this study, evaluating single-year data, yield stability via Shukla’s stability variance is analysed by estimating a genotype-specific σ G L 2 . Genotypes with stability variance σ i 2 = 0 do not show any genotype-by-location interaction and are therefore considered perfectly stable across locations.

The top-performing genotypes for seed yield differed between locations showing a crossing-over between most locations ( Figure 2 b). In accordance with this, stability rankings were quite different from the performance rankings of the genotypes for seed yield across all locations ( Table S1 ). However, the correlation between stability variance and seed yield was significant ( p < 0.05) and negative at a moderate level ( R = −0.34) ( Figure 4 a). This indicates a trend of higher-yielding genotypes tending to show higher yield stability across locations. However, the highest-yielding genotype G35 only reached the rank of 28 out of 50 in stability. Genotype G37, by contrast, ranked third in mean seed yield and second in stability, resulting in G37 being the closest to an ideal genotype for seed yield in our study. It should be noted that the rankings by performance were not identical to the respective GGE biplot ( Figure 3 b) but showed minor rank changes. This is likely due to the different methods of the metan package [ 22 ] to estimate the mean performance of each genotype.

Similar to the results for seed yield, protein content and its stability variance were moderately negatively correlated at R = −0.31 ( Figure 4 b). Stability variances for both traits were not significantly correlated with each other ( Figure 4 c). For protein content, only two top-performing genotypes were found, G38 in all locations but MOS, where G42 had the highest protein content ( Figure 2 c,d). Stability analysis according to Shukla’s stability variance, however, classified G38 as highly unstable with the rank 49 out of 50. Following the concept of dynamic stability, stability variance measures how the genotype performance changes relative to the environmental mean (Equation (1)). Genotype G38 did not follow this pattern. In consequence, despite being the superior genotype in almost all locations and in addition to showing high performance in MOS, G38 is not classified as stable. Notably, genotypes that vary around the location means will be classified as unstable, but also a genotype that performs better than the mean at every location. So from a practical point of view, classifying G38 as unstable is misleading since the main interest of a breeder lies in the consistently high performance of a genotype. Accordingly, it should be noted that trait performance (especially yield) is the most important breeding objective while stability measures can only support selection decisions. Notably, if the number of environments in a multi-environment trial is large enough, the estimation of mean performance across environments will automatically penalise unstable genotypes that vary around the environment means [ 23 ]. It should be noted here, that this study focuses on the genotypic stability across locations in one test year and therefore cannot provide conclusions about the trait stability across years.

3.3. Expansion of Soybean Cultivation to Higher Latitudes Is Possible

The mean seed yield across locations within STH was 18.6% higher than across locations within NTH ( Table 2 ). However, the significant ( p < 0.05) positive correlation ( R = 0.59) between seed yield in NTH and STH ( Figure 5 c) suggests that genotypes with high performance in STH tended to also show high seed yield relative to the other genotypes in NTH. This is further highlighted and supported by the genotype mean comparison shown in Figure 5 b. Despite the rank changes due to crossing-over genotype-by-location interaction between NTH and STH, high-ranking genotypes in STH tended to also rank higher in NTH in many cases. Generally, crossing-over type genotype-by-environment interaction can be met by different selection decisions between target regions or mega environments [ 17 , 23 , 24 ]. However, within mega environments or target regions, genotype-by-environment interaction should be minimised. In our case, this goal was reached with the environmental target regions, as the ratio of genotype-by-location interaction to genotypic variance ( σ G L 2 : σ G 2 ) decreased from 0.40 across all locations to 0.28 within NTH and 0.22 within STH.

When evaluating genotypes in both environmental target regions for high seed yield and yield stability in the year 2021, it appears that different selection decisions should be taken for NTH and STH ( Figure 5 d,e). The highest-performing genotype (G35) was the same in both target regions. However, for example, G8 was just as suitable for STH as G35, while it merely showed seed yield close to the mean in NTH. Consequently, despite the single-year analysis, the results indicate that selection decisions made for one environmental target region cannot confidently be transferred directly to another. To further explore this, we compared the mean seed yield of the five highest-ranking genotypes across all locations to the mean of the five highest-ranking genotypes within each of the two regions. Selecting for specific adaptation led to a 0.40% yield increase in STH compared to selection across all locations. For NTH, by contrast, the yield increase was higher at 2.32%. This analysis must, however, be treated with caution, due to the higher number of locations in STH. This is also an issue when analysing mean seed yield and stability across all seven locations ( Figure 3 b). Rankings and stability depictions are similar to those observed for STH. As an environmental target region consisting of five locations compared to NTH with only two locations, STH is overrepresented in the analysis across all locations. Keeping this in mind, our results indicate that separate analysis of environmental target regions enables more accurate selection against genotype-by-location interaction. However, short environment vectors for GÜL and BID ( Figure 3 a), proportional to the smaller variation within those locations ( Figure 2 b), show that both locations were not discriminating. This means that GÜL and BID were less suitable as test locations for the specific set of genotypes in this study compared to the other five locations [ 25 , 26 ]. One possible explanation could be less favourable environmental conditions during the growth period ( Table 1 ), as especially in GÜL many genotypes did not fully mature before the field had to be harvested. However, we could not observe any single environmental factor explaining the substantially lower yields in NTH ( Table 1 ). In consequence, the lower yields likely resulted from a combination of environmental conditions, including latitude.

The cause of genotype-by-environment interaction and trait instability are different sensitivities to environmental factors. In soybean, yield changes are often due to differences in tolerances to cooler temperatures [ 27 , 28 , 29 ] and photoperiod [ 11 , 30 ]. The aim of the differentiation between NTH and STH as environmental target regions was to mainly capture the effect of different latitudes on seed yield. In the north of Germany, due to typically lower temperatures and longer day length, earlier maturing genotypes are required compared to southern Germany in order to maximise the yield potential [ 29 , 30 , 31 ]. This generally indicates the need for selecting different genotypes for different latitudinal regions, which warrants further research. The majority of the genotypes evaluated in this study stem from a breeding program developed in the south of Germany at the location of EWE. Consequently, they were previously selected for high performance in that region. A question of practical relevance is therefore, whether we can successfully select soybeans in the south for cultivation in the north of Germany. Despite the correlation between the two regions, we found that evaluations made in the south did not sufficiently reflect the performance at higher latitudes. Consequently, specific selection decisions made separately in each of those regions appear necessary in soybean breeding. Nevertheless, the trend highlighted in this study needs to be validated by further research, also taking the effect of different growing seasons into account. Further, in practice, selection within mega environments or environmental target regions is only appropriate if the response to selection can be substantially increased [ 24 ]. Alternatively, an in part shared selection, for example, for simple traits in early generations, may be recommended.

3.4. Conclusions for Soybean Breeding at Higher Latitudes

This study illustrates that soybean can be cultivated in different regions in Germany, covering a wide range of latitudes from south to north. Notably, however, the suitability of genotypes for cultivation in southern Germany and more northern regions differed in part substantially. Consequently, the trend observed for the year 2021 indicates that different selection in breeding programs specifically targeting different regions appears beneficial. Key environmental influences affecting performance and therefore also genotype-by-environment interaction and stability in soybean are latitude, temperature, and precipitation [ 27 , 28 , 29 , 30 , 31 , 32 ]. Notably, to reliably classify mega environments, knowledge about the stability of the genotype-by-location interaction patterns over several years is necessary [ 17 ]. Patterns that are not repeatable across years cannot be utilised [ 23 ] and are therefore of no use for specific cultivar selection. Therefore, further research evaluating multi-environment trials comprising several years is crucial for a deeper understanding of these trends. Consequently, in addition to further exploring the environmental target regions defined in this study, it would be interesting to investigate the patterns of genotype-by-location interaction more thoroughly by evaluating a larger number of locations in different regions covering a wide range of latitudes across several years. An unconventional but promising approach to do so is citizen science which allows for a dense coverage of locations across an entire country or even across countries [ 33 ]. This would also offer further insights into the significance of latitude and specific weather conditions in the formation of environmental target regions and mega environments relevant to soybean breeding at higher latitudes.

4. Materials and Methods

4.1. Plant Material and Field Trials

We evaluated a total of 50 genotypes with favourable tofu quality, including eight registered cultivars ( Table S3 ). The field trials were conducted in 2021 at seven locations in Germany: EWE, HOH, MOS, NIT, LDH, BID, and GÜL ( Figure 1 , Table 1 ). Further two trial locations in Beckum (BEC) and Kiel (KIL) could not provide data, as at these locations the soybeans did not reach maturity due to bird pests and unfavourable weather conditions throughout the year. Weather data for the locations MOS and NIT were obtained from Agrarmeteorologie Bayern [ 34 ]. Weather data for the locations in Baden-Württemberg (EWE and HOH) were obtained from Agrarmeteorologie Baden-Württemberg [ 35 ]. For LDH, BID, and GÜL weather data were obtained from the Climate Data Center [ 36 ].

The trials were designed in an alpha lattice design with two replications, ten blocks per replication, and five genotypes per block. The soybeans were sown in four rows in plots of 6 m length and 1.5 m width. Standard agricultural practices for soybean production were followed. Seed yield in dt ha −1 was measured on the combine harvester and corrected for seed moisture content. Protein and oil content in percent were measured on dried seed (9–10% moisture content) in the laboratory with near-infrared spectroscopy (NIRS) using a Polytec PSS 2120 diode spectrometer and the corresponding software PSS-HOP (Polytec GmbH, Waldenbronn, Germany). Days to maturity were assessed as days between sowing and maturity determined by dry pods rustling when touched (growth stage R8). Further, plant height in cm and kernel dry matter in percent were measured.

4.2. Phenotypic Analysis

BLUPs and variance components were calculated across all locations and across locations within environmental target regions with the model

where y ijkm describes the phenotypic observation, µ is the intercept, g i is the effect of genotype i , l j is the effect of location j , gl ij is the interaction between genotype i and location j , r jk is the effect of replication k within location j , b jkm is the effect of block m within replication k in location j , and ε ijkm is the corresponding residual error term. BLUPs and variance components within locations were calculated with the model

where y ijk describes the phenotypic observation, µ is the intercept, g i is the effect of genotype i , r j is the effect of replication j , b jk is the effect of block k within replication j , and ε ijk is the corresponding residual error term. For both models, the variance components were tested for significance using likelihood ratio tests. Heritability ( h ²) and repeatability ( w ²) were calculated as proposed by Cullis et al. [ 37 ] and described by Piepho and Möhring [ 38 ]:

where ϑ B L U P ¯ describes the mean variance of the difference between two BLUPs and σ g 2 is the genotypic variance. The sole difference between these values is the calculation across environments ( h ²) and the calculation within environments ( w ²). Standardised values were calculated following the min-max normalisation. Correlations were calculated as Pearson’s product-moment correlation between pairwise complete observations.

4.3. Genotype-by-Location Interaction and Stability Analysis

Environmental target regions were defined by their latitude as an important component of soybean adaptation. For this, locations were separated into two groups by their geographic location: GÜL and BID in the northern half of Germany (NTH) and EWE, HOH, MOS, NIT, and LDH with rather similar latitudes in the south of Germany (STH).

Genotype plus Genotype-by-Environment (GGE) biplot is an approach proposed by Yan et al. [ 39 ]. It aids the interpretation of genotype-by-environment interaction by visualising genotype-by-environment interaction patterns in the form of a principal component analysis (PCA) biplot [ 26 , 39 ]. GGE biplot can be used to investigate the relationships between environments and the differences between genotypes across environments [ 25 , 26 ]. Further, different GGE biplots can be used to visualise the stability of genotypes and aid in selection decisions by combining performance and stability parameters within a single plot.

Genotype-by-location interaction was analysed using GGE biplot analysis [ 25 , 26 ] based on the equation

where y ij is the mean performance of genotype i in location j , λ 1 , and λ 2 are the singular values, ξ i 1 , and ξ i 2 are the eigenvectors of genotype i , and η 1 j and η 2 j are the eigenvectors of location j for the two largest principal components (PC1 and PC2), respectively, and ε ij is the corresponding error term. All plots were made without scaling and using environment-centering. To evaluate the relationships between locations, singular values were partitioned into the locations (SVP = 2). In order to compare genotypes, singular values were partitioned completely into the genotypes (SVP = 1).

Dynamic stability of the genotypes across locations was measured by estimating Shukla’s stability variance ( σ i 2 ) [ 40 ] based on the equation

where G is the number of genotypes, E is the number of locations, and y ij is the mean performance of genotype i in location j . Negative estimated values for σ i 2 were set to 0 [ 40 ].

All analyses were conducted with R within the RStudio [ 41 ] environment. BLUPs and variance components were calculated using the package ASReml-R [ 42 ]. GGE biplots, genotype-by-location heatmaps, and stability were calculated and produced using the package metan [ 22 ]. All graphs were drawn with the package ggplot2 [ 43 ]. Correlation matrices were produced with the extension package ggcorrplot2 [ 44 ].

Acknowledgments

We would like to thank the technical staff of the State Plant Breeding Institute for their efforts in the fieldwork and data collection. In addition, we would like to thank the staff at the additional test locations.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/plants12040756/s1 , Figure S1: Map of (a) locations and plots of (b) daily mean temperature and (c) precipitation summed up for each week from sowing to harvest for each location. Colours of locations in (a) indicate the locations in plots (b) and (c). Crosses in (a) show locations that could not be harvested; Figure S2: Correlation matrix for seed yield between seven locations; Figure S3: Correlation matrix for protein content between six locations; Figure S4: Correlation matrix for six traits across all locations; Table S1: Shukla’s stability measure and mean ranks of 50 genotypes for seed yield across seven locations; Table S2: Shukla’s stability measure and mean ranks of 50 genotypes for protein content across six locations; Table S3: List of genotypes evaluated in this study; data and data information.

Click here for additional data file.

Author Contributions

Conceptualization, C.A.D., V.H., W.L.L. and T.W.; methodology, C.A.D.; formal analysis, C.A.D.; investigation, C.A.D. and V.H.; data curation, C.A.D.; writing—original draft preparation, C.A.D.; writing—review and editing, C.A.D., V.H., W.L.L. and T.W.; visualization, C.A.D.; funding acquisition, V.H., W.L.L. and T.W. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data is contained in the Supplementary Materials .

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Funding Statement

This research was funded by the German Federal Ministry for Education and Research (BMBF) under grant number 031B1046A and handled by the project management organization Projektträger Jülich (PtJ). The responsibility for the contents lies with the authors. The funders were not involved in the design of the study, data collection, analysis, or preparation of the manuscript.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

References

1. Nielsen Z. Lebensmittel Zeitung. 2021. [(accessed on 19 March 2021)]. p. 41. Available online: www.lebensmittelzeitung.net .

2. Ahrens S. Umsatz mit Milchalternativen im Lebensmitteleinzelhandel nach Produktgruppe 2020. [(accessed on 22 August 2022)]. Available online: https://de.statista.com/statistik/daten/studie/1222972/umfrage/umsatz-mit-veganen-alternativen-im-lebensmitteleinzelhandel-nach-produktgruppe/

3. The Vegan Society A Growing Plant Milk Market. [(accessed on 22 August 2022)]. Available online: https://www.vegansociety.com/news/market-insights/plant-milk-market .

4. Qiu L.J., Chang R.Z. The origin and history of soybean. In: Guriqbal Singh G.S., editor. The Soybean: Botany, Production and Uses. CABI; Wallingford, UK: 2010. pp. 1–23.

5. Roßberg D., Recknagel J. Untersuchungen zur Anbaueignung von Sojabohnen in Deutschland. J. Cultiv. Plants. 2017;69:137–145. doi: 10.5073/JFK.2017.04.02.

6. Wilcox J.R. World Distribution and Trade of Soybean. In: Shibles R.M., Harper J.E., Wilson R.F., Shoemaker R.C., editors. Soybeans: Improvement, Production, and Uses. 3rd ed. American Socienty of Agronomy, Crop Science Society of America, and Soil Science Society of America; Madison, WI, USA: 2004. pp. 1–14.

7. Kumudini S. Soybean Growth and Development. In: Guriqbal Singh G.S., editor. The Soybean: Botany, Production and Uses. CABI; Wallingford, UK: 2010. pp. 48–73.

8. Panthee D.R. Varietal Improvement in Soybean. In: Guriqbal Singh G.S., editor. The Soybean: Botany, Production and Uses. CABI; Wallingford, UK: 2010. pp. 92–112.

9. Statistisches Bundesamt Landwirtschaftliche Betriebe, Fläche: Bundesländer, Jahre, Bodennutzungsarten. [(accessed on 1 July 2022)]. Available online: https://www-genesis.destatis.de/genesis/online?operation=abruftabelleBearbeiten&levelindex=2&levelid=1656687314161&auswahloperation=abruftabelleAuspraegungAuswaehlen&auswahlverzeichnis=ordnungsstruktur&auswahlziel=werteabruf&code=41271-0012&auswahltext=&wertauswahl=823&wertauswahl=275&werteabruf=Werteabruf#abreadcrumb .

10. Hamner K.C. Interrelation of Light and Darkness in Photoperiodic Induction. Bot. Gaz. 1940;101:658–687.

11. Zhang L.-X., Liu W., Tsegaw M., Xu X., Qi Y.-P., Sapey E., Liu L.-P., Wu T.-T., Sun S., Han T.-F. Principles and practices of the photo-thermal adaptability improvement in soybean. J. Integr. Agric. 2020;19:295–310. doi: 10.1016/S2095-3119(19)62850-9.

12. Food and Agricultural Organization of the United Nations FAOSTAT Statistical Database. [(accessed on 23 August 2022)]. Available online: https://www.fao.org/faostat/en/#data/TCL .

13. Fearnside P.M. Soybean cultivation as a threat to the environment in Brazil. Environ. Conserv. 2001;28:23–38. doi: 10.1017/S0376892901000030.

14. Bundesanstalt für Landwirtschaft und Ernährung . Ackerbohne, Erbse & Co.: Die Eiweißpflanzenstrategie des Bundesministeriums für Ernährung und Landwirtschaft zur Förderung des Leguminosenanbaus in Deutschland. Bundesministerium für Ernährung und Landwirtschaft; Berlin, Germany: 2020.

15. Lin C.S., Binns M.R., Lefkovitch L.P. Stability Analysis: Where Do We Stand? 1. Crop Sci. 1986;26:894–900. doi: 10.2135/cropsci1986.0011183X002600050012x.

16. Becker H.C., Léon J. Stability Analysis in Plant Breeding. Plant Breed. 1988;101:1–23. doi: 10.1111/j.1439-0523.1988.tb00261.x.

17. Annicchiarico P. Genotype x Environment Interaction: Challenges and Opportunities for Plant Breeding and Cultivar Recommendations. Food and Agriculture Organization of the United Nations; Rome, Italy: 2002.

18. Miedaner T. Grundlagen der Pflanzenzüchtung. DLG-Verlag; Frankfurt am Main, Germany: 2010.

19. Becker H. Pflanzenzüchtung, 2. Aufl. Ulmer; Stuttgart, Germany: 2011. UTB GmbH.

20. Matus-Cádiz M.A., Hucl P., Perron C.E., Tyler R.T. Genotype × Environment Interaction for Grain Color in Hard White Spring Wheat. Crop Sci. 2003;43:219. doi: 10.2135/cropsci2003.2190.

21. Fan X.-M., Kang M.S., Chen H., Zhang Y., Tan J., Xu C. Yield Stability of Maize Hybrids Evaluated in Multi-Environment Trials in Yunnan, China. Agron. J. 2007;99:220–228. doi: 10.2134/agronj2006.0144.

22. Olivoto T., L’ucio A.D. ’C. metan: An R package for multi-environment trial analysis. Methods Ecol. Evol. 2020;11:783–789. doi: 10.1111/2041-210X.13384.

23. Yan W. Analysis and Handling of G × E in a Practical Breeding Program. Crop Sci. 2016;56:2106–2118. doi: 10.2135/cropsci2015.06.0336.

24. Atlin G.N., Baker R.J., McRae K.B., Lu X. Selection Response in Subdivided Target Regions. Crop Sci. 2000;40:7–13. doi: 10.2135/cropsci2000.4017.

25. Yan W., Tinker N.A. Biplot analysis of multi-environment trial data: Principles and applications. Can. J. Plant Sci. 2006;86:623–645. doi: 10.4141/P05-169.

26. Yan W., Kang M.S. GGE Biplot Analysis: A Graphical Tool for Breeders, Geneticists, and Agronomists. CRC Press; Boca Raton, FL, USA: 2003.

27. Balko C., Hahn V., Ordon F. Kühletoleranz bei der Sojabohne (Glycine max (L.) Merr.)–Voraussetzung für die Ausweitung des Sojaanbaus in Deutschland. J. Für Kult. 2014;66:378–388. doi: 10.5073/JFK.2014.11.02.

28. Jähne F., Balko C., Hahn V., Würschum T., Leiser W.L. Cold stress tolerance of soybeans during flowering: QTL mapping and efficient selection strategies under controlled conditions. Plant Breed. 2019;138:708–720. doi: 10.1111/pbr.12734.

29. Karges K., Bellingrath-Kimura S.D., Watson C.A., Stoddard F.L., Halwani M., Reckling M. Agro-economic prospects for expanding soybean production beyond its current northerly limit in Europe. Eur. J. Agron. 2022;133:126415. doi: 10.1016/j.eja.2021.126415.

30. Kurasch A.K., Hahn V., Leiser W.L., Vollmann J., Schori A., Bétrix C.-A., Mayr B., Winkler J., Mechtler K., Aper J., et al. Identification of mega-environments in Europe and effect of allelic variation at maturity E loci on adaptation of European soybean. Plant Cell Environ. 2017;40:765–778. doi: 10.1111/pce.12896.

31. Sobko O., Stahl A., Hahn V., Zikeli S., Claupein W., Gruber S. Environmental Effects on Soybean (Glycine max (L.) Merr) Production in Central and South Germany. Agronomy. 2020;10:1847. doi: 10.3390/agronomy10121847.

32. Toleikiene M., Slepetys J., Sarunaite L., Lazauskas S., Deveikyte I., Kadziuliene Z. Soybean Development and Productivity in Response to Organic Management above the Northern Boundary of Soybean Distribution in Europe. Agronomy. 2021;11:214. doi: 10.3390/agronomy11020214.

33. Würschum T., Leiser W.L., Jähne F., Bachteler K., Miersch M., Hahn V. The soybean experiment ‘1000 Gardens’: A case study of citizen science for research, education, and beyond. Theor. Appl. Genet. 2019;132:617–626. doi: 10.1007/s00122-018-3134-2.

34. Agrarmeteorologie Bayern. [(accessed on 30 March 2022)]. Available online: www.wetter-by.de .

35. Agrarmeteorologie Baden-Württemberg. [(accessed on 30 March 2022)]. Available online: www.wetter-bw.de .

36. Climate Data Center. [(accessed on 30 March 2022)]. Available online: https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/daily/kl/historical/

37. Cullis B.R., Smith A.B., Coombes N.E. On the design of early generation variety trials with correlated data. JABES. 2006;11:381–393. doi: 10.1198/108571106X154443.

38. Piepho H.-P., Möhring J. Computing heritability and selection response from unbalanced plant breeding trials. Genetics. 2007;177:1881–1888. doi: 10.1534/genetics.107.074229.

39. Yan W., Hunt L.A., Sheng Q., Szlavnics Z. Cultivar Evaluation and Mega-Environment Investigation Based on the GGE Biplot. Crop Sci. 2000;40:597–605. doi: 10.2135/cropsci2000.403597x.

40. Shukla G.K. Some Statistical Aspects of Partitioning Genotype-Environmental Components of Variability. Heredity. 1972;29:237–345. doi: 10.1038/hdy.1972.87.

41. R Core Team . R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing; Vienna, Austria: 2021.

42. Butler D. Asreml: Fits the Linear Mixed Model. VSN International; Hemel Hempstead, UK: 2020.

43. Wickham H. Ggplot2: Elegant Graphics for Data Analysis. 2nd ed. Springer; New York, NY, USA: 2016.

44. Jun C., Granville M., Samson L.D. Ggcorrplot2: Visualize a Correlation Matrix using Ggplot2. R Package Version 0.1.2. [(accessed on 1 September 2022)]. Available online: https://github.com/caijun/ggcorrplot2 .

Associated Data

Supplementary Materials

Click here for additional data file.

Data Availability Statement

Data is contained in the Supplementary Materials .