# Envirotyping within a multi-environment trial allowed identifying genetic determinants of winter oilseed rape yield stability

Erwan Corlouer, Christopher Sauvage, Magalie Leveugle, Nathalie Nesi, Anne Laperche

► To cite this version:

Erwan Corlouer, Christopher Sauvage, Magalie Leveugle, Nathalie Nesi, Anne Laperche. Envirotyping within a multi-environment trial allowed identifying genetic determinants of winter oilseed rape yield stability. TAG Theoretical and Applied Genetics, 2024, 137 (7), pp.164. 10.1007/s00122-024-04664-3. hal-04618173

HAL Id: hal-04618173

https://hal.science/hal-04618173

Submitted on 20 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Envirotyping within a multi-environment trial allowed identifying genetic determinants of winter oilseed rape yield stability

Erwan Corlouer<fnref n="1" /> ID • Christopher Sauvage<fnref n="2" /> ID • Magalie Leveugle<fnref n="3" /> • Nathalie Nesi<fnref n="1" /> ID • Anne Laperche<fnref n="1" /> ID

Received: 20 December 2023 / Accepted: 28 May 2024

© The Author(s) 2024

## Abstract

**Key message** A comprehensive environmental characterization allowed identifying stable and interactive QTL for seed yield: QA09 and QC09a were detected across environments; whereas QA07a was specifically detected on the most stressed environments.

**Abstract** A main challenge for rapeseed consists in maintaining seed yield while adapting to climate changes and contributing to environmental-friendly cropping systems. Breeding for cultivar adaptation is one of the keys to meet this challenge. Therefore, we propose to identify the genetic determinant of seed yield stability for winter oilseed rape using GWAS coupled with a multi-environmental trial and to interpret them in the light of environmental characteristics. Due to a comprehensive characterization of a multi-environmental trial using 79 indicators, four contrasting envirotypes were defined and used to identify interactive and stable seed yield QTL. A total of four QTLs were detected, among which, QA09 and QC09a, were stable (detected at the multi-environmental trial scale or for different envirotypes and environments); and one, QA07a, was specifically detected into the most stressed envirotype. The analysis of the molecular diversity at QA07a showed a lack of genetic diversity within modern lines compared to older cultivars bred before the selection for low glucosinolate content. The results were discussed in comparison with other studies and methods as well as in the context of breeding programs.

## Introduction

The main challenge for agriculture consists in ensuring food security while adapting to climate changes and contributing to eco-friendly farming systems. Both will be met through cultural practices adaptations and the development of new crop varieties (Lobell et al. 2008). Indeed, crops must face the effects of climate change characterized by an increase in temperature and CO₂ concentration but also by the increase of intense climatic phenomena (e.g., droughts, floods or frosts; Bell et al. 2004). In addition, agroecological practices have emerged as a way to produce more food in a sustainable manner by enhancing ecology-based practices (Wezel et al.

2014) and reducing the use of pesticide and chemical inputs. These practices include the use of natural biological control of pests, cover crops, intercropping or cultivar mixtures, in addition to new soil management practices (reduced tillage, ...) that increases the complexity of plant-environment interactions (Lamichhane et al. 2018).

In the face of these profound changes, crops need to perform under fluctuating agricultural and climatic conditions. The question of breeding new varieties adapted to a wide range of pedoclimatic and environmental conditions, as well as varieties specifically dedicated to targeted environmental conditions, has arisen. Broadly-adapted varieties correspond to those with stable yields under all environmental conditions, while specifically-adapted varieties can present contrasting yields, depending on the characteristics of the environment.

The wide *versus* specific adaptation can be represented by the seed yield variation observed for a dedicated genotype under a variety of environmental conditions (El-Soda et al. 2014). Yield variation differences exist among genotypes indicating the possibility to breed for stability; however, this remains complex. Indeed, the observed phenotype can be expressed as the sum of a genetic effect (G), an

Communicated by Matthew N Nelson.

<figure><img src="image_1.png" /></figure>

Anne Laperche
anne.laperche@institut-agro.fr

<footnote n="1">1 IGEPP, INRAE, Institut Agro, Université de Rennes,
35650 Le Rheu, France</footnote>

<footnote n="2">2 Syngenta SA France, 1228 Chemin de L'Hobit,
31790 Saint-Sauveur, France</footnote>

<footnote n="3">3 Limagrain, 28 Route d'Ennezat, 63720 Chappes, France</footnote>
# Envirotyping within a multi-environment trial allowed identifying genetic determinants of winter oilseed rape yield stability

Erwan Corlouer<fnref n="1" /> ID • Christopher Sauvage<fnref n="2" /> ID • Magalie Leveugle<fnref n="3" /> • Nathalie Nesi<fnref n="1" /> ID • Anne Laperche<fnref n="1" /> ID

Received: 20 December 2023 / Accepted: 28 May 2024

© The Author(s) 2024

## Abstract

**Key message** A comprehensive environmental characterization allowed identifying stable and interactive QTL for seed yield: QA09 and QC09a were detected across environments; whereas QA07a was specifically detected on the most stressed environments.

**Abstract** A main challenge for rapeseed consists in maintaining seed yield while adapting to climate changes and contributing to environmental-friendly cropping systems. Breeding for cultivar adaptation is one of the keys to meet this challenge. Therefore, we propose to identify the genetic determinant of seed yield stability for winter oilseed rape using GWAS coupled with a multi-environmental trial and to interpret them in the light of environmental characteristics. Due to a comprehensive characterization of a multi-environmental trial using 79 indicators, four contrasting envirotypes were defined and used to identify interactive and stable seed yield QTL. A total of four QTLs were detected, among which, QA09 and QC09a, were stable (detected at the multi-environmental trial scale or for different envirotypes and environments); and one, QA07a, was specifically detected into the most stressed envirotype. The analysis of the molecular diversity at QA07a showed a lack of genetic diversity within modern lines compared to older cultivars bred before the selection for low glucosinolate content. The results were discussed in comparison with other studies and methods as well as in the context of breeding programs.

## Introduction

The main challenge for agriculture consists in ensuring food security while adapting to climate changes and contributing to eco-friendly farming systems. Both will be met through cultural practices adaptations and the development of new crop varieties (Lobell et al. 2008). Indeed, crops must face the effects of climate change characterized by an increase in temperature and CO₂ concentration but also by the increase of intense climatic phenomena (e.g., droughts, floods or frosts; Bell et al. 2004). In addition, agroecological practices have emerged as a way to produce more food in a sustainable manner by enhancing ecology-based practices (Wezel et al.

2014) and reducing the use of pesticide and chemical inputs. These practices include the use of natural biological control of pests, cover crops, intercropping or cultivar mixtures, in addition to new soil management practices (reduced tillage, ...) that increases the complexity of plant-environment interactions (Lamichhane et al. 2018).

In the face of these profound changes, crops need to perform under fluctuating agricultural and climatic conditions. The question of breeding new varieties adapted to a wide range of pedoclimatic and environmental conditions, as well as varieties specifically dedicated to targeted environmental conditions, has arisen. Broadly-adapted varieties correspond to those with stable yields under all environmental conditions, while specifically-adapted varieties can present contrasting yields, depending on the characteristics of the environment.

The wide *versus* specific adaptation can be represented by the seed yield variation observed for a dedicated genotype under a variety of environmental conditions (El-Soda et al. 2014). Yield variation differences exist among genotypes indicating the possibility to breed for stability; however, this remains complex. Indeed, the observed phenotype can be expressed as the sum of a genetic effect (G), an

Communicated by Matthew N Nelson.

<figure><img src="image_2.png" /></figure>

Anne Laperche
anne.laperche@institut-agro.fr

<footnote n="1">1 IGEPP, INRAE, Institut Agro, Université de Rennes,
35650 Le Rheu, France</footnote>

<footnote n="2">2 Syngenta SA France, 1228 Chemin de L'Hobit,
31790 Saint-Sauveur, France</footnote>

<footnote n="3">3 Limagrain, 28 Route d'Ennezat, 63720 Chappes, France</footnote>
environmental effect (E) defined here as the combination of pedoclimatic conditions and cultural practices, and the G×E interaction, corresponding to the modification of the phenotypic response to environment between genotypes (Becker and Leon 1988; Malosetti et al. 2016; van Eeuwijk et al. 2016). Breeding programs must therefore consider SY stability to breed for new varieties better adapted to fluctuating environments (Cooper et al. 2020; Snowdon et al. 2021). Classically, breeding programs rely on average performance of genotypes grown under multi-environmental and multi-annual field trials. Although this allows passive selection for a small adaptive effect (Snowdon et al. 2021), this methodology does not allow a direct access to the genetic determinism of seed yield stability (Garin et al. 2020).

A first approach to unravel the genetic determinism of SY stability consists in a comparison of QTL detected using a population trialed across a multi-environmental network. This method is efficient to access to the QTL stability or specificity, but the results remain difficult to interpret (Garin et al. 2020). To gain power in detecting QTL associated to plant adaptation, stability indicators have been developed and used in genetic analyses of stability. For instance, ecovalence indices (Wricke 1962; dos Santos Silva et al. 2021) or AMMI stability values are used to describe genotypic contribution to G×E and to characterize the adaptability versus stability in different environments (Purchase et al. 2000; Bouchet et al. 2016; Lozada and Carter, 2020; Has-sani et al. 2023). Genotypic reaction norms to environmen-tal gradient are also commonly used; they correspond for instance to Finlay and Wilkinson's regression slope (1963) (see Diouf et al. 2020; Xavier et al. 2018; Mangin et al. 2017 for application examples). More recently, linear mixed models and factor analytic models have been developed to consider simultaneously all environments and genotypes in a multi-environment trial (MET) (van Eeuwijk et al. 2010, 2016; Smith et al 2021) and have opened new avenues to detect G×E interactions and access its genetic determinism (Malosetti et al. 2008; Happ et al. 2021; Chidzanga et al 2022; Chaves et al. 2023). However, without a comprehen-sive environmental characterization of the MET, the QTL involved in stability is difficult to interpret in terms of the underlying mechanisms of plant adaptation.

Rapeseed is a major worldwide crop cultivated mainly for its seeds oil and meal production, presenting an estimated annual production of 71.3 Mt in 2021 (FAOSTAT 2023). Farmers have more and more difficulties to maintain rape-seed seed yield under adverse environmental conditions (e.g., insect damages during fall, nutritional constraints and abiotic stresses). This leads to a reduction of the cultivated surfaces (– 21% between 2016 and 2022, AGRESTE 2023). Thus, improving seed yield stability is of major concern for oilseed rape (Zandberg et al 2022; Cowling et al 2023). Seed yield (SY) is a complex trait defined by multiple components as

plant population density, the number of pods per plant, the number of seeds per pods or the seed weight (Diepenbrock 2000). The potential seed yield of winter oilseed rape is determined since the end of the fall but depends on many factors and stresses occurring all along the crop life cycle as abiotic stresses (temperature, water and radiation), biotic stresses (pest, pathogens and weeds) or nutritional stresses, particularly with nitrogen (Rathke et al. 2006). All these constraints on winter oilseed rape seed yield result in an important environmental effect and G×E interaction that explained around 10% of the seed yield variation under French conditions (Bouchet et al. 2016; Corlouer et al. 2019). Numerous studies reported the genetic determinant of SY in rapeseed (reviewed by Delourme et al. 2018) and reported a high number of QTL (Shi et al. 2009; Raboanatahiry et al. 2018), thus confirming SY as highly polygenic. The comparison between QTL detected in different environments revealed some QTL being characterized as “stable” (i.e., detected across all environments) and QTL being characterized as “interactive” (i.e., specific to an environmental condition, Bouchet et al. 2016). However, the QTL specificity observed in a MET is rarely associated with the identification of the environmental features causing the observed adaptation.

In this context, this study aims at identifying the genetic determinants of seed yield stability in winter oilseed rape and at interpreting them in the light of environmental characteristics. We based our strategy on the analysis of winter oilseed rape accessions experimented across a multi-environment trial (MET-47) consisting in 47 environments representing the diversity of French growing conditions. First, we developed a comprehensive characterization of the MET-47 to identify the limiting factors that occurred. Then, to identify the genetic determinant of SY stability, a panel of 173 accessions was experimented in a sub-MET of 22 environments out of the 47 (MET-22). Envirotypes were defined among the MET-22 and corresponded to the clustering of the 22 environments according to their limiting factors pattern. Then, GWAS analyses were carried out using BLUE (Best Linear Unbiased Estimator) obtained for each genotype and each envirotype to identify QTL specific to environmental patterns. Finally, a genetic diversity analysis was carried out to decipher the potential impact of breeding for seed quality on the reduction of genetic diversity at those detected QTL that may limit adaptation to specific environmental conditions.

# Materials and methods

## Field network description

Field experiments were run in a multi-environment trial (MET-47) consisting of 47 environments defined as
environmental effect (E) defined here as the combination of pedoclimatic conditions and cultural practices, and the G×E interaction, corresponding to the modification of the phenotypic response to environment between genotypes (Becker and Leon 1988; Malosetti et al. 2016; van Eeuwijk et al. 2016). Breeding programs must therefore consider SY stability to breed for new varieties better adapted to fluctuating environments (Cooper et al. 2020; Snowdon et al. 2021). Classically, breeding programs rely on average performance of genotypes grown under multi-environmental and multi-annual field trials. Although this allows passive selection for a small adaptive effect (Snowdon et al. 2021), this methodology does not allow a direct access to the genetic determinism of seed yield stability (Garin et al. 2020).

A first approach to unravel the genetic determinism of SY stability consists in a comparison of QTL detected using a population trialed across a multi-environmental network. This method is efficient to access to the QTL stability or specificity, but the results remain difficult to interpret (Garin et al. 2020). To gain power in detecting QTL associated to plant adaptation, stability indicators have been developed and used in genetic analyses of stability. For instance, ecovalence indices (Wricke 1962; dos Santos Silva et al. 2021) or AMMI stability values are used to describe genotypic contribution to G×E and to characterize the adaptability versus stability in different environments (Purchase et al. 2000; Bouchet et al. 2016; Lozada and Carter, 2020; Has-sani et al. 2023). Genotypic reaction norms to environmen-tal gradient are also commonly used; they correspond for instance to Finlay and Wilkinson's regression slope (1963) (see Diouf et al. 2020; Xavier et al. 2018; Mangin et al. 2017 for application examples). More recently, linear mixed models and factor analytic models have been developed to consider simultaneously all environments and genotypes in a multi-environment trial (MET) (van Eeuwijk et al. 2010, 2016; Smith et al 2021) and have opened new avenues to detect G×E interactions and access its genetic determinism (Malosetti et al. 2008; Happ et al. 2021; Chidzanga et al 2022; Chaves et al. 2023). However, without a comprehen-sive environmental characterization of the MET, the QTL involved in stability is difficult to interpret in terms of the underlying mechanisms of plant adaptation.

Rapeseed is a major worldwide crop cultivated mainly for its seeds oil and meal production, presenting an estimated annual production of 71.3 Mt in 2021 (FAOSTAT 2023). Farmers have more and more difficulties to maintain rape-seed seed yield under adverse environmental conditions (e.g., insect damages during fall, nutritional constraints and abiotic stresses). This leads to a reduction of the cultivated surfaces (– 21% between 2016 and 2022, AGRESTE 2023). Thus, improving seed yield stability is of major concern for oilseed rape (Zandberg et al 2022; Cowling et al 2023). Seed yield (SY) is a complex trait defined by multiple components as

plant population density, the number of pods per plant, the number of seeds per pods or the seed weight (Diepenbrock 2000). The potential seed yield of winter oilseed rape is determined since the end of the fall but depends on many factors and stresses occurring all along the crop life cycle as abiotic stresses (temperature, water and radiation), biotic stresses (pest, pathogens and weeds) or nutritional stresses, particularly with nitrogen (Rathke et al. 2006). All these constraints on winter oilseed rape seed yield result in an important environmental effect and G×E interaction that explained around 10% of the seed yield variation under French conditions (Bouchet et al. 2016; Corlouer et al. 2019). Numerous studies reported the genetic determinant of SY in rapeseed (reviewed by Delourme et al. 2018) and reported a high number of QTL (Shi et al. 2009; Raboanatahiry et al. 2018), thus confirming SY as highly polygenic. The comparison between QTL detected in different environments revealed some QTL being characterized as “stable” (i.e., detected across all environments) and QTL being characterized as “interactive” (i.e., specific to an environmental condition, Bouchet et al. 2016). However, the QTL specificity observed in a MET is rarely associated with the identification of the environmental features causing the observed adaptation.

In this context, this study aims at identifying the genetic determinants of seed yield stability in winter oilseed rape and at interpreting them in the light of environmental characteristics. We based our strategy on the analysis of winter oilseed rape accessions experimented across a multi-environment trial (MET-47) consisting in 47 environments representing the diversity of French growing conditions. First, we developed a comprehensive characterization of the MET-47 to identify the limiting factors that occurred. Then, to identify the genetic determinant of SY stability, a panel of 173 accessions was experimented in a sub-MET of 22 environments out of the 47 (MET-22). Envirotypes were defined among the MET-22 and corresponded to the clustering of the 22 environments according to their limiting factors pattern. Then, GWAS analyses were carried out using BLUE (Best Linear Unbiased Estimator) obtained for each genotype and each envirotype to identify QTL specific to environmental patterns. Finally, a genetic diversity analysis was carried out to decipher the potential impact of breeding for seed quality on the reduction of genetic diversity at those detected QTL that may limit adaptation to specific environmental conditions.

# Materials and methods

## Field network description

Field experiments were run in a multi-environment trial (MET-47) consisting of 47 environments defined as

<figure><img src="image_3.png" /></figure>
combinations of ‘year × location × nitrogen (N) fertilization’ in France between 2011 and 2018 (Supplemental Data 1). Each individual trial was conducted using classical crop management for winter oilseed rape (WOSR) with comprehensive protection against weeds, pests and pathogens. Optimal N fertilization was estimated using the balance sheet method (Rémy and Hébert 1977; Parnaudeau et al. 2009) for a target yield of 35 q ha⁻¹ and applied in a subset of 19 environments defined as high N (N⁺). In contrast, a low N fertilization regime was applied in the remaining 28 environments defined as low N (N⁻) that corresponded to the N⁺ regime lowered by 80–100 kg ha⁻¹ of N (Supplemental Data 1). Each environment was designed as a randomized complete block design with two to four replicates depending on the environment, with an individual plot area ranging from 6.75 to 14 m². Mature dry seeds were harvested when the vegetative parts were fully senescent and the seeds were dark and hard. The targeted traits were the seed yield (SY in q ha⁻¹) determined for each genotype in each trial and adjusted to 0% water content and 0% impurities, as well as the seed number (SN in seeds m⁻²) calculated according to the SY and the thousand seed weight (SN = SY × 100 000/TSW).

## Plant material

The WOSR genotypes “Aviso” and “Montego” were trialed over the whole MET-47 and therefore considered as probe genotypes for environmental characterization as reported by Corlouer et al. (2019). A diversity panel of 173 WOSR accessions (hereafter referred to as P173) was scored for SY and SN over a subset of 22 environments (11 N⁺ and 11 N⁻, called hereafter MET-22) run over the 2013–2014 and 2014–2015 growing seasons. The P173 accessions originated primarily from Western Europe and mostly represented commercial varieties released from 1959 to 2010 (Supplemental Data 2) with 31 accessions of the ‘++’ type (high contents in both erucic acid and glucosinolates; high C22:1, high GSL), 15 accessions of the ‘0+’ type (low C22:1, high GSL), 1 accession of the ‘+0’ type (high C22:1, low GSL) and 126 accessions of the ‘00’ type (low C22:1, low GSL) of which 46 were elite lines. Seeds were provided by the BrACySol Biological Resource Center or private seed companies. All accessions of the P173 population were produced in the same nursery at Le Rheu, France in 2010–2011 growth season to avoid seed lot biases. The P173 population was genotyped using the Brassica 60K Illumina Infinium™ array (Clarke et al. 2016) and an exome sequence capture assay (Leveugle et al. 2015). A total of 217,805 SNP were scored and validated for the current study using a threshold of 2.5% for the minor allele frequency (MAF) and of 10% for the missing values. The missing genotyping data were inferred using Beagle v3 software (Browning and Browning

2009). All the SNPs were physically anchored onto the latest *Brassica napus* reference genome of Darmor-bzh (Rousseau-Gueutin et al. 2020).

## Environmental characterization and envirotyping

Pedoclimatic indicators were calculated as described by Corlouer et al. (2019) and estimated for each of the key periods of the winter oilseed rape growing cycle. Briefly, these periods were fixed on meteorological data and the phenology of the probe genotypes. They covered fall (F), climatic winter (CW), bolting (B), seed number fixation (P300), reserve allocation to the pod (P600) and seed growth (P1000) stages. The 70 indicators correspond to four main categories: water stress, temperature, radiation and vernalization conditions (Corlouer et al. 2019). In addition, new indicators were developed to consider the contrasting N nutrition status (Supplemental Data 3). These included indicators related to the plant N status such as the nitrogen nutrition index (NNI; Colnenne et al. 1998) as well as the aerial dry biomass (ADB) and the aerial nitrogen quantity (AN), all scored at bolting stage (BBCH50, Lancashire et al. 1991) with a minimum delay of 2 weeks after the latest N supply. These values were measured on plants collected in the field on a surface of 1 m² and are expressed in quantity per hectare. On the other hand, indicators related to the N fertilization management were also considered such as the total amount of N supplied to the crop (Ntotal, in kg ha⁻¹), and the water regime at the time of N supply was assessed considering the maximum number of days without rainfall during the 10 days preceding the N supply (DBI, dryness before input), the number of days without rainfall over 10 days after the N supply (DAI, dryness after input), the sum of the rainfalls over 10 days after the N supply (RF, rainfall), the number of days of leaching over 10 days after the N supply (RO, run off) and the mean value of the water soil content at the time of the N supply (from 3 days before up to 10 days after the input) expressed in percentage of the maximal water soil content (WSC_I). Given that N fertilization was provided in one, two or three applications depending on the local pedoclimatic conditions, we defined a unique indicator for each trial that considers the number of N supplies, calculated as following:

$$ \text{Final indicator} = \sum_{i=0}^{n} \text{Indicator}_i \times \frac{N_i}{N_{\text{total}}} \quad (1) $$

where Indicatorᵢ is the given indicator (DBI, DAI, RF, RO or WSC_I) calculated at the N supply i, Nᵢ is the amount of N brought to the environment at the application i and Ntotal the total N amount considering all the N applications. When no N fertilization was applied, the final indicators were fixed to 0. Missing indicators were imputed using the missMDA package (Josse and Husson 2016).
combinations of ‘year × location × nitrogen (N) fertilization’ in France between 2011 and 2018 (Supplemental Data 1). Each individual trial was conducted using classical crop management for winter oilseed rape (WOSR) with comprehensive protection against weeds, pests and pathogens. Optimal N fertilization was estimated using the balance sheet method (Rémy and Hébert 1977; Parnaudeau et al. 2009) for a target yield of 35 q ha⁻¹ and applied in a subset of 19 environments defined as high N (N⁺). In contrast, a low N fertilization regime was applied in the remaining 28 environments defined as low N (N⁻) that corresponded to the N⁺ regime lowered by 80–100 kg ha⁻¹ of N (Supplemental Data 1). Each environment was designed as a randomized complete block design with two to four replicates depending on the environment, with an individual plot area ranging from 6.75 to 14 m². Mature dry seeds were harvested when the vegetative parts were fully senescent and the seeds were dark and hard. The targeted traits were the seed yield (SY in q ha⁻¹) determined for each genotype in each trial and adjusted to 0% water content and 0% impurities, as well as the seed number (SN in seeds m⁻²) calculated according to the SY and the thousand seed weight (SN = SY × 100 000/TSW).

## Plant material

The WOSR genotypes “Aviso” and “Montego” were trialed over the whole MET-47 and therefore considered as probe genotypes for environmental characterization as reported by Corlouer et al. (2019). A diversity panel of 173 WOSR accessions (hereafter referred to as P173) was scored for SY and SN over a subset of 22 environments (11 N⁺ and 11 N⁻, called hereafter MET-22) run over the 2013–2014 and 2014–2015 growing seasons. The P173 accessions originated primarily from Western Europe and mostly represented commercial varieties released from 1959 to 2010 (Supplemental Data 2) with 31 accessions of the ‘++’ type (high contents in both erucic acid and glucosinolates; high C22:1, high GSL), 15 accessions of the ‘0+’ type (low C22:1, high GSL), 1 accession of the ‘+0’ type (high C22:1, low GSL) and 126 accessions of the ‘00’ type (low C22:1, low GSL) of which 46 were elite lines. Seeds were provided by the BrACySol Biological Resource Center or private seed companies. All accessions of the P173 population were produced in the same nursery at Le Rheu, France in 2010–2011 growth season to avoid seed lot biases. The P173 population was genotyped using the Brassica 60K Illumina Infinium™ array (Clarke et al. 2016) and an exome sequence capture assay (Leveugle et al. 2015). A total of 217,805 SNP were scored and validated for the current study using a threshold of 2.5% for the minor allele frequency (MAF) and of 10% for the missing values. The missing genotyping data were inferred using Beagle v3 software (Browning and Browning

2009). All the SNPs were physically anchored onto the latest *Brassica napus* reference genome of Darmor-bzh (Rousseau-Gueutin et al. 2020).

## Environmental characterization and envirotyping

Pedoclimatic indicators were calculated as described by Corlouer et al. (2019) and estimated for each of the key periods of the winter oilseed rape growing cycle. Briefly, these periods were fixed on meteorological data and the phenology of the probe genotypes. They covered fall (F), climatic winter (CW), bolting (B), seed number fixation (P300), reserve allocation to the pod (P600) and seed growth (P1000) stages. The 70 indicators correspond to four main categories: water stress, temperature, radiation and vernalization conditions (Corlouer et al. 2019). In addition, new indicators were developed to consider the contrasting N nutrition status (Supplemental Data 3). These included indicators related to the plant N status such as the nitrogen nutrition index (NNI; Colnenne et al. 1998) as well as the aerial dry biomass (ADB) and the aerial nitrogen quantity (AN), all scored at bolting stage (BBCH50, Lancashire et al. 1991) with a minimum delay of 2 weeks after the latest N supply. These values were measured on plants collected in the field on a surface of 1 m² and are expressed in quantity per hectare. On the other hand, indicators related to the N fertilization management were also considered such as the total amount of N supplied to the crop (Ntotal, in kg ha⁻¹), and the water regime at the time of N supply was assessed considering the maximum number of days without rainfall during the 10 days preceding the N supply (DBI, dryness before input), the number of days without rainfall over 10 days after the N supply (DAI, dryness after input), the sum of the rainfalls over 10 days after the N supply (RF, rainfall), the number of days of leaching over 10 days after the N supply (RO, run off) and the mean value of the water soil content at the time of the N supply (from 3 days before up to 10 days after the input) expressed in percentage of the maximal water soil content (WSC_I). Given that N fertilization was provided in one, two or three applications depending on the local pedoclimatic conditions, we defined a unique indicator for each trial that considers the number of N supplies, calculated as following:

$$ \text{Final indicator} = \sum_{i=0}^{n} \text{Indicator}_i \times \frac{N_i}{N_{\text{total}}} \quad (1) $$

where Indicatorᵢ is the given indicator (DBI, DAI, RF, RO or WSC_I) calculated at the N supply i, Nᵢ is the amount of N brought to the environment at the application i and Ntotal the total N amount considering all the N applications. When no N fertilization was applied, the final indicators were fixed to 0. Missing indicators were imputed using the missMDA package (Josse and Husson 2016).
A total of 79 indicators (Supplemental Data 3) were considered in the present study. For a more precise analysis of the SY stability, an envirotyping was carried out according to the methodology proposed in Corlouer et al. (2019). Briefly, a univariate Partial Least Square (PLS) regression analysis was run at the MET-47 scale by regressing the mean seed yield of Aviso and Montego for each of the 47 environments on the 79 environmental indicators. A variable selection was performed to keep the reduced set of indicators that best explained seed yield variation, hereafter referred as limiting factors. Then, a hierarchical clustering was performed on the environments of the MET-22 to define envirotypes according to their pattern of limiting factors.

## Phenotypic data analyses

Trait analysis was run using three different mixed linear models that were fitted using the “lme4” (Bates et al. 2015) and ‘lmerTest’ packages (Kuznetsova et al. 2017), and corresponded to the three scales of analysis used in this study, i.e., the whole MET-22 scale, the envirotype scale (cluster of different environments presenting the same pattern of limiting factors) or the single environment scale. In the different models, the random effects are underlined in the equation, the others effects are considered as fixed.

A first linear mixed model (2) was fitted at the MET-22 scale including an effect of the envirotype to evaluate the impact of the envirotyping on the variance repartition.

$$Y_{ijkl} = \mu + G_i + En_l + En_{l(k)} + G_i \times En_l + G_i \times En_{l(k)} + \frac{R_{j(l \times k)}}{\mu} + \frac{\varepsilon_{ijkl}}{\mu} \quad (2)$$

where $Y_{ijkl}$ is the phenotypic value, $\mu$ is the population mean, $G_i$ stands for the effect of genotype i, $En_l$ for the envirotype l, $En_{l(k)}$ for the environment k nested in the envirotype l, $R_{j(l \times k)}$ for the replicate j, nested in the environment k and envirotype l; and $\varepsilon_{ijkl}$ is the residual.

Then, at the single environment scale, models (3) and (4) were used, and at the envirotype and at the MET-22 scale, models (5) and (6) were used to study the genotypic effect, to calculate the corresponding best linear unbiased estimators (BLUE), and to estimate trait heritability. These models are presented below.

Model (3) was fitted at the scale of each environment:

$$Y_{ij} = \mu + G_i + R_j + \varepsilon_{ij} \quad (3)$$

where $Y_{ij}$ is the phenotypic value, $\mu$ is the population mean, $G_i$ stands for the effect of genotype i, $R_j$ for the replicate j and $\varepsilon_{ij}$ is the residual. The replicate effect and the residual were declared as random. The heritability was calculated as:

$$h^2 = \frac{\sigma_G^2}{\sigma_G^2 + \frac{\sigma_r^2}{r}} \quad (4)$$

where $\sigma_G^2$ is the genetic variance, $\sigma_r^2$ the residual variance and $r$ the number of replicates per genotype.

Model (5) was applied at the envirotype or at the MET-22 scale:

$$Y_{ijk} = \mu + G_i + E_k + G_i \times E_k + \frac{R_{j(k)}}{\mu} + \frac{\varepsilon_{ijk}}{\mu} \quad (5)$$

where $Y_{ijk}$ is the phenotypic value, $\mu$ is the population mean, $G_i$ stands for the effect of genotype i, $E_k$ for the environment k, $R_{j(k)}$ for the replicate j, nested in the environment k and $\varepsilon_{ijk}$ is the residual. The corresponding heritability was defined as follow with all terms declared as random in Eq. (5):

$$h^2 = \frac{\sigma_G^2}{\sigma_G^2 + \frac{\sigma_{G \times E}^2}{t} + \frac{\sigma_r^2}{r \times e}} \quad (6)$$

where $\sigma_G^2$ is the genetic variance, $\sigma_{G \times E}^2$ the $G \times E$ interaction variance, $\sigma_r^2$ the residual variance, $e$ the number of environment and $r$ the number of replicates per genotype.

## GWAS analyses

The BLUE defined at each scale (environment, envirotype, MET-22) using models (3) and (5) for SY and SN was used to perform the genetic analyses. GWAS analyses were conducted using the FastLMM algorithm (Lippert et al. 2011) using a kinship matrix calculated for each linkage group as described by Rincent et al. (2014) and following the Astle algorithm (Astle and Balding 2009). A Bonferroni detection threshold was set to 10%, based on a corrected SNP population as proposed by the simpleM method developed by Gao et al. (2008). This method was based on the composite linkage disequilibrium (CLD) correlation between SNP. The CLD was used to calculate the effective number of independent tests ($M_{eff}$). In this study, the $M_{eff}$ was calculated and fixed to 19,886 SNP. We defined the threshold of highly significant associated SNP as:

$$t = - \log_{10} \left( \frac{\alpha_G}{M_{eff}} \right) \quad (7)$$

where $t$ is the value of the value of the threshold, $\alpha_G$ is the threshold of the $p$-value and was fixed to 0.1 and $M_{eff}$ was the effective population of SNP. The final threshold used was $t=5.30$.

## Genetic diversity analysis

A genetic diversity was conducted at the associated loci detected through GWAS. This evaluated the potential impact
A total of 79 indicators (Supplemental Data 3) were considered in the present study. For a more precise analysis of the SY stability, an envirotyping was carried out according to the methodology proposed in Corlouer et al. (2019). Briefly, a univariate Partial Least Square (PLS) regression analysis was run at the MET-47 scale by regressing the mean seed yield of Aviso and Montego for each of the 47 environments on the 79 environmental indicators. A variable selection was performed to keep the reduced set of indicators that best explained seed yield variation, hereafter referred as limiting factors. Then, a hierarchical clustering was performed on the environments of the MET-22 to define envirotypes according to their pattern of limiting factors.

## Phenotypic data analyses

Trait analysis was run using three different mixed linear models that were fitted using the “lme4” (Bates et al. 2015) and ‘lmerTest’ packages (Kuznetsova et al. 2017), and corresponded to the three scales of analysis used in this study, i.e., the whole MET-22 scale, the envirotype scale (cluster of different environments presenting the same pattern of limiting factors) or the single environment scale. In the different models, the random effects are underlined in the equation, the others effects are considered as fixed.

A first linear mixed model (2) was fitted at the MET-22 scale including an effect of the envirotype to evaluate the impact of the envirotyping on the variance repartition.

$$Y_{ijkl} = \mu + G_i + En_l + En_{l(k)} + G_i \times En_l + G_i \times En_{l(k)} + \frac{R_{j(l \times k)}}{\sum R_{j(l \times k)}} + \epsilon_{ijkl} \quad (2)$$

where $Y_{ijkl}$ is the phenotypic value, $\mu$ is the population mean, $G_i$ stands for the effect of genotype i, $En_l$ for the envirotype l, $En_{l(k)}$ for the environment k nested in the envirotype l, $R_{j(l \times k)}$ for the replicate j, nested in the environment k and envirotype l; and $\epsilon_{ijkl}$ is the residual.

Then, at the single environment scale, models (3) and (4) were used, and at the envirotype and at the MET-22 scale, models (5) and (6) were used to study the genotypic effect, to calculate the corresponding best linear unbiased estimators (BLUE), and to estimate trait heritability. These models are presented below.

Model (3) was fitted at the scale of each environment:

$$Y_{ij} = \mu + G_i + R_j + \epsilon_{ij} \quad (3)$$

where $Y_{ij}$ is the phenotypic value, $\mu$ is the population mean, $G_i$ stands for the effect of genotype i, $R_j$ for the replicate j and $\epsilon_{ij}$ is the residual. The replicate effect and the residual were declared as random. The heritability was calculated as:

$$h^2 = \frac{\sigma_G^2}{\sigma_G^2 + \frac{\sigma_r^2}{r}} \quad (4)$$

where $\sigma_G^2$ is the genetic variance, $\sigma_r^2$ the residual variance and $r$ the number of replicates per genotype.

Model (5) was applied at the envirotype or at the MET-22 scale:

$$Y_{ijk} = \mu + G_i + E_k + G_i \times E_k + \frac{R_{j(k)}}{\sum R_{j(k)}} + \epsilon_{ijk} \quad (5)$$

where $Y_{ijk}$ is the phenotypic value, $\mu$ is the population mean, $G_i$ stands for the effect of genotype i, $E_k$ for the environment k, $R_{j(k)}$ for the replicate j, nested in the environment k and $\epsilon_{ijk}$ is the residual. The corresponding heritability was defined as follow with all terms declared as random in Eq. (5):

$$h^2 = \frac{\sigma_G^2}{\sigma_G^2 + \frac{\sigma_{G \times E}^2}{t} + \frac{\sigma_r^2}{r * e}} \quad (6)$$

where $\sigma_G^2$ is the genetic variance, $\sigma_{G \times E}^2$ the $G \times E$ interaction variance, $\sigma_r^2$ the residual variance, $e$ the number of environment and $r$ the number of replicates per genotype.

## GWAS analyses

The BLUE defined at each scale (environment, envirotype, MET-22) using models (3) and (5) for SY and SN was used to perform the genetic analyses. GWAS analyses were conducted using the FastLMM algorithm (Lippert et al. 2011) using a kinship matrix calculated for each linkage group as described by Rincent et al. (2014) and following the Astle algorithm (Astle and Balding 2009). A Bonferroni detection threshold was set to 10%, based on a corrected SNP population as proposed by the simpleM method developed by Gao et al. (2008). This method was based on the composite linkage disequilibrium (CLD) correlation between SNP. The CLD was used to calculate the effective number of independent tests ($M_{eff}$). In this study, the $M_{eff}$ was calculated and fixed to 19,886 SNP. We defined the threshold of highly significant associated SNP as:

$$t = - \log_{10} \left( \frac{\alpha_G}{M_{eff}} \right) \quad (7)$$

where $t$ is the value of the value of the threshold, $\alpha_G$ is the threshold of the $p$-value and was fixed to 0.1 and $M_{eff}$ was the effective population of SNP. The final threshold used was $t=5.30$.

## Genetic diversity analysis

A genetic diversity was conducted at the associated loci detected through GWAS. This evaluated the potential impact
of a major breeding event following the selection of varieties without glucosinolates in seeds. Therefore, two populations were confronted: a first population (GSL+) composed by genotypes with high contents of glucosinolate in seeds (> 18 µmol g⁻¹) (46 WOSR, “++” and “0+”) and a second population (GSL-) composed by the remaining accessions of the P173 that present low glucosinolate content (< 18 µmol g⁻¹) (127 WOSR, “+0” and “00”) (Supplemental Data 2). The nucleotide divergence statistic (π) and the mean Fst (Weir and Cockerham 1984) for the whole P173 and for the GSL+ population and the GSL- population were calculated using VCFtools v0.1.13 software (Danecek et al. 2011) for a window of 10 kb around each genomic region of interest.

## Results

### N stress was not the only factor impacting seed yield

A first overlook of the MET-47 was carried out by analyzing the relationship between the mean values of seed yield and NNI for Aviso and Montego per environment (Fig. 1). When focusing on environments defined as N⁻, it is to be noticed than NNI values ranged from 0.63 to 1.26, with NNI values exceeding 1.0 for 7 N⁻-environments. Moreover, in two N⁻ environments (Lou18_N⁻ and Chr18_N⁻), the

observed SY was even higher than the targeted value for N⁺ conditions. Thus, these results demonstrated the difficulty of implementing a nitrogen stress under field conditions. When considering the N⁺ environments, it can be observed that the targeted seed yield of 35 q ha⁻¹ was achieved in 13 environments out of 19. However, in these 13 environments, NNI was always higher than 1 (except for LR13_N⁺), indicating the presence of other environmental factors that did impact seed yield. These results illustrated the need to consider a comprehensive environmental description of the trials to detect the factors that indeed affected SY.

### Identification of 13 environmental indicators limiting seed yield

MET-47 offered a unique chance to get insights into the SY most limiting factors. For that purpose, the set of 79 indicators was used to run a PLS regression analysis on the mean SY variation of Aviso and Montego for each of the 47 environments. SY-limiting factors corresponded to a set of 13 indicators (see Table 1 for a precise description). Predicted seed yields using the PLS regression model were strongly correlated to the observed seed yields (R²=0.87, RMSE=2.6), thus validating the set of indicators. These 13 indicators covered the whole crop growing cycle from the climatic winter (CW) to the reserve allocation to the pod period (P600). Six indicators referred to the temperature

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Relationship between mean values of seed yield (SY, Y-axis) and nitrogen nutrition index (NNI, X-axis) observed for the probe genotypes, Aviso and Montego, across the MET-47. Empty and plain circles correspond, respectively, to the N⁻ and N⁺ nitrogen fertilization regimes. A vertical line was added and corresponds to a NNI value of 1, above which environment are considered as not impacted by N stress. The horizontal bar corresponds to a SY of 35 q ha⁻¹ which is the targeted yield for N⁺ conditions</figcaption></figure>
of a major breeding event following the selection of varieties
without glucosinolates in seeds. Therefore, two populations
were confronted: a first population (GSL+) composed by
genotypes with high contents of glucosinolate in seeds (> 18
µmol g⁻¹) (46 WOSR, “++” and “0+”) and a second popula-
tion (GSL-) composed by the remaining accessions of the
P173 that present low glucosinolate content (< 18 µmol g⁻¹)
(127 WOSR, “+0” and “00”) (Supplemental Data 2). The
nucleotide divergence statistic (π) and the mean Fst (Weir
and Cockerham 1984) for the whole P173 and for the GSL+
population and the GSL- population were calculated using
VCFtools v0.1.13 software (Danecek et al. 2011) for a win-
dow of 10 kb around each genomic region of interest.

Results

N stress was not the only factor impacting seed yield

A first overlook of the MET-47 was carried out by analyz-
ing the relationship between the mean values of seed yield
and NNI for Aviso and Montego per environment (Fig. 1).
When focusing on environments defined as N⁻, it is to be
noticed than NNI values ranged from 0.63 to 1.26, with NNI
values exceeding 1.0 for 7 N⁻ -environments. Moreover,
in two N⁻ environments (Lou18_N⁻ and Chr18_N⁻), the

observed SY was even higher that the targeted value for N⁺
conditions. Thus, these results demonstrated the difficulty of
implementing a nitrogen stress under field conditions. When
considering the N⁺ environments, it can be observed that the
targeted seed yield of 35 q ha⁻¹ was achieved in 13 environ-
ments out of 19. However, in these 13 environments, NNI
was always higher than 1 (except for LR13_N⁺), indicating
the presence of other environmental factors that did impact
seed yield. These results illustrated the need to consider a
comprehensive environmental description of the trials to
detect the factors that indeed affected SY.

Identification of 13 environmental indicators
limiting seed yield

MET-47 offered a unique chance to get insights into the
SY most limiting factors. For that purpose, the set of 79
indicators was used to run a PLS regression analysis on
the mean SY variation of Aviso and Montego for each of
the 47 environments. SY-limiting factors corresponded to
a set of 13 indicators (see Table 1 for a precise description).
Predicted seed yields using the PLS regression model were
strongly correlated to the observed seed yields (R²=0.87,
RMSE=2.6), thus validating the set of indicators. These
13 indicators covered the whole crop growing cycle from
the climatic winter (CW) to the reserve allocation to the
pod period (P600). Six indicators referred to the temperature

**Fig. 1** Relationship between mean values of seed yield (SY, Y-axis) and nitrogen nutrition index (NNI, X-axis) observed for the probe genotypes, Aviso and Montego, across the MET-47. Empty and plain circles correspond, respectively, to the N⁻ and N⁺ nitrogen fertilization regimes. A vertical line was added and corresponds to a NNI value of 1, above which environment are considered as not impacted by N stress. The horizontal bar corresponds to a SY of 35 q ha⁻¹ which is the targeted yield for N⁺ conditions

<figure><img src="image_4.png" /></figure>
<table><caption>**Table 1** Seed yield limiting indicators across the MET-47</caption><thead><tr><th>Indicators</th><th>Description</th><th>Correlated indicators*</th></tr></thead><tbody><tr><td>TMN_CW</td><td>Mean temperature during CW</td><td>–</td></tr><tr><td>TMAX_B</td><td>Maximal temperature registered during B</td><td>–</td></tr><tr><td>TMAX_FLO</td><td>Maximal temperature registered during FLO</td><td>SSR_CW (0.67); LT_B (-0.69); LT_FLO (-0.66)</td></tr><tr><td>HT_FLO</td><td>Number of days with high temperature (> 25 °C) at FLO</td><td>TMIN_CW (-0.7); TMN_FLO (0.66); LGDD_FLO (0.66); TMN_P600 (0.7)</td></tr><tr><td>TMAX_P300</td><td>Maximal temperature registered during P300</td><td>LGDD_CW (0.68); SSR_CW (0.67); LSR_CW (0.7); LT_FLO (-0.68); HT_P300 (0.68)</td></tr><tr><td>TMIN_P600</td><td>Minimal temperature registered during P600</td><td>TMIN_FLO (0.73); TMN_FLO (0.7); LT_FLO (-0.67); LGDD_FLO (0.7); TMN_P600 (0.78)</td></tr><tr><td>WD_P600</td><td>Number of days with WSC=0 during P600</td><td>HT_F (0.66); WS_FLO (0.7); WSC_P300 (-0.68); WSC_P600 (-0.66)</td></tr><tr><td>WSC_MAX</td><td>Maximal water soil capacity</td><td>TMN_F (-0.67)</td></tr><tr><td>SSR_FLO</td><td>Sum of solar radiation during FLO</td><td>LSR_FLO (-0.67); TMAX_P1000 (0.69); HT_P1000 (0.65)</td></tr><tr><td>VERN</td><td>Optimal vernalization treatment</td><td>LT_CW (0.65)</td></tr><tr><td>NNI</td><td>Nitrogen nutrition index</td><td>AN (0.65); N<sub>total</sub> (0.66)</td></tr><tr><td>N<sub>total</sub></td><td>Nitrogen input</td><td>NNI (0.66)</td></tr><tr><td>DAI</td><td>Number of days of dryness after N supply</td><td>–</td></tr></tbody></table>

(TMN_CW, TMAX_B, TMAX_FLO, HT_FLO, TMAX_P300, and TMIN_P600), two to the water status (WSC_MAX and WD_P600), one to the radiation (SSR_FLO), one to vernalization conditions (VERN) and three to nitrogen status (NNI, N<sub>total</sub> and DAI; Table 2). Finally, the study of the Pearson correlations (Table 2, Supplemental Fig. 1) showed important correlations between the 13 identified indicators and indicators related to the reproductive periods of winter oilseed rape (FLO to P1000) and to the climatic winter period (Table 1), highlighting the key role of these periods on the elaboration of the final seed yield.

## MET-22 included four contrasting envirotypes

The MET-22 is a subset of the MET-47 where 173 genotypes were trialed. The mean seed yield of the MET-22 was 30.35 q ha⁻¹ (29.40 q ha⁻¹ for the MET-47), and the mean NNI

across the MET-22 was 1.06 (1.02 for the MET-47). The MET-22 description using the limiting indicators previously identified presented a similar profile to the MET-47 (Supplemental Fig. 2). A hierarchical clustering was performed on the environments of the MET-22 to define envirotypes according to their pattern of limiting indicators. Four envirotypes were defined (EA, EB, EC and ED; Fig. 2, Supplemental Data 1). Envirotype EA is composed of three N⁺ environments and one N⁻ environment. The mean seed yield observed across EA was 37.6 q ha⁻¹, and the mean NNI was 1.29. This envirotype can therefore be defined as non N stressed and high yielding. The analysis of the limiting indicators pattern showed that EA, when compared to the mean MET-22, is characterized by lower temperature during winter, lower heat stress and solar radiation at flowering (HT_FLO), a higher temperature during the seed number fixation period (TMAX_P300) and a higher NNI (Fig. 3).

<table><caption>**Table 2** Results of the mixed linear model applied on the MET-22 to evaluate the impact of the envirotyping on variance repartition</caption><thead><tr><th rowspan="2">Trait</th><th>Models</th><th colspan="5">% Variance</th></tr><tr><th>Models (5)</th><th>G</th><th>E</th><th>G×E</th><th>R</th><th>ε</th></tr></thead><tbody><tr><td>SY</td><td></td><td>24.7</td><td>55.5</td><td>5.8</td><td>3.8</td><td>10.1</td></tr><tr><td></td><td>Model (2)</td><td>G</td><td>En</td><td>En×E</td><td>G×En</td><td>G×En×E</td><td>R</td><td>ε</td></tr><tr><td>SY</td><td></td><td>20</td><td>44.4</td><td>19.2</td><td>1.2</td><td>3.9</td><td>3.1</td><td>8.2</td></tr></tbody></table>

Data had been analyzed using linear model (2): $Y_{ijkl} = \mu + G_i + En_i + En_{l(k)} + G_i \times En_i + G_i \times En_{l(k)} + R_{l(k)} + \epsilon_{ijkl}$; with the genotype (G), the envirotype (En), the environment nested in the envirotype (En<sub>l(k)</sub> or En×E), the genotype by envirotype interaction (G×En), the genotype by environment by envirotype interaction (G<sub>i</sub>×En<sub>l(k)</sub>) or G×En×E) and the block effect nested in the environment, nested in the envirotype (R) and compared with the results of linear model (5): $Y_{ijk} = \mu + G_i + E_k + G_i \times E_k + R_{j(k)} + \epsilon_{ijk}$, with the genotype (G), the environment (E), the genotype by environment interaction (G×E) and the block (R) nested in the environment

<footnote>*Value of the Pearson correlation coefficient was indicated between brackets</footnote>
<table><caption>**Table 1** Seed yield limiting indicators across the MET-47</caption><thead><tr><th>Indicators</th><th>Description</th><th>Correlated indicators*</th></tr></thead><tbody><tr><td>TMN_CW</td><td>Mean temperature during CW</td><td>–</td></tr><tr><td>TMAX_B</td><td>Maximal temperature registered during B</td><td>–</td></tr><tr><td>TMAX_FLO</td><td>Maximal temperature registered during FLO</td><td>SSR_CW (0.67); LT_B (-0.69); LT_FLO (-0.66)</td></tr><tr><td>HT_FLO</td><td>Number of days with high temperature (> 25 °C) at FLO</td><td>TMIN_CW (-0.7); TMN_FLO (0.66); LGDD_FLO (0.66); TMN_P600 (0.7)</td></tr><tr><td>TMAX_P300</td><td>Maximal temperature registered during P300</td><td>LGDD_CW (0.68); SSR_CW (0.67); LSR_CW (0.7); LT_FLO (-0.68); HT_P300 (0.68)</td></tr><tr><td>TMIN_P600</td><td>Minimal temperature registered during P600</td><td>TMIN_FLO (0.73); TMN_FLO (0.7); LT_FLO (-0.67); LGDD_FLO (0.7); TMN_P600 (0.78)</td></tr><tr><td>WD_P600</td><td>Number of days with WSC=0 during P600</td><td>HT_F (0.66); WS_FLO (0.7); WSC_P300 (-0.68); WSC_P600 (-0.66)</td></tr><tr><td>WSC_MAX</td><td>Maximal water soil capacity</td><td>TMN_F (-0.67)</td></tr><tr><td>SSR_FLO</td><td>Sum of solar radiation during FLO</td><td>LSR_FLO (-0.67); TMAX_P1000 (0.69); HT_P1000 (0.65)</td></tr><tr><td>VERN</td><td>Optimal vernalization treatment</td><td>LT_CW (0.65)</td></tr><tr><td>NNI</td><td>Nitrogen nutrition index</td><td>AN (0.65); N<sub>total</sub> (0.66)</td></tr><tr><td>N<sub>total</sub></td><td>Nitrogen input</td><td>NNI (0.66)</td></tr><tr><td>DAI</td><td>Number of days of dryness after N supply</td><td>–</td></tr></tbody></table>

(TMN_CW, TMAX_B, TMAX_FLO, HT_FLO, TMAX_P300, and TMIN_P600), two to the water status (WSC_MAX and WD_P600), one to the radiation (SSR_FLO), one to vernalization conditions (VERN) and three to nitrogen status (NNI, N<sub>total</sub> and DAI; Table 2). Finally, the study of the Pearson correlations (Table 2, Supplemental Fig. 1) showed important correlations between the 13 identified indicators and indicators related to the reproductive periods of winter oilseed rape (FLO to P1000) and to the climatic winter period (Table 1), highlighting the key role of these periods on the elaboration of the final seed yield.

## MET-22 included four contrasting envirotypes

The MET-22 is a subset of the MET-47 where 173 genotypes were trialed. The mean seed yield of the MET-22 was 30.35 q ha⁻¹ (29.40 q ha⁻¹ for the MET-47), and the mean NNI

across the MET-22 was 1.06 (1.02 for the MET-47). The MET-22 description using the limiting indicators previously identified presented a similar profile to the MET-47 (Supplemental Fig. 2). A hierarchical clustering was performed on the environments of the MET-22 to define envirotypes according to their pattern of limiting indicators. Four envirotypes were defined (EA, EB, EC and ED; Fig. 2, Supplemental Data 1). Envirotype EA is composed of three N⁺ environments and one N⁻ environment. The mean seed yield observed across EA was 37.6 q ha⁻¹, and the mean NNI was 1.29. This envirotype can therefore be defined as non N stressed and high yielding. The analysis of the limiting indicators pattern showed that EA, when compared to the mean MET-22, is characterized by lower temperature during winter, lower heat stress and solar radiation at flowering (HT_FLO), a higher temperature during the seed number fixation period (TMAX_P300) and a higher NNI (Fig. 3).

<table><caption>**Table 2** Results of the mixed linear model applied on the MET-22 to evaluate the impact of the envirotyping on variance repartition</caption><thead><tr><th rowspan="2">Trait</th><th>Models</th><th colspan="5">% Variance</th></tr><tr><th>Models (5)</th><th>G</th><th>E</th><th>G×E</th><th>R</th><th>ε</th></tr></thead><tbody><tr><td>SY</td><td></td><td>24.7</td><td>55.5</td><td>5.8</td><td>3.8</td><td>10.1</td></tr><tr><td></td><td>Model (2)</td><td>G</td><td>En</td><td>En×E</td><td>G×En</td><td>G×En×E</td><td>R</td><td>ε</td></tr><tr><td>SY</td><td></td><td>20</td><td>44.4</td><td>19.2</td><td>1.2</td><td>3.9</td><td>3.1</td><td>8.2</td></tr></tbody></table>

Data had been analyzed using linear model (2): $Y_{ijkl} = \mu + G_i + En_l + En_{l(k)} + G_i \times En_l + G_i \times En_{l(k)} + R_{l(k)} + \varepsilon_{ijkl}$; with the genotype (G), the envirotype (En), the environment nested in the envirotype (En<sub>l(k)</sub> or En×E), the genotype by envirotype interaction (G×En), the genotype by environment by envirotype interaction (G<sub>i</sub>×En<sub>l(k)</sub>) or G×En×E) and the block effect nested in the environment, nested in the envirotype (R) and compared with the results of linear model (5): $Y_{ijk} = \mu + G_i + E_k + G_i \times E_k + R_{j(k)} + \varepsilon_{ijk}$, with the genotype (G), the environment (E), the genotype by environment interaction (G×E) and the block (R) nested in the environment

<footnote>*Value of the Pearson correlation coefficient was indicated between brackets</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Envirotyping of the MET-22. **a** Dendrogram tree of the 22 environments of the MET-22 based on the PLS regression results. Mean seed yield (SY) and mean nitrogen nutrition index (NNI) values are indicated for each envirotype. **b** Variation of the inertia gain depending on the number of groups chosen for the hierarchical clustering. The dashed line represents the choice of groups number</figcaption></figure>

In contrast, envirotype ED was characterized by the lowest mean seed yield (18.5 q ha⁻¹) and a mean NNI of 0.91 suggesting that the environments of this envirotype have been more stressed. ED gathered three N⁻ environments and one N⁺ environment. The limiting indicators pattern of envirotype ED is characterized by poor vernalization conditions (VERN), higher temperature during winter (TMN_CW), lower SSR_FLO and TMAX_P300 than the mean MET-22 (Fig. 3).

Finally, the two remaining envirotypes, EB and EC, presented mean seed yields of 28.4 and 27.2 q ha⁻¹, respectively, while the mean MET-22 seed yield was 30.34 q ha⁻¹. Envirotype EB was composed of eight environments (three N⁻ and five N⁺) and envirotype EC of six environments (four N⁻ and two N⁺). Envirotypes EC and EB distinguished from each other according to NNI (1.22 for EB and 0.92 for EC), Ntotal, VRN, and conditions at flowering (HT_FLO and SSR_FLO). The profiles of limiting indicators showed that TMAX_FLO, HT_FLO, SSR_FLO, TMAX_P300, WSC_MAX and TMN_CW were lower in envirotype EB compared to the MET-22. Envirotype EC was characterized by higher HT_FLO, SSR_FLO, WSC_MAX and TMN_CW while DAI, NNI, Ntotal, TMAX_B and VERN were lower than observed for the MET-22 (Fig. 3).

Our results demonstrated that the 22 environments of the MET-22 could be classified in four envirotypes: envirotype

EA that can be considered as not stressed, envirotype ED considered as the most stressed, while envirotypes EB and EC showed moderate stress and mean SY similar to the mean SY observed at the MET-22 scale, but with a different pattern of limiting factors.

## Envirotyping explained up to 70% of environmental variation and 24.6% of G×E

To evaluate the impact of envirotyping on the variance distribution of the E and G×E effects, a linear model (2) was applied to SY. At the MET-22 scale, the envirotype effect (En) explained a main part of the environmental effect (70% for SY) as shown by the comparison of models (2) and (5) (Table 2). Genotype by envirotype interaction (G×En) explained 24.6% of the G×E interaction observed at the MET-22 scale. Within each envirotype, the G×E variation was reduced in envirotype EA (2.3% of total variation) when compared to the MET-22 (5.8%), but increased for envirotypes EB (8.7%) and particularly for ED (15.9%) (Table 3). As SY and SN were highly correlated ($r=0.94$, p-value < 0.001), the envirotyping based on SY was used to study SN variation and its repartition. Similar results were observed and are reported on Table 4.

Heritabilities were calculated for SY and SN for each environment, for each envirotype and for the mean
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Envirotyping of the MET-22. **a** Dendrogram tree of the 22 environments of the MET-22 based on the PLS regression results. Mean seed yield (SY) and mean nitrogen nutrition index (NNI) values are indicated for each envirotype. **b** Variation of the inertia gain depending on the number of groups chosen for the hierarchical clustering. The dashed line represents the choice of groups number</figcaption></figure>

In contrast, envirotype ED was characterized by the lowest mean seed yield (18.5 q ha⁻¹) and a mean NNI of 0.91 suggesting that the environments of this envirotype have been more stressed. ED gathered three N⁻ environments and one N⁺ environment. The limiting indicators pattern of envirotype ED is characterized by poor vernalization conditions (VERN), higher temperature during winter (TMN_CW), lower SSR_FLO and TMAX_P300 than the mean MET-22 (Fig. 3).

Finally, the two remaining envirotypes, EB and EC, presented mean seed yields of 28.4 and 27.2 q ha⁻¹, respectively, while the mean MET-22 seed yield was 30.34 q ha⁻¹. Envirotype EB was composed of eight environments (three N⁻ and five N⁺) and envirotype EC of six environments (four N⁻ and two N⁺). Envirotypes EC and EB distinguished from each other according to NNI (1.22 for EB and 0.92 for EC), Ntotal, VRN, and conditions at flowering (HT_FLO and SSR_FLO). The profiles of limiting indicators showed that TMAX_FLO, HT_FLO, SSR_FLO, TMAX_P300, WSC_MAX and TMN_CW were lower in envirotype EB compared to the MET-22. Envirotype EC was characterized by higher HT_FLO, SSR_FLO, WSC_MAX and TMN_CW while DAI, NNI, Ntotal, TMAX_B and VERN were lower than observed for the MET-22 (Fig. 3).

Our results demonstrated that the 22 environments of the MET-22 could be classified in four envirotypes: envirotype

EA that can be considered as not stressed, envirotype ED considered as the most stressed, while envirotypes EB and EC showed moderate stress and mean SY similar to the mean SY observed at the MET-22 scale, but with a different pattern of limiting factors.

## Envirotyping explained up to 70% of environmental variation and 24.6% of G×E

To evaluate the impact of envirotyping on the variance distribution of the E and G×E effects, a linear model (2) was applied to SY. At the MET-22 scale, the envirotype effect (En) explained a main part of the environmental effect (70% for SY) as shown by the comparison of models (2) and (5) (Table 2). Genotype by envirotype interaction (G×En) explained 24.6% of the G×E interaction observed at the MET-22 scale. Within each envirotype, the G×E variation was reduced in envirotype EA (2.3% of total variation) when compared to the MET-22 (5.8%), but increased for envirotypes EB (8.7%) and particularly for ED (15.9%) (Table 3). As SY and SN were highly correlated ($r=0.94$, p-value < 0.001), the envirotyping based on SY was used to study SN variation and its repartition. Similar results were observed and are reported on Table 4.

Heritabilities were calculated for SY and SN for each environment, for each envirotype and for the mean
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Envirotype description based on the seed yield limiting factors. The profile of each limiting factor is represented on a 0–10 scale for each envirotype as well as for the MET-22. Abbreviations of the indicators are as following: *TMN_CW* mean temperature at climatic winter, *VERN* vernalization condition, *TMAX_B* maximal temperature at bolting, *HT_FLO* number of days with high temperature (> 25 °C)</figcaption></figure>

at flowering, *SSR_FLO* sum of solar radiation at flowering, *TMAX_P300* maximal temperature during the seed number fixation period, *TMIN_P600* minimal temperature during the pod growth, *WD_P600* water deficit during the pod growth, *NNI* nitrogen nutrition index, *N_total* amount of nitrogen bring to the environment, *DAI* number of dry days after the N⁻ input, *WSC_MAX* maximal water soil capacity

<table><caption>**Table 3** Results of the mixed linear model applied for the MET-22 and for each envirotype considering all trials as confounded [linear model (5)] and heritabilities estimation</caption><thead><tr><th rowspan="2">Trait</th><th rowspan="2">Group</th><th rowspan="2">Nb of environment<sup>b</sup></th><th colspan="5">% Variance<sup>a</sup></th><th rowspan="2">h<sup>2</sup></th></tr><tr><th>G</th><th>E</th><th>G×E</th><th>R</th><th>ε</th></tr></thead><tbody><tr><td rowspan="5">SY</td><td>MET-22</td><td>22</td><td>24.7</td><td>55.5</td><td>5.8</td><td>3.8</td><td>10.1</td><td>0.98</td></tr><tr><td>Envirotype EA</td><td>4</td><td>37.7</td><td>48.3</td><td>2.3</td><td>3.8</td><td>7.8</td><td>0.96</td></tr><tr><td>Envirotype EB</td><td>8</td><td>38.5</td><td>34</td><td>8.7</td><td>4.8</td><td>14.1</td><td>0.95</td></tr><tr><td>Envirotype EC</td><td>6</td><td>45.7</td><td>8.7</td><td>5.9</td><td>10.8</td><td>29</td><td>0.94</td></tr><tr><td>Envirotype ED</td><td>4</td><td>28.9</td><td>34.8</td><td>15.9</td><td>5.3</td><td>15</td><td>0.83</td></tr><tr><td rowspan="5">SN</td><td>MET-22</td><td>18</td><td>33</td><td>40.4</td><td>8.1</td><td>4.9</td><td>13.5</td><td>0.98</td></tr><tr><td>Envirotype EA</td><td>2</td><td>46.9</td><td>34.1</td><td>8.1</td><td>1.1</td><td>9.9</td><td>0.88</td></tr><tr><td>Envirotype EB</td><td>7</td><td>49</td><td>27</td><td>9.4</td><td>2.1</td><td>12.5</td><td>0.96</td></tr><tr><td>Envirotype EC</td><td>6</td><td>43.8</td><td>6.9</td><td>9.1</td><td>13.7</td><td>26.5</td><td>0.93</td></tr><tr><td>Envirotype ED</td><td>3</td><td>29.6</td><td>32.5</td><td>11.5</td><td>7.4</td><td>19</td><td>0.81</td></tr></tbody></table>

Linear model (5): $Y_{ijk} = \mu + G_i + E_k + G_i \times E_k + R_{j(k)} + \varepsilon_{ijk}$, with the genotype (G), the environment (E), the genotype by environment interaction ($G \times E$) and the block (R) nested in the trial. Heritabilities ($h^2$) were estimated using Eq. (6): $h^2 = \frac{\sigma_G^2}{\sigma_G^2 + \frac{\sigma_{G \times E}^2}{r+\varepsilon}}$

MET-22. Into more details, SY and SN were highly heritable (0.98 for the MET-22 for both traits) and heritabilities ranged from 0.83 to 0.96 for SY and from 0.81 to 0.86 for SN depending on the envirotype (Table 3). At the

environment scale, heritabilities ranged from 0.63 to 0.96 for SY and from 0.53 to 0.86 for SN (Table 4). For most of the environments, SY and SN were highly heritable.

<footnote><sup>a</sup>Percentage of variance explained by each parameter</footnote>

<footnote><sup>b</sup>Number of trials considered for the evaluation of the different effects and estimation of heritabilities</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Envirotype description based on the seed yield limiting factors. The profile of each limiting factor is represented on a 0–10 scale for each envirotype as well as for the MET-22. Abbreviations of the indicators are as following: *TMN_CW* mean temperature at climatic winter, *VERN* vernalization condition, *TMAX_B* maximal temperature at bolting, *HT_FLO* number of days with high temperature (> 25 °C)</figcaption></figure>

at flowering, *SSR_FLO* sum of solar radiation at flowering, *TMAX_P300* maximal temperature during the seed number fixation period, *TMIN_P600* minimal temperature during the pod growth, *WD_P600* water deficit during the pod growth, *NNI* nitrogen nutrition index, *N_total* amount of nitrogen bring to the environment, *DAI* number of dry days after the N⁻ input, *WSC_MAX* maximal water soil capacity

<table><caption>**Table 3** Results of the mixed linear model applied for the MET-22 and for each envirotype considering all trials as confounded [linear model (5)] and heritabilities estimation</caption><thead><tr><th rowspan="2">Trait</th><th rowspan="2">Group</th><th rowspan="2">Nb of environment<sup>b</sup></th><th colspan="5">% Variance<sup>a</sup></th><th rowspan="2">h<sup>2</sup></th></tr><tr><th>G</th><th>E</th><th>G×E</th><th>R</th><th>ε</th></tr></thead><tbody><tr><td rowspan="5">SY</td><td>MET-22</td><td>22</td><td>24.7</td><td>55.5</td><td>5.8</td><td>3.8</td><td>10.1</td><td>0.98</td></tr><tr><td>Envirotype EA</td><td>4</td><td>37.7</td><td>48.3</td><td>2.3</td><td>3.8</td><td>7.8</td><td>0.96</td></tr><tr><td>Envirotype EB</td><td>8</td><td>38.5</td><td>34</td><td>8.7</td><td>4.8</td><td>14.1</td><td>0.95</td></tr><tr><td>Envirotype EC</td><td>6</td><td>45.7</td><td>8.7</td><td>5.9</td><td>10.8</td><td>29</td><td>0.94</td></tr><tr><td>Envirotype ED</td><td>4</td><td>28.9</td><td>34.8</td><td>15.9</td><td>5.3</td><td>15</td><td>0.83</td></tr><tr><td rowspan="5">SN</td><td>MET-22</td><td>18</td><td>33</td><td>40.4</td><td>8.1</td><td>4.9</td><td>13.5</td><td>0.98</td></tr><tr><td>Envirotype EA</td><td>2</td><td>46.9</td><td>34.1</td><td>8.1</td><td>1.1</td><td>9.9</td><td>0.88</td></tr><tr><td>Envirotype EB</td><td>7</td><td>49</td><td>27</td><td>9.4</td><td>2.1</td><td>12.5</td><td>0.96</td></tr><tr><td>Envirotype EC</td><td>6</td><td>43.8</td><td>6.9</td><td>9.1</td><td>13.7</td><td>26.5</td><td>0.93</td></tr><tr><td>Envirotype ED</td><td>3</td><td>29.6</td><td>32.5</td><td>11.5</td><td>7.4</td><td>19</td><td>0.81</td></tr></tbody></table>

Linear model (5): $Y_{ijk} = \mu + G_i + E_k + G_i \times E_k + R_{j(k)} + \epsilon_{ijk}$, with the genotype (G), the environment (E), the genotype by environment interaction ($G \times E$) and the block (R) nested in the trial. Heritabilities ($h^2$) were estimated using Eq. (6): $h^2 = \frac{\sigma_G^2}{\sigma_G^2 + \frac{\sigma_{G \times E}^2}{r+\epsilon}}$

MET-22. Into more details, SY and SN were highly heritable (0.98 for the MET-22 for both traits) and heritabilities ranged from 0.83 to 0.96 for SY and from 0.81 to 0.86 for SN depending on the envirotype (Table 3). At the

environment scale, heritabilities ranged from 0.63 to 0.96 for SY and from 0.53 to 0.86 for SN (Table 4). For most of the environments, SY and SN were highly heritable.

<footnote><sup>a</sup>Percentage of variance explained by each parameter</footnote>

<footnote><sup>b</sup>Number of trials considered for the evaluation of the different effects and estimation of heritabilities</footnote>
<table><caption>**Table 4** Heritabilities estimated for each trait and each environment of the MET-22</caption><thead><tr><th rowspan="2">Envirotype</th><th rowspan="2">Environment</th><th colspan="2">h<sup>2</sup></th></tr><tr><th>SY</th><th>SN</th></tr></thead><tbody><tr><td rowspan="4">EA</td><td>Dij15_N<sup>-</sup></td><td>0.90</td><td>-</td></tr><tr><td>Dij15_N<sup>+</sup></td><td>0.91</td><td>-</td></tr><tr><td>Ver15_N<sup>+</sup></td><td>0.91</td><td>-</td></tr><tr><td>Yeb15_N<sup>+</sup></td><td>0.92</td><td>0.93</td></tr><tr><td rowspan="8">EB</td><td>Ch14_N<sup>+</sup></td><td>0.93</td><td>0.93</td></tr><tr><td>Pre15_N<sup>-</sup></td><td>0.81</td><td>0.85</td></tr><tr><td>Pre15_N<sup>+</sup></td><td>0.80</td><td>0.85</td></tr><tr><td>Sel14_N<sup>-</sup></td><td>0.94</td><td>0.93</td></tr><tr><td>Sel14_N<sup>+</sup></td><td>0.95</td><td>0.93</td></tr><tr><td>Sel15_N<sup>+</sup></td><td>0.63</td><td>-</td></tr><tr><td>Ver14_N<sup>+</sup></td><td>0.91</td><td>0.9</td></tr><tr><td>Ver15_N<sup>-</sup></td><td>0.86</td><td>-</td></tr><tr><td rowspan="8">EC</td><td>LR15_N<sup>-</sup></td><td>0.85</td><td>0.88</td></tr><tr><td>LR15_N<sup>+</sup></td><td>0.87</td><td>0.9</td></tr><tr><td>Md15_N<sup>-</sup></td><td>0.67</td><td>0.53</td></tr><tr><td>Md15_N<sup>+</sup></td><td>0.74</td><td>0.65</td></tr><tr><td>Ver14_N<sup>-</sup></td><td>0.77</td><td>0.71</td></tr><tr><td>Yeb15_N<sup>-</sup></td><td>0.95</td><td>0.96</td></tr><tr><td>Ch14_N<sup>-</sup></td><td>0.81</td><td>0.78</td></tr><tr><td>Pre14_N<sup>-</sup></td><td>0.68</td><td>0.75</td></tr><tr><td rowspan="3">ED</td><td>Pre14_N<sup>+</sup></td><td>0.87</td><td>0.87</td></tr><tr><td>Sel15_N<sup>-</sup></td><td>0.96</td><td>-</td></tr></tbody></table>

Components of the heritabilities were obtained using the linear model (3): $Y_{ij} = \mu + G_i + R_j + \epsilon_{ij}$; Heritabilities ($h^2$) were calculated according to Eq. (4): $h^2 = \frac{\sigma_G^2}{\sigma_G^2 + \sigma_T^2}$

## Most of the SY-related QTL detected at the environment scale were specific of a single environment

Comparison of SY-QTL detected for each environment is a first step to characterize the genetic determinism of SY stability. Therefore, the BLUE obtained for each genotype and each environment using model (3) was used as input for GWAS analysis. A total of 87 SNP were detected showing a significant association with SY or SN considering all the analyses performed for each environment (Supplemental Data 4). These SNP were grouped into 11 QTL related to SY (QA03, QA05a, QA07a, QA07b, QA09, QC02, QC03a, QC04, QC08, QC09a and QC09b) and four QTL related to SN (QA05b, QA07a, QA09 and QA10) (Table 5). For almost all loci, the most frequent allele in the population was the favorable one except for QA03, QA07a and QC09b. Considering both traits, it is worth noting that all QTLs were specific to a single environment

except QA07a, QA09, QA10, QC02, QC04 and QC09a. QA07a, QA10, QC02 and QC04 were only detected in two environments. The QA09 and QCA09a were detected in a large range of environments with six environments for QA09 and seven for QC09a, respectively.

## Envirotyping highlighted five QTLs

BLUE calculated for each genotype and for each of the four envirotypes was used for GWAS analysis. In this way, the detected QTL will be linked to the envirotype characterization ran upon through their profile of limiting indicators.

QTL QC09a (Table 5) was detected for SY and SN and for three envirotypes (EA, EB and EC). It was also detected for SY in seven environments as previously mentioned (Table 5, Supplemental Data 4). It explained between 9.5 and 16.7% of SY variation depending on environment or envirotype; and 20.4% of SN variation for envirotype EA and 18.8% of SN variation for envirotype EB. The major allele was favorable to increase both traits (Supplemental Data 4).

The SN QTL QC03b (Table 5) was specifically detected in envirotype EC. It was not detected at the environmental scale. This QTL explained 21% of the SN variation ($R^2=0.21$) in envirotype EC, and the major allele was the favorable one. However, confidence in this QTL is low as it is composed of a single SNP (Supplemental Data 4).

QA07a was detected for SN in envirotype ED on chromosome A07 (Table 5). This QTL was also detected in environment LR15_N<sup>+</sup> (EC) and environment Pre14_N<sup>-</sup> (ED) for SN and for SY in LR15_N<sup>+</sup> (EC) (Table 5). This QTL explained 19.5% of the SN variation in envirotype ED, and the minor allele of the population was the favorable one. (Supplemental Data 4). QA07a could be considered as specific to stressed conditions.

Lastly, the QA09 and QC09a QTL, that were detected for a wide range of environments, were also detected for several envirotypes: EA and ED for QA09 and EA, EB and EC for QC09a.

## QA09 was the unique QTL detected across the MET-22

BLUE obtained for each genotype at the whole MET-22 scale was used to detect consistent QTL controlling SY and SN across environments. A set of 15 SNP was detected for SY on chromosome A09 (Supplemental Data 4) and consisted in the QTL QA09 (Table 5). QA09 explained 12% of the SY variation, and the favorable allele was the major allele. As previously shown, QA09 was also detected for SY at the envirotype scale for EA and ED (Table 5), and at the environment scale in six environments that belonged to EA
<table><caption>**Table 4** Heritabilities estimated for each trait and each environment of the MET-22</caption><thead><tr><th rowspan="2">Envirotype</th><th rowspan="2">Environment</th><th colspan="2">h<sup>2</sup></th></tr><tr><th>SY</th><th>SN</th></tr></thead><tbody><tr><td rowspan="4">EA</td><td>Dij15_N<sup>-</sup></td><td>0.90</td><td>-</td></tr><tr><td>Dij15_N<sup>+</sup></td><td>0.91</td><td>-</td></tr><tr><td>Ver15_N<sup>+</sup></td><td>0.91</td><td>-</td></tr><tr><td>Yeb15_N<sup>+</sup></td><td>0.92</td><td>0.93</td></tr><tr><td rowspan="8">EB</td><td>Ch14_N<sup>+</sup></td><td>0.93</td><td>0.93</td></tr><tr><td>Pre15_N<sup>-</sup></td><td>0.81</td><td>0.85</td></tr><tr><td>Pre15_N<sup>+</sup></td><td>0.80</td><td>0.85</td></tr><tr><td>Sel14_N<sup>-</sup></td><td>0.94</td><td>0.93</td></tr><tr><td>Sel14_N<sup>+</sup></td><td>0.95</td><td>0.93</td></tr><tr><td>Sel15_N<sup>+</sup></td><td>0.63</td><td>-</td></tr><tr><td>Ver14_N<sup>+</sup></td><td>0.91</td><td>0.9</td></tr><tr><td>Ver15_N<sup>-</sup></td><td>0.86</td><td>-</td></tr><tr><td rowspan="8">EC</td><td>LR15_N<sup>-</sup></td><td>0.85</td><td>0.88</td></tr><tr><td>LR15_N<sup>+</sup></td><td>0.87</td><td>0.9</td></tr><tr><td>Md15_N<sup>-</sup></td><td>0.67</td><td>0.53</td></tr><tr><td>Md15_N<sup>+</sup></td><td>0.74</td><td>0.65</td></tr><tr><td>Ver14_N<sup>-</sup></td><td>0.77</td><td>0.71</td></tr><tr><td>Yeb15_N<sup>-</sup></td><td>0.95</td><td>0.96</td></tr><tr><td>Ch14_N<sup>-</sup></td><td>0.81</td><td>0.78</td></tr><tr><td>Pre14_N<sup>-</sup></td><td>0.68</td><td>0.75</td></tr><tr><td rowspan="3">ED</td><td>Pre14_N<sup>+</sup></td><td>0.87</td><td>0.87</td></tr><tr><td>Sel15_N<sup>-</sup></td><td>0.96</td><td>-</td></tr></tbody></table>

Components of the heritabilities were obtained using the linear model (3): $Y_{ij} = \mu + G_i + R_j + \epsilon_{ij}$; Heritabilities ($h^2$) were calculated according to Eq. (4): $h^2 = \frac{\sigma_G^2}{\sigma_G^2 + \sigma_T^2}$

## Most of the SY-related QTL detected at the environment scale were specific of a single environment

Comparison of SY-QTL detected for each environment is a first step to characterize the genetic determinism of SY stability. Therefore, the BLUE obtained for each genotype and each environment using model (3) was used as input for GWAS analysis. A total of 87 SNP were detected showing a significant association with SY or SN considering all the analyses performed for each environment (Supplemental Data 4). These SNP were grouped into 11 QTL related to SY (QA03, QA05a, QA07a, QA07b, QA09, QC02, QC03a, QC04, QC08, QC09a and QC09b) and four QTL related to SN (QA05b, QA07a, QA09 and QA10) (Table 5). For almost all loci, the most frequent allele in the population was the favorable one except for QA03, QA07a and QC09b. Considering both traits, it is worth noting that all QTLs were specific to a single environment

except QA07a, QA09, QA10, QC02, QC04 and QC09a. QA07a, QA10, QC02 and QC04 were only detected in two environments. The QA09 and QCA09a were detected in a large range of environments with six environments for QA09 and seven for QC09a, respectively.

## Envirotyping highlighted five QTLs

BLUE calculated for each genotype and for each of the four envirotypes was used for GWAS analysis. In this way, the detected QTL will be linked to the envirotype characterization ran upon through their profile of limiting indicators.

QTL QC09a (Table 5) was detected for SY and SN and for three envirotypes (EA, EB and EC). It was also detected for SY in seven environments as previously mentioned (Table 5, Supplemental Data 4). It explained between 9.5 and 16.7% of SY variation depending on environment or envirotype; and 20.4% of SN variation for envirotype EA and 18.8% of SN variation for envirotype EB. The major allele was favorable to increase both traits (Supplemental Data 4).

The SN QTL QC03b (Table 5) was specifically detected in envirotype EC. It was not detected at the environmental scale. This QTL explained 21% of the SN variation ($R^2=0.21$) in envirotype EC, and the major allele was the favorable one. However, confidence in this QTL is low as it is composed of a single SNP (Supplemental Data 4).

QA07a was detected for SN in envirotype ED on chromosome A07 (Table 5). This QTL was also detected in environment LR15_N<sup>+</sup> (EC) and environment Pre14_N<sup>-</sup> (ED) for SN and for SY in LR15_N<sup>+</sup> (EC) (Table 5). This QTL explained 19.5% of the SN variation in envirotype ED, and the minor allele of the population was the favorable one. (Supplemental Data 4). QA07a could be considered as specific to stressed conditions.

Lastly, the QA09 and QC09a QTL, that were detected for a wide range of environments, were also detected for several envirotypes: EA and ED for QA09 and EA, EB and EC for QC09a.

## QA09 was the unique QTL detected across the MET-22

BLUE obtained for each genotype at the whole MET-22 scale was used to detect consistent QTL controlling SY and SN across environments. A set of 15 SNP was detected for SY on chromosome A09 (Supplemental Data 4) and consisted in the QTL QA09 (Table 5). QA09 explained 12% of the SY variation, and the favorable allele was the major allele. As previously shown, QA09 was also detected for SY at the envirotype scale for EA and ED (Table 5), and at the environment scale in six environments that belonged to EA
<table><caption>**Table 5** Description of the QTL detected for seed yield and seed number stability/instability</caption><thead><tr><th>Name</th><th>Chromosome</th><th>Positions</th><th>Traits</th><th>Environments</th><th>Envirotypes</th><th>MET-22</th></tr></thead><tbody><tr><td>QA03</td><td>A03</td><td>6,452,989</td><td>SY</td><td>Pre15_N<sup>+(EB)</sup></td><td>—</td><td>—</td></tr><tr><td>QA05a</td><td>A05</td><td>2,078,385–2,080,519</td><td>SY</td><td>Chr14_N<sup>+(EB)</sup></td><td>—</td><td>—</td></tr><tr><td>QA05b</td><td>A05</td><td>3,630,228</td><td>SN</td><td>Pre14_N<sup>+(ED)</sup></td><td>—</td><td>—</td></tr><tr><td rowspan="2">QA07a</td><td rowspan="2">A07</td><td>138,830–1,009,271</td><td>SN</td><td>LR15_N<sup>+(EC)</sup>/Pre14_N<sup>-(ED)</sup></td><td rowspan="2">ED</td><td rowspan="2">—</td></tr><tr><td>138,830–1,009,271</td><td>SY</td><td>LR15_N<sup>+(EC)</sup></td></tr><tr><td>QA07b</td><td>A07</td><td>20,033,023–20,033,025</td><td>SY</td><td>LR15_N<sup>+(EC)</sup></td><td>—</td><td>—</td></tr><tr><td>QA09</td><td>A09</td><td>3,933,991–4,504,817</td><td>SY/SN</td><td>Yeb15_N<sup>+(EA)</sup>/Ver14_N<sup>+(EB)</sup> Ch14_N<sup>-(ED)</sup>/Ch14_N<sup>+(EB)</sup> Dij15_N<sup>-(EA)</sup>/LR15_N<sup>+(EC)</sup></td><td>EA/ED</td><td>Yes</td></tr><tr><td>QA10</td><td>A10</td><td>2,515,307</td><td>SN</td><td>Pre14_N<sup>-(ED)</sup>/Ver14_N<sup>-(EC)</sup></td><td>—</td><td>—</td></tr><tr><td>QC02</td><td>C02</td><td>60,063,320–60,615,209</td><td>SY</td><td>Ver14_N<sup>+(EB)</sup>/Yeb15_N<sup>+(EA)</sup></td><td>—</td><td>—</td></tr><tr><td>QC03a</td><td>C03</td><td>20,140,168–20,809,663</td><td>SY</td><td>Sel15_N<sup>+(EB)</sup></td><td>—</td><td>—</td></tr><tr><td>QC03b</td><td>C03</td><td>23,645,808</td><td>SN</td><td>—</td><td>EC</td><td>—</td></tr><tr><td>QC04</td><td>C04</td><td>6,033,108–6,289,953</td><td>SY</td><td>Dij15_N<sup>+(EA)</sup>/Pre15_N<sup>-(EB)</sup></td><td>—</td><td>—</td></tr><tr><td>QC08</td><td>C08</td><td>37,326,511–37,326,814</td><td>SY</td><td>Ver14_N<sup>+(EB)</sup></td><td>—</td><td>—</td></tr><tr><td>QC09a</td><td>C09</td><td>2,184,690–4,511,229</td><td>SN</td><td>—</td><td>EA/EB</td><td>—</td></tr><tr><td rowspan="2">QC09a</td><td rowspan="2">C09</td><td rowspan="2">2,184,690–4,511,229</td><td rowspan="2">SY</td><td>Ver14_N<sup>+(EB)</sup>/Sel14_N<sup>-(EB)</sup> Md15_N<sup>+(EC)</sup>/Dij15_N<sup>-(EA)</sup></td><td rowspan="2">EA/EB/EC</td><td rowspan="2">—</td></tr><tr><td>Yeb15_N<sup>+(EA)</sup>/Ch14_N<sup>+(EB)</sup> LR15_N<sup>+(EC)</sup></td></tr><tr><td>QC09b</td><td>C09</td><td>51,377,310</td><td>SY</td><td>Ver14_N<sup>-(EC)</sup></td><td>—</td><td>—</td></tr></tbody></table>

A QTL is defined as a region that can gather different individual QTLs detected at different scales (MET, enviotype, environment; Supplemental Data 4) with overlapping positions. The range of the individual positions observed for a dedicated QTL is indicated in the “Positions” column. The positions refer to the reference genome Darmor-bzh v10 version (Rousseau-Gueutin et al 2020). The “Traits” column refers to Seed Yield (SY) or Seed Number (SN), the environments, envirotypes and MET-22 indicated the different scales for which the considered QTL was detected

(2), EB (2), EC (1) and ED (1) (Table 5, Supplemental Data 4). The same analysis carried out for SN revealed that locus QA09 was also detected, but only in enviotype EA. Thus, QA09 was qualified as a stable genomic region controlling SY and SY components across a wide range of environmental conditions.

## Two stable QTL (QA09 and QC09a) and an interactive QTL (QA07a) were revealed using a multi-scale QTL detection (environment, enviotype, MET)

Finally, the QTL analysis revealed two loci (QA09 and QC09a) that can be characterized as stable and a locus (QA07a) that can be characterized as interactive. QA09 was detected across the MET-22 in six environments and two envirotypes, and QC09a was detected in three envirotypes and seven individual environments. QA07a was specifically detected for the most stressed enviotype (ED) and could be a good candidate for breeding programs. The allelic diversity underlying this QTL must be further investigated to evaluate the genetic diversity available at this locus for breeding purposes. Indeed, the 46 elite lines of the panel already fixed the favorable allele for the two stable QTL QA09 and QC09a, but not for the interactive one QA07a (Supp Data 4). Ten other

QTL were detected at the environment scale but could not be linked to a limiting factor profile.

To characterize the three main QTL detected, a mixed linear model was fitted to estimate the QTL×E interaction using the data obtained for the whole MET-22 and both for SY and SN. (Table 6). The QA09×E interaction explained 31.9 and 34.5% of the total variance for SY and SN, respectively, whereas the genotype effect explained only 13.2 and 20.1% of the total variance for SY and SN, respectively. The QC09a×E interaction explained 31.9 and 38.5% of the total variance for SY and SN, respectively, whereas the genotype effect explained only 12.3 and 18.7% of the total variance for SY and SN, respectively. For the QA07a, that was qualified as interactive, the QTL×E effect explained less than 1% of variance for both SY and SN, this can be related to the overall small effect of this QTL at the MET-22 scale.

## Diversity analysis at QA07a revealed a lack of favorable allele in the modern-grown varieties

The composition of the panel P173 (GSL+ vs GSL- lines) allowed highlighting the recent history of winter oilseed rape breeding, including rapid selection for low glucosinolate contents in seeds, and estimating its impact on genetic
<table><caption>**Table 5** Description of the QTL detected for seed yield and seed number stability/instability</caption><thead><tr><th>Name</th><th>Chromosome</th><th>Positions</th><th>Traits</th><th>Environments</th><th>Envirotypes</th><th>MET-22</th></tr></thead><tbody><tr><td>QA03</td><td>A03</td><td>6,452,989</td><td>SY</td><td>Pre15_N<sup>+(EB)</sup></td><td>–</td><td>–</td></tr><tr><td>QA05a</td><td>A05</td><td>2,078,385–2,080,519</td><td>SY</td><td>Chr14_N<sup>+(EB)</sup></td><td>–</td><td>–</td></tr><tr><td>QA05b</td><td>A05</td><td>3,630,228</td><td>SN</td><td>Pre14_N<sup>+(ED)</sup></td><td>–</td><td>–</td></tr><tr><td rowspan="2">QA07a</td><td rowspan="2">A07</td><td>138,830–1,009,271</td><td>SN</td><td>LR15_N<sup>+(EC)</sup>/Pre14_N<sup>-(ED)</sup></td><td rowspan="2">ED</td><td rowspan="2">–</td></tr><tr><td>138,830–1,009,271</td><td>SY</td><td>LR15_N<sup>+(EC)</sup></td></tr><tr><td>QA07b</td><td>A07</td><td>20,033,023–20,033,025</td><td>SY</td><td>LR15_N<sup>+(EC)</sup></td><td>–</td><td>–</td></tr><tr><td>QA09</td><td>A09</td><td>3,933,991–4,504,817</td><td>SY/SN</td><td>Yeb15_N<sup>+(EA)</sup>/Ver14_N<sup>+(EB)</sup> Ch14_N<sup>-(ED)</sup>/Ch14_N<sup>+(EB)</sup> Dij15_N<sup>-(EA)</sup>/LR15_N<sup>+(EC)</sup></td><td>EA/ED</td><td>Yes</td></tr><tr><td>QA10</td><td>A10</td><td>2,515,307</td><td>SN</td><td>Pre14_N<sup>-(ED)</sup>/Ver14_N<sup>-(EC)</sup></td><td>–</td><td>–</td></tr><tr><td>QC02</td><td>C02</td><td>60,063,320–60,615,209</td><td>SY</td><td>Ver14_N<sup>+(EB)</sup>/Yeb15_N<sup>+(EA)</sup></td><td>–</td><td>–</td></tr><tr><td>QC03a</td><td>C03</td><td>20,140,168–20,809,663</td><td>SY</td><td>Sel15_N<sup>+(EB)</sup></td><td>–</td><td>–</td></tr><tr><td>QC03b</td><td>C03</td><td>23,645,808</td><td>SN</td><td>–</td><td>EC</td><td>–</td></tr><tr><td>QC04</td><td>C04</td><td>6,033,108–6,289,953</td><td>SY</td><td>Dij15_N<sup>+(EA)</sup>/Pre15_N<sup>-(EB)</sup></td><td>–</td><td>–</td></tr><tr><td>QC08</td><td>C08</td><td>37,326,511–37,326,814</td><td>SY</td><td>Ver14_N<sup>+(EB)</sup></td><td>–</td><td>–</td></tr><tr><td>QC09a</td><td>C09</td><td>2,184,690–4,511,229</td><td>SN</td><td>–</td><td>EA/EB</td><td>–</td></tr><tr><td rowspan="2">QC09a</td><td rowspan="2">C09</td><td rowspan="2">2,184,690–4,511,229</td><td rowspan="2">SN</td><td>Ver14_N<sup>+(EB)</sup>/Sel14_N<sup>-(EB)</sup> Md15_N<sup>+(EC)</sup>/Dij15_N<sup>-(EA)</sup></td><td rowspan="2">EA/EB/EC</td><td rowspan="2">–</td></tr><tr><td>Yeb15_N<sup>+(EA)</sup>/Ch14_N<sup>+(EB)</sup> LR15_N<sup>+(EC)</sup></td></tr><tr><td>QC09b</td><td>C09</td><td>51,377,310</td><td>SY</td><td>Ver14_N<sup>-(EC)</sup></td><td>–</td><td>–</td></tr></tbody></table>

A QTL is defined as a region that can gather different individual QTLs detected at different scales (MET, enviotype, environment; Supplemental Data 4) with overlapping positions. The range of the individual positions observed for a dedicated QTL is indicated in the “Positions” column. The positions refer to the reference genome Darmor-bzh v10 version (Rousseau-Gueutin et al 2020). The “Traits” column refers to Seed Yield (SY) or Seed Number (SN), the environments, envirotypes and MET-22 indicated the different scales for which the considered QTL was detected

(2), EB (2), EC (1) and ED (1) (Table 5, Supplemental Data 4). The same analysis carried out for SN revealed that locus QA09 was also detected, but only in enviotype EA. Thus, QA09 was qualified as a stable genomic region controlling SY and SY components across a wide range of environmental conditions.

## Two stable QTL (QA09 and QC09a) and an interactive QTL (QA07a) were revealed using a multi-scale QTL detection (environment, enviotype, MET)

Finally, the QTL analysis revealed two loci (QA09 and QC09a) that can be characterized as stable and a locus (QA07a) that can be characterized as interactive. QA09 was detected across the MET-22 in six environments and two envirotypes, and QC09a was detected in three envirotypes and seven individual environments. QA07a was specifically detected for the most stressed enviotype (ED) and could be a good candidate for breeding programs. The allelic diversity underlying this QTL must be further investigated to evaluate the genetic diversity available at this locus for breeding purposes. Indeed, the 46 elite lines of the panel already fixed the favorable allele for the two stable QTL QA09 and QC09a, but not for the interactive one QA07a (Supp Data 4). Ten other

QTL were detected at the environment scale but could not be linked to a limiting factor profile.

To characterize the three main QTL detected, a mixed linear model was fitted to estimate the QTL×E interaction using the data obtained for the whole MET-22 and both for SY and SN. (Table 6). The QA09×E interaction explained 31.9 and 34.5% of the total variance for SY and SN, respectively, whereas the genotype effect explained only 13.2 and 20.1% of the total variance for SY and SN, respectively. The QC09a×E interaction explained 31.9 and 38.5% of the total variance for SY and SN, respectively, whereas the genotype effect explained only 12.3 and 18.7% of the total variance for SY and SN, respectively. For the QA07a, that was qualified as interactive, the QTL×E effect explained less than 1% of variance for both SY and SN, this can be related to the overall small effect of this QTL at the MET-22 scale.

## Diversity analysis at QA07a revealed a lack of favorable allele in the modern-grown varieties

The composition of the panel P173 (GSL+ vs GSL- lines) allowed highlighting the recent history of winter oilseed rape breeding, including rapid selection for low glucosinolate contents in seeds, and estimating its impact on genetic
<table><caption>**Table 6** Results of the mixed linear model applied at the MET-22 scale to test the effect of the QTL×E or the QTL×En interaction</caption><thead><tr><th rowspan="2"></th><th colspan="3">Seed yield (SY)</th><th colspan="3">Seed number (SN)</th></tr><tr><th>QA09 (%)</th><th>QC09a (%)</th><th>QA07a (%)</th><th>QA09 (%)</th><th>QC09a (%)</th><th>QA07a (%)</th></tr></thead><tbody><tr><td>Genotype</td><td>13.2</td><td>12.3</td><td>24.3</td><td>20.1</td><td>18.7</td><td>31.9</td></tr><tr><td>Trial</td><td>35.7</td><td>35.7</td><td>56</td><td>19.4</td><td>16.7</td><td>41.4</td></tr><tr><td>Block×Trial</td><td>4.1</td><td>4.2</td><td>3.7</td><td>5.4</td><td>5.3</td><td>4.9</td></tr><tr><td>QTL×Trial</td><td>31.9</td><td>31.9</td><td>0.3</td><td>34.5</td><td>38.5</td><td>0.9</td></tr><tr><td>residual</td><td>15.2</td><td>15.9</td><td>15.6</td><td>20.7</td><td>20.8</td><td>20.9</td></tr></tbody></table>

The three QTL QA09, QC09a and QA07a (described in Table 5) were considered. The following mixed Linear model was used $Y_{ijkl} = \mu + G_i + E_k + marker_j \times E_k + R_{l(k)} + \epsilon_{ijkl}$, with the genotype (G), the environment (E), the QTL (marker), the genotype by environment interaction effect (G×E) and the block effect (R) nested in the trial., as proposed by Happ et al. (2021). The QTL effect is tested using the marker presenting the highest -log(p-value) within its confidence interval. In each case, the percentage of variance explained by the considered effect is indicated

<table><caption>**Table 7** Mean nucleotide diversity (π) and mean F<sub>ST</sub> statistics per chromosome and per population</caption><thead><tr><th rowspan="2"></th><th colspan="3">Mean π</th><th rowspan="2">Mean F<sub>ST</sub><br>GSL+ vs GSL-</th></tr><tr><th>P173</th><th>GSL+</th><th>GSL-</th></tr></thead><tbody><tr><td>Whole genome</td><td>1.69e<sup>-04</sup></td><td>1.77e<sup>-04</sup></td><td>1.61e<sup>-04</sup></td><td>0.078</td></tr><tr><td>Sub-genome A</td><td>1.81e<sup>-04</sup></td><td>1.92e<sup>-04</sup></td><td>1.70e<sup>-04</sup></td><td>0.070</td></tr><tr><td>Sub-genome C</td><td>1.58e<sup>-04</sup></td><td>1.62e<sup>-04</sup></td><td>1.51e<sup>-04</sup></td><td>0.086</td></tr><tr><td>A01</td><td>1.34e<sup>-04</sup></td><td>1.63e<sup>-04</sup></td><td>1.19e<sup>-04</sup></td><td>0.058</td></tr><tr><td>A02</td><td>1.55e<sup>-04</sup></td><td>1.67e<sup>-04</sup></td><td>1.47e<sup>-04</sup></td><td>0.055</td></tr><tr><td>A03</td><td>2.03e<sup>-04</sup></td><td>2.26e<sup>-04</sup></td><td>1.89e<sup>-04</sup></td><td>0.057</td></tr><tr><td>A04</td><td>1.65e<sup>-04</sup></td><td>1.69e<sup>-04</sup></td><td>1.60e<sup>-04</sup></td><td>0.060</td></tr><tr><td>A05</td><td>1.99e<sup>-04</sup></td><td>2.06e<sup>-04</sup></td><td>1.88e<sup>-04</sup></td><td>0.064</td></tr><tr><td>A06</td><td>2.17e<sup>-04</sup></td><td>2.28e<sup>-04</sup></td><td>2.06e<sup>-04</sup></td><td>0.052</td></tr><tr><td>A07</td><td>2.03e<sup>-04</sup></td><td>2.10e<sup>-04</sup></td><td>1.98e<sup>-04</sup></td><td>0.046</td></tr><tr><td>A08</td><td>1.96e<sup>-04</sup></td><td>1.91e<sup>-04</sup></td><td>1.81e<sup>-04</sup></td><td>0.159</td></tr><tr><td>A09</td><td>1.45e<sup>-04</sup></td><td>1.46e<sup>-04</sup></td><td>1.30e<sup>-04</sup></td><td>0.128</td></tr><tr><td>A10</td><td>1.76e<sup>-04</sup></td><td>1.86e<sup>-04</sup></td><td>1.69e<sup>-04</sup></td><td>0.039</td></tr><tr><td>C01</td><td>1.73e<sup>-04</sup></td><td>1.22e<sup>-04</sup></td><td>1.85e<sup>-04</sup></td><td>0.056</td></tr><tr><td>C02</td><td>1.31e<sup>-04</sup></td><td>1.46e<sup>-04</sup></td><td>1.21e<sup>-04</sup></td><td>0.111</td></tr><tr><td>C03</td><td>1.85e<sup>-04</sup></td><td>1.64e<sup>-04</sup></td><td>1.85e<sup>-04</sup></td><td>0.107</td></tr><tr><td>C04</td><td>1.43e<sup>-04</sup></td><td>1.82e<sup>-04</sup></td><td>1.23e<sup>-04</sup></td><td>0.087</td></tr><tr><td>C05</td><td>1.40e<sup>-04</sup></td><td>1.68e<sup>-04</sup></td><td>1.33e<sup>-04</sup></td><td>0.045</td></tr><tr><td>C06</td><td>1.60e<sup>-04</sup></td><td>1.85e<sup>-04</sup></td><td>1.42e<sup>-04</sup></td><td>0.098</td></tr><tr><td>C07</td><td>1.75e<sup>-04</sup></td><td>1.68e<sup>-04</sup></td><td>1.70e<sup>-04</sup></td><td>0.075</td></tr><tr><td>C08</td><td>1.59e<sup>-04</sup></td><td>1.77e<sup>-04</sup></td><td>1.54e<sup>-04</sup></td><td>0.046</td></tr><tr><td>C09</td><td>1.29e<sup>-04</sup></td><td>1.49e<sup>-04</sup></td><td>1.16e<sup>-04</sup></td><td>0.125</td></tr></tbody></table>

P173 corresponds to the whole diversity set (173 accessions). GSL+ corresponds to the 46 accessions of the P173 with high contents of glucosinolate in seeds (> 18 µmol g<sup>-1</sup>). GSL- corresponds to the 127 accessions of the P173 with low contents of glucosinolate in seeds (< 18 µmol g<sup>-1</sup>)

diversity at the whole genome scale, as well as at the scale of previously detected QTL. The F<sub>ST</sub> analysis showed a slight differentiation between the GSL+ and GSL- populations (F<sub>ST</sub> value of 0.078 at the whole genome scale, Table 7).

However, the genetic differentiation was higher when considering specific chromosomes such as A09 and C09, already known to carry genes controlling seed glucosinolate pathway, but also for chromosomes A08, C02 and C03 (F<sub>ST</sub> values > 0.1) (Table 7). No specific pattern was observed on the A07 chromosome harboring the interactive QA07a QTL for SN.

The nucleotide diversity index π was calculated at the whole genome scale and for each chromosome for both panel GSL+ and GSL-. Special attention was given to GSL- since this germplasm is more connected to the elite germplasm currently used by breeders. A scan of π values was also carried out targeting detected QTL regions. On average, lower nucleotide diversity was observed in the GSL+ lines than in the GSL- lines excepted for chromosomes A07, A09, A10, C01, C07 and C08 (Table 7). However, the study of the π index at QTL QA07a showed a higher diversity in the GSL+ lines than in the GSL- (Fig. 4). These results showed a higher nucleotide diversity in GSL+ and consequently a lack of diversity in the GSL- germplasm, corresponding here to a deficit of the favorable allele.. In contrast, at the QTL QA09, a higher π value was observed in the GSL- illustrating a gain of nucleotide diversity induced by breeding at this locus (Fig. 4).

## Discussion

Based on a comprehensive environmental characterization of a MET composed of 22 locations, we identified a group of environments (ED) characterized by a combination of stresses (poor vernalization conditions, N stress and low temperature and radiation during flowering and grain filling period) that drastically impacted seed yield. In addition, we identified a QTL specific to these conditions (QA07a). This QTL was characterized by a reduced genetic diversity in the modern germplasm, when compared to older cultivars. Our analysis also highlighted two stable QTL controlling seed
<table><caption>**Table 6** Results of the mixed linear model applied at the MET-22 scale to test the effect of the QTL×E or the QTL×En interaction</caption><thead><tr><th rowspan="2"></th><th colspan="3">Seed yield (SY)</th><th colspan="3">Seed number (SN)</th></tr><tr><th>QA09 (%)</th><th>QC09a (%)</th><th>QA07a (%)</th><th>QA09 (%)</th><th>QC09a (%)</th><th>QA07a (%)</th></tr></thead><tbody><tr><td>Genotype</td><td>13.2</td><td>12.3</td><td>24.3</td><td>20.1</td><td>18.7</td><td>31.9</td></tr><tr><td>Trial</td><td>35.7</td><td>35.7</td><td>56</td><td>19.4</td><td>16.7</td><td>41.4</td></tr><tr><td>Block×Trial</td><td>4.1</td><td>4.2</td><td>3.7</td><td>5.4</td><td>5.3</td><td>4.9</td></tr><tr><td>QTL×Trial</td><td>31.9</td><td>31.9</td><td>0.3</td><td>34.5</td><td>38.5</td><td>0.9</td></tr><tr><td>residual</td><td>15.2</td><td>15.9</td><td>15.6</td><td>20.7</td><td>20.8</td><td>20.9</td></tr></tbody></table>

The three QTL QA09, QC09a and QA07a (described in Table 5) were considered. The following mixed Linear model was used $Y_{ijkl} = \mu + G_i + E_k + marker_j \times E_k + R_{l(k)} + \epsilon_{ijkl}$, with the genotype (G), the environment (E), the QTL (marker), the genotype by environment interaction effect (G×E) and the block effect (R) nested in the trial., as proposed by Happ et al. (2021). The QTL effect is tested using the marker presenting the highest -log(p-value) within its confidence interval. In each case, the percentage of variance explained by the considered effect is indicated

<table><caption>**Table 7** Mean nucleotide diversity (π) and mean F<sub>ST</sub> statistics per chromosome and per population</caption><thead><tr><th rowspan="2"></th><th colspan="3">Mean π</th><th rowspan="2">Mean F<sub>ST</sub><br>GSL+ vs GSL-</th></tr><tr><th>P173</th><th>GSL+</th><th>GSL-</th></tr></thead><tbody><tr><td>Whole genome</td><td>1.69e<sup>-04</sup></td><td>1.77e<sup>-04</sup></td><td>1.61e<sup>-04</sup></td><td>0.078</td></tr><tr><td>Sub-genome A</td><td>1.81e<sup>-04</sup></td><td>1.92e<sup>-04</sup></td><td>1.70e<sup>-04</sup></td><td>0.070</td></tr><tr><td>Sub-genome C</td><td>1.58e<sup>-04</sup></td><td>1.62e<sup>-04</sup></td><td>1.51e<sup>-04</sup></td><td>0.086</td></tr><tr><td>A01</td><td>1.34e<sup>-04</sup></td><td>1.63e<sup>-04</sup></td><td>1.19e<sup>-04</sup></td><td>0.058</td></tr><tr><td>A02</td><td>1.55e<sup>-04</sup></td><td>1.67e<sup>-04</sup></td><td>1.47e<sup>-04</sup></td><td>0.055</td></tr><tr><td>A03</td><td>2.03e<sup>-04</sup></td><td>2.26e<sup>-04</sup></td><td>1.89e<sup>-04</sup></td><td>0.057</td></tr><tr><td>A04</td><td>1.65e<sup>-04</sup></td><td>1.69e<sup>-04</sup></td><td>1.60e<sup>-04</sup></td><td>0.060</td></tr><tr><td>A05</td><td>1.99e<sup>-04</sup></td><td>2.06e<sup>-04</sup></td><td>1.88e<sup>-04</sup></td><td>0.064</td></tr><tr><td>A06</td><td>2.17e<sup>-04</sup></td><td>2.28e<sup>-04</sup></td><td>2.06e<sup>-04</sup></td><td>0.052</td></tr><tr><td>A07</td><td>2.03e<sup>-04</sup></td><td>2.10e<sup>-04</sup></td><td>1.98e<sup>-04</sup></td><td>0.046</td></tr><tr><td>A08</td><td>1.96e<sup>-04</sup></td><td>1.91e<sup>-04</sup></td><td>1.81e<sup>-04</sup></td><td>0.159</td></tr><tr><td>A09</td><td>1.45e<sup>-04</sup></td><td>1.46e<sup>-04</sup></td><td>1.30e<sup>-04</sup></td><td>0.128</td></tr><tr><td>A10</td><td>1.76e<sup>-04</sup></td><td>1.86e<sup>-04</sup></td><td>1.69e<sup>-04</sup></td><td>0.039</td></tr><tr><td>C01</td><td>1.73e<sup>-04</sup></td><td>1.22e<sup>-04</sup></td><td>1.85e<sup>-04</sup></td><td>0.056</td></tr><tr><td>C02</td><td>1.31e<sup>-04</sup></td><td>1.46e<sup>-04</sup></td><td>1.21e<sup>-04</sup></td><td>0.111</td></tr><tr><td>C03</td><td>1.85e<sup>-04</sup></td><td>1.64e<sup>-04</sup></td><td>1.85e<sup>-04</sup></td><td>0.107</td></tr><tr><td>C04</td><td>1.43e<sup>-04</sup></td><td>1.82e<sup>-04</sup></td><td>1.23e<sup>-04</sup></td><td>0.087</td></tr><tr><td>C05</td><td>1.40e<sup>-04</sup></td><td>1.68e<sup>-04</sup></td><td>1.33e<sup>-04</sup></td><td>0.045</td></tr><tr><td>C06</td><td>1.60e<sup>-04</sup></td><td>1.85e<sup>-04</sup></td><td>1.42e<sup>-04</sup></td><td>0.098</td></tr><tr><td>C07</td><td>1.75e<sup>-04</sup></td><td>1.68e<sup>-04</sup></td><td>1.70e<sup>-04</sup></td><td>0.075</td></tr><tr><td>C08</td><td>1.59e<sup>-04</sup></td><td>1.77e<sup>-04</sup></td><td>1.54e<sup>-04</sup></td><td>0.046</td></tr><tr><td>C09</td><td>1.29e<sup>-04</sup></td><td>1.49e<sup>-04</sup></td><td>1.16e<sup>-04</sup></td><td>0.125</td></tr></tbody></table>

P173 corresponds to the whole diversity set (173 accessions). GSL+ corresponds to the 46 accessions of the P173 with high contents of glucosinolate in seeds (> 18 µmol g<sup>-1</sup>). GSL- corresponds to the 127 accessions of the P173 with low contents of glucosinolate in seeds (< 18 µmol g<sup>-1</sup>)

diversity at the whole genome scale, as well as at the scale of previously detected QTL. The F<sub>ST</sub> analysis showed a slight differentiation between the GSL+ and GSL- populations (F<sub>ST</sub> value of 0.078 at the whole genome scale, Table 7).

However, the genetic differentiation was higher when considering specific chromosomes such as A09 and C09, already known to carry genes controlling seed glucosinolate pathway, but also for chromosomes A08, C02 and C03 (F<sub>ST</sub> values > 0.1) (Table 7). No specific pattern was observed on the A07 chromosome harboring the interactive QA07a QTL for SN.

The nucleotide diversity index π was calculated at the whole genome scale and for each chromosome for both panel GSL+ and GSL-. Special attention was given to GSL- since this germplasm is more connected to the elite germplasm currently used by breeders. A scan of π values was also carried out targeting detected QTL regions. On average, lower nucleotide diversity was observed in the GSL+ lines than in the GSL- lines excepted for chromosomes A07, A09, A10, C01, C07 and C08 (Table 7). However, the study of the π index at QTL QA07a showed a higher diversity in the GSL+ lines than in the GSL- (Fig. 4). These results showed a higher nucleotide diversity in GSL+ and consequently a lack of diversity in the GSL- germplasm, corresponding here to a deficit of the favorable allele.. In contrast, at the QTL QA09, a higher π value was observed in the GSL- illustrating a gain of nucleotide diversity induced by breeding at this locus (Fig. 4).

## Discussion

Based on a comprehensive environmental characterization of a MET composed of 22 locations, we identified a group of environments (ED) characterized by a combination of stresses (poor vernalization conditions, N stress and low temperature and radiation during flowering and grain filling period) that drastically impacted seed yield. In addition, we identified a QTL specific to these conditions (QA07a). This QTL was characterized by a reduced genetic diversity in the modern germplasm, when compared to older cultivars. Our analysis also highlighted two stable QTL controlling seed
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Nucleotide diversity ($\pi$) under QTL QA07a and QA09. The nucleotide diversity was calculated for all P173 population (blue curve), the GSL+ (yellow curve) and GSL- (pink curve) accessions. The associated SNP for each QTL are indicated in black. The $\pi$ index is calculated for windows of 10 Kbp (color figure online)</figcaption></figure>

yields (QA09 and QC09a) that were expressed whatever the environmental conditions, as well as ten other QTL specifically detected in one or two single environments of the MET.

The question of crop adaptation and the underlying genetic determinism of crop stability is clearly a major challenge for agriculture for the coming decades. To answer this question, tools are being developed and the number of related studies raised drastically during the last years (van Eeuwijk et al. 2010, 2016). However, to our knowledge, only few studies dedicated to the analysis of *Brassica napus* genetic determinism of yield stability have been reported yet. The recent studies of yield stability specifically targeted the response of seed yield to water deficit (Zandberg et al. 2022; Raman et al. 2023) but did not address combination of different agro-pedoclimatic limiting factors. To address yield stability, authors usually test the QTL detected across different environmental conditions (Li et al. 2016; Lu et al. 2017; Zou et al. 2022). However, the proposed experimental designs are often limited to 4–6 environments, and $G \times E$ and stability genetic determinism are addressed by the confrontation of QTL detected for single environment to QTL detected across all the environments (Wang et al. 2018; Zheng et al. 2017; Lu et al. 2017; Li et al. 2020; Arifuzzaman et al. 2019; Pal et al. 2021; Sun et al. 2016; Gajardo et al. 2015) or by the meta-analysis of QTL detected for each environment using linkage analyses (Xie et al. 2020; Deng et al. 2019).

Within the MET-22, 25% of SY variation was due to a genetic effect, 55% to an environmental effect and 6% to $G \times E$. This predominance of the E and G effects over the

$G \times E$ effect supports the choice of grouping environments according to SY-limiting factors. This method made it possible to explain the E effect using envirotypes, and to examine the genetic determinant of SY in relation to environment characteristics. Indeed, the agro-pedoclimatic description of the MET-22 environments allowed explaining up to 70% of the environmental effect and only 24.6% of the $G \times E$ effect affecting seed yield. At the envirotype scale, the $G \times E$ part was reduced in the non-limiting envirotype EA (2.3%), but not in the most stressed envirotype ED (15.9%), when compared to MET-22 (5.8%). This result may be linked to the fact that the limiting factors used to group MET-22 were detected at the MET-47 scale. Indeed, although MET-22 and MET-47 present similar profiles in terms of limiting factors, slight differences were observed with regards to TMN-CW, $N_{total}$ and DAI indicators. From an analytical point of view, this resulted in a loss of power in the decomposition of the $G \times E$ effect at the MET-22 scale. However, from an agro-nomic point of view, this resulted in a more representative list of limiting factors at the MET scale representing Western European rapeseed growing areas, leading to envirotypes more representative of growing conditions.

Overall, the envirotyping-based methodology allowed prioritizing three QTLs involved in the genetic determinism of seed yield stability across the 16 QTL detected. These three QTLs have also been reported by Bouchet et al. (2016) genetic germplasm related to the P173 population used here. The stable QTL QA09 was also reported in different studies for seed yield-related traits as summed up by Raboanatahiry
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Nucleotide diversity ($\pi$) under QTL QA07a and QA09. The nucleotide diversity was calculated for all P173 population (blue curve), the GSL+ (yellow curve) and GSL- (pink curve) accessions. The associated SNP for each QTL are indicated in black. The $\pi$ index is calculated for windows of 10 Kbp (color figure online)</figcaption></figure>

yields (QA09 and QC09a) that were expressed whatever the environmental conditions, as well as ten other QTL specifically detected in one or two single environments of the MET.

The question of crop adaptation and the underlying genetic determinism of crop stability is clearly a major challenge for agriculture for the coming decades. To answer this question, tools are being developed and the number of related studies raised drastically during the last years (van Eeuwijk et al. 2010, 2016). However, to our knowledge, only few studies dedicated to the analysis of *Brassica napus* genetic determinism of yield stability have been reported yet. The recent studies of yield stability specifically targeted the response of seed yield to water deficit (Zandberg et al. 2022; Raman et al. 2023) but did not address combination of different agro-pedoclimatic limiting factors. To address yield stability, authors usually test the QTL detected across different environmental conditions (Li et al. 2016; Lu et al. 2017; Zou et al. 2022). However, the proposed experimental designs are often limited to 4–6 environments, and $G \times E$ and stability genetic determinism are addressed by the confrontation of QTL detected for single environment to QTL detected across all the environments (Wang et al. 2018; Zheng et al. 2017; Lu et al. 2017; Li et al. 2020; Arifuzzaman et al. 2019; Pal et al. 2021; Sun et al. 2016; Gajardo et al. 2015) or by the meta-analysis of QTL detected for each environment using linkage analyses (Xie et al. 2020; Deng et al. 2019).

Within the MET-22, 25% of SY variation was due to a genetic effect, 55% to an environmental effect and 6% to $G \times E$. This predominance of the E and G effects over the

$G \times E$ effect supports the choice of grouping environments according to SY-limiting factors. This method made it possible to explain the E effect using envirotypes, and to examine the genetic determinant of SY in relation to environment characteristics. Indeed, the agro-pedoclimatic description of the MET-22 environments allowed explaining up to 70% of the environmental effect and only 24.6% of the $G \times E$ effect affecting seed yield. At the envirotype scale, the $G \times E$ part was reduced in the non-limiting envirotype EA (2.3%), but not in the most stressed envirotype ED (15.9%), when compared to MET-22 (5.8%). This result may be linked to the fact that the limiting factors used to group MET-22 were detected at the MET-47 scale. Indeed, although MET-22 and MET-47 present similar profiles in terms of limiting factors, slight differences were observed with regards to TMN-CW, $N_{total}$ and DAI indicators. From an analytical point of view, this resulted in a loss of power in the decomposition of the $G \times E$ effect at the MET-22 scale. However, from an agro-nomic point of view, this resulted in a more representative list of limiting factors at the MET scale representing Western European rapeseed growing areas, leading to envirotypes more representative of growing conditions.

Overall, the envirotyping-based methodology allowed prioritizing three QTLs involved in the genetic determinism of seed yield stability across the 16 QTL detected. These three QTLs have also been reported by Bouchet et al. (2016) genetic germplasm related to the P173 population used here. The stable QTL QA09 was also reported in different studies for seed yield-related traits as summed up by Raboanatahiry
et al. (2018), whereas no-colocalization was found in the literature for QC09a. The interactive QTL QA07a may correspond to the QTL detected previously for branching number (Zhao et al. 2016) or thousand seed weight and plant height (Quijada et al. 2006; Udall et al. 2006). The 13 remaining SY-related QTLs were only detected in one to two environments. This pattern of QTL specificity to a unique environment is widely reported in literature (Bouchet et al. 2016; Garin et al. 2020) and offers a first approach to QTL × E interactions, but the robustness of these specific QTL is also questionable. The methodology that we proposed in this study was helpful to prioritize QTL and link the QTL with agro-pedoclimatic scenarii.

The low level of coincidence between our results and the QTL reported in the literature can result from differences in terms of genetic material used for GWAS or linkage studies. Indeed, most of studies reporting QTL for SY-related traits were carried out using spring type germplasm experimented in the field under short growth cycle conditions, unlike the winter oilseed rape germplasm used for this study. Winter and spring germplasms went through separated breeding history that can explain the differences of QTL detected, driven through selection for different breeding targets. Thus, broadening the genetic diversity considered for GWAS may help identifying additional seed yield-related QTL. For example, adding semi-winter accessions, more likely to withstand western Europe growing conditions, is a promising way to increase the population resolution by reducing the extent of linkage disequilibrium. A second way to improve the power of GWAS consists in increasing the number of genomic markers. However, recent genomics advances already led to a common use of resolutive genotyping resources such as the 60K Illumina array in most of the genetic studies (Clarke et al. 2016). Here, we developed and used a novel and dense genomic resource based on whole genome exome capture to characterize the population (Leveugle et al. 2015). This resource provided 217,805 SNP covering the entire genome and therefore increased the resolution of the QTL regions as a higher number SNP were detected per genomic region. Due to this genomic resource, the genome can be decomposed into 6763 haploblocks that presented a mean size of 48kb and that were in average composed of 28 SNPs (Supp Data 5). We also could precisely describe and capture the LD pattern of the P173 WOSR panel (Supp Fig. 3). This consists in an important resource, especially when dealing with complex traits such as seed yield.

To facilitate interpretation of QTL stability, approaches have already been proposed based on clustering of environments within a MET for other species. This involves grouping individual environments into mega-environments such that genotypes exhibit similar behaviors in a mega-environment and may differ between mega-environments. Thus, Moreau et al. (2004) proposed a clustering of environments

according to the G×E interaction matrix and were able to explain the specificity of QTL according to climatic or water stress conditions. More recently, the use of factor analytic models has been proposed and is specifically devoted to GxE decomposition among MET (Smith et al 2021). This method is highly resolutive in explaining GxE and in building cluster that minimize within GxE. The use of this methodology could have allowed to detect more “interactive” QTL; it is to say QTL presenting higher QTLxE effect. However, from an agronomic point of view, it is more difficult to link the cluster characteristics to environmental stresses experienced by plants. Millet et al. (2016) and Touzy et al. (2019) proposed to cluster the environments of a MET according to drought and/or heat scenarii and then to study the pattern of QTL effects from one scenario to another. Most of SY-related QTL presented significant interaction with the climatic scenario and for some of them, the favorable allele changed according to the scenario (Millet et al. 2016). Our study supports the interest of clustering the environments to identify and interpret QTL effects. In this study, we identified QTL (QA07a, QC03b) present in given envirotypes (ED and EC, respectively) and absent for others, a pattern also revealed in maize by Millet et al. (2016) and Touzy et al. (2019) when opposing contrasting scenarii. However, as opposed to Millet et al. (2016), we did not record any change of the favorable allele at a given QTL, depending on the envirotype. These two last studies focused their environments clustering on a priori defined stresses (water and temperature) and developed a targeted characterization of the environments of the MET. For winter oilseed rape, the 11-month growth cycle makes it more difficult to focus on a dedicated stress, and we therefore choose an alternative method that reports a posteriori the main combinations of environmental factors that did impact seed yield. Even if nitrogen input was managed to be one of the main limiting factors, it was clearly shown that it had to be considered in combination with others stress (radiation, temperature, ...) occurring during the crop cycle, thus making ineffective an a priori clustering based on nitrogen indicators only. Moreover, Ravier et al. (2017) showed that N deficiencies, even intense, do not always affect seed yield, especially if they occur early during the wheat growth cycle. This method was also successfully used to identify limiting factors occurring over a MET for barley (Beillouin et al. 2018) and highlights the need for indicators to account for potential stresses occurring in the field.

In this study, the envirotype ED corresponded to the combination of limiting factors that impacted the most seed yield. Its limiting factors targeted different phases of the crop cycle: winter (fulfillment of vernalization requirements and high temperatures during winter), bolting (N fertilization), flowering (lack of solar radiation) and grain filling (lower temperature). This high impacting combination of stress was observed for the four environments of envirotype ED.
et al. (2018), whereas no-colocalization was found in the literature for QC09a. The interactive QTL QA07a may correspond to the QTL detected previously for branching number (Zhao et al. 2016) or thousand seed weight and plant height (Quijada et al. 2006; Udall et al. 2006). The 13 remaining SY-related QTLs were only detected in one to two environments. This pattern of QTL specificity to a unique environment is widely reported in literature (Bouchet et al. 2016; Garin et al. 2020) and offers a first approach to QTL × E interactions, but the robustness of these specific QTL is also questionable. The methodology that we proposed in this study was helpful to prioritize QTL and link the QTL with agro-pedoclimatic scenarii.

The low level of coincidence between our results and the QTL reported in the literature can result from differences in terms of genetic material used for GWAS or linkage studies. Indeed, most of studies reporting QTL for SY-related traits were carried out using spring type germplasm experimented in the field under short growth cycle conditions, unlike the winter oilseed rape germplasm used for this study. Winter and spring germplasms went through separated breeding history that can explain the differences of QTL detected, driven through selection for different breeding targets. Thus, broadening the genetic diversity considered for GWAS may help identifying additional seed yield-related QTL. For example, adding semi-winter accessions, more likely to withstand western Europe growing conditions, is a promising way to increase the population resolution by reducing the extent of linkage disequilibrium. A second way to improve the power of GWAS consists in increasing the number of genomic markers. However, recent genomics advances already led to a common use of resolutive genotyping resources such as the 60K Illumina array in most of the genetic studies (Clarke et al. 2016). Here, we developed and used a novel and dense genomic resource based on whole genome exome capture to characterize the population (Leveugle et al. 2015). This resource provided 217,805 SNP covering the entire genome and therefore increased the resolution of the QTL regions as a higher number SNP were detected per genomic region. Due to this genomic resource, the genome can be decomposed into 6763 haploblocks that presented a mean size of 48kb and that were in average composed of 28 SNPs (Supp Data 5). We also could precisely describe and capture the LD pattern of the P173 WOSR panel (Supp Fig. 3). This consists in an important resource, especially when dealing with complex traits such as seed yield.

To facilitate interpretation of QTL stability, approaches have already been proposed based on clustering of environments within a MET for other species. This involves grouping individual environments into mega-environments such that genotypes exhibit similar behaviors in a mega-environment and may differ between mega-environments. Thus, Moreau et al. (2004) proposed a clustering of environments

according to the G×E interaction matrix and were able to explain the specificity of QTL according to climatic or water stress conditions. More recently, the use of factor analytic models has been proposed and is specifically devoted to GxE decomposition among MET (Smith et al 2021). This method is highly resolutive in explaining GxE and in building cluster that minimize within GxE. The use of this methodology could have allowed to detect more “interactive” QTL; it is to say QTL presenting higher QTLxE effect. However, from an agronomic point of view, it is more difficult to link the cluster characteristics to environmental stresses experienced by plants. Millet et al. (2016) and Touzy et al. (2019) proposed to cluster the environments of a MET according to drought and/or heat scenarii and then to study the pattern of QTL effects from one scenario to another. Most of SY-related QTL presented significant interaction with the climatic scenario and for some of them, the favorable allele changed according to the scenario (Millet et al. 2016). Our study supports the interest of clustering the environments to identify and interpret QTL effects. In this study, we identified QTL (QA07a, QC03b) present in given envirotypes (ED and EC, respectively) and absent for others, a pattern also revealed in maize by Millet et al. (2016) and Touzy et al. (2019) when opposing contrasting scenarii. However, as opposed to Millet et al. (2016), we did not record any change of the favorable allele at a given QTL, depending on the envirotype. These two last studies focused their environments clustering on a priori defined stresses (water and temperature) and developed a targeted characterization of the environments of the MET. For winter oilseed rape, the 11-month growth cycle makes it more difficult to focus on a dedicated stress, and we therefore choose an alternative method that reports a posteriori the main combinations of environmental factors that did impact seed yield. Even if nitrogen input was managed to be one of the main limiting factors, it was clearly shown that it had to be considered in combination with others stress (radiation, temperature, ...) occurring during the crop cycle, thus making ineffective an a priori clustering based on nitrogen indicators only. Moreover, Ravier et al. (2017) showed that N deficiencies, even intense, do not always affect seed yield, especially if they occur early during the wheat growth cycle. This method was also successfully used to identify limiting factors occurring over a MET for barley (Beillouin et al. 2018) and highlights the need for indicators to account for potential stresses occurring in the field.

In this study, the envirotype ED corresponded to the combination of limiting factors that impacted the most seed yield. Its limiting factors targeted different phases of the crop cycle: winter (fulfillment of vernalization requirements and high temperatures during winter), bolting (N fertilization), flowering (lack of solar radiation) and grain filling (lower temperature). This high impacting combination of stress was observed for the four environments of envirotype ED.
Individually, each of these factors has already been shown to impact seed yield: thermal stress during flowering affects flower fertility, pod number and seed number (Morrison 1993; Angadi et al. 2000; Young et al. 2004); a poor verbalization conditions could lead to delayed or no flowering (Ferreira et al. 1995; Chandler et al. 2005) and N limitation is also known to impact seed yield (Rathke et al. 2006). Here, we were able to identify that these different stresses co-occurred in the field and truly impact seed yield. Excepting N stress, the main limiting factors in ED corresponded to climatic factors.

We suggest that combining envirotyping and QTL analysis must be considered an effective approach, enabling the identification of both agro-pedoclimatic indicators and QTL involved in yield stability. QTL with the highest effects and qualified as "stable" may nevertheless present significant QTL × E interactions. However, they are detected in most environments and envirotypes, with the favorable allele being the same whatever the environment considered, which underlines their interest in improving yield for a large range of environmental conditions. Specific QTLs such as QA07a presented smaller effect and were not significantly involved in QTL × E at the whole MET-22 scale. However, this QTL presented a specific interest for a specific combination of limiting agro-pedoclimatic conditions that was observed in the enviotype ED. The analyses performed at the MET-22 scale or directly at the single environment scale were not consistent to highlight this genomic region. The molecular diversity as this locus indicated that the GSL cultivars did not fix the favorable allele, demonstrating the interest of the GSL+ cultivars as valuable source of genetic diversity for improving seed yield and its stability. This was confirmed by the fact that the elite lines already fixed the favorable allele for QTL QA09 and QC09 but not for QTL QA07a (supp Data 4). It is to be noticed that the QA07a does not colocalize with the QTL for glucosinolate content on the A07 chromosome (Tang et al 2023). This QTL indeed colocalizes with the QA07b that was detected for SY only in LR15⁺. Therefore, at the QA07a region, the favorable allele is more frequent in the “++” germplasm but the underlying genes might not be related to the glucosinolate pathway.

For further breeding programs, there is a particular interest in validating the four envirotypes described in the present study for a wider range of agro-pedoclimatic conditions. Indeed, a posteriori analysis of larger climatic datasets at the same locations, coupled with crop physiological models, such as AZODYN-colza (Jeuffroy et al. 2006), may conduct to an estimation of the frequency of these four specific envirotypes across growing seasons. Envirotype ED can consist in a new target environment for breeding if it occurs rather frequently within the rapeseed production area. This strategy could also help redesigning multi-environment trials for WOSR breeding, for instance by discarding redundant

locations, reducing experimental costs and by maximizing the opportunities of desired envirotypes/pedoclimatic scenario within a MET. The envirotyping approach can also lead to an estimation of a similarity matrix of a MET locations, according to their limiting factors pattern. In silico experiments will be useful to test a wide range of genotypes and environmental conditions such as described by Wang et al. (2023) to validate our results. This can be a clue to identify accurate match between dedicated genotypes and environments (Resende et al. 2021) leading, notably, to better product placement for the seed industry.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00122-024-04664-3.

**Acknowledgements** The authors would like to thank the technical staff of the UMR IGEPP for collecting the field data as well as for management of samples at Le Rheu (Elise Alix, Bernard Moulin, Solenn Guichard, Alina Tollenaere, Tiffany Bourlet), as well as Cécile Baron who handled the genotyping data, and Mathieu Rousseau-Gueutin who provided genomic support concerning the Darmor-bzh v10 genome version. The authors would also like to thank the “Domaine de la Motte” Experimental Unit (INRAE Bretagne Normandie, Domaine de la Motte, 35650 Le Rheu) for the provision of the experimental plots, the cultural interventions and the agri-environmental data recorded and used for this study. The authors are also grateful to the partners of the RAPSODYN project (Innolea, Limagrain Europe, Lidea, MAS seeds, Syngenta, RAGT, Terres Inovia) that provided the field data at the MET scale.

**Author contribution statement** AL and NN conceived and designed the analyses, planned the experiments, collected the data, contributed to the interpretation of the results and supervised the work. EC analyzed the data, interpreted the results and wrote the manuscript with inputs from all co-authors. CS contributed to the data collection and ran the diversity analyses. ML ran the genomic work and analyzed the exome capture data. All authors discussed the results, contributed, edited and validated the final manuscript.

**Funding** This research was supported by two national collaborative projects untitled GENERGY (ANR-07-GPLA-016) funded by the French National Research Agency (ANR) and RAPSODYN (ANR-11-BTBR-0004) funded by the program "Investments for the Future".

**Data availability** The datasets generated and analyzed in this study are available using the following link (https://entrepot.recherche.data.gouv.fr/dataset.xhtml?persistentId=doi:10.57745/0TWHMA).

## Declarations

**Conflict of interest** The authors declare no conflicts of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will
Individually, each of these factors has already been shown to impact seed yield: thermal stress during flowering affects flower fertility, pod number and seed number (Morrison 1993; Angadi et al. 2000; Young et al. 2004); a poor verbalization conditions could lead to delayed or no flowering (Ferreira et al. 1995; Chandler et al. 2005) and N limitation is also known to impact seed yield (Rathke et al. 2006). Here, we were able to identify that these different stresses co-occurred in the field and truly impact seed yield. Excepting N stress, the main limiting factors in ED corresponded to climatic factors.

We suggest that combining envirotyping and QTL analysis must be considered an effective approach, enabling the identification of both agro-pedoclimatic indicators and QTL involved in yield stability. QTL with the highest effects and qualified as "stable" may nevertheless present significant QTL×E interactions. However, they are detected in most environments and envirotypes, with the favorable allele being the same whatever the environment considered, which underlines their interest in improving yield for a large range of environmental conditions. Specific QTLs such as QA07a presented smaller effect and were not significantly involved in QTL×E at the whole MET-22 scale. However, this QTL presented a specific interest for a specific combination of limiting agro-pedoclimatic conditions that was observed in the enviotype ED. The analyses performed at the MET-22 scale or directly at the single environment scale were not consistent to highlight this genomic region. The molecular diversity as this locus indicated that the GSL cultivars did not fix the favorable allele, demonstrating the interest of the GSL+ cultivars as valuable source of genetic diversity for improving seed yield and its stability. This was confirmed by the fact that the elite lines already fixed the favorable allele for QTL QA09 and QC09 but not for QTL QA07a (supp Data 4). It is to be noticed that the QA07a does not colocalize with the QTL for glucosinolate content on the A07 chromosome (Tang et al 2023). This QTL indeed colocalizes with the QA07b that was detected for SY only in LR15⁺. Therefore, at the QA07a region, the favorable allele is more frequent in the “++” germplasm but the underlying genes might not be related to the glucosinolate pathway.

For further breeding programs, there is a particular interest in validating the four envirotypes described in the present study for a wider range of agro-pedoclimatic conditions. Indeed, a posteriori analysis of larger climatic datasets at the same locations, coupled with crop physiological models, such as AZODYN-colza (Jeuffroy et al. 2006), may conduct to an estimation of the frequency of these four specific envirotypes across growing seasons. Envirotype ED can consist in a new target environment for breeding if it occurs rather frequently within the rapeseed production area. This strategy could also help redesigning multi-environment trials for WOSR breeding, for instance by discarding redundant

locations, reducing experimental costs and by maximizing the opportunities of desired envirotypes/pedoclimatic scenario within a MET. The envirotyping approach can also lead to an estimation of a similarity matrix of a MET locations, according to their limiting factors pattern. In silico experiments will be useful to test a wide range of genotypes and environmental conditions such as described by Wang et al. (2023) to validate our results. This can be a clue to identify accurate match between dedicated genotypes and environments (Resende et al. 2021) leading, notably, to better product placement for the seed industry.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00122-024-04664-3.

**Acknowledgements** The authors would like to thank the technical staff of the UMR IGEPP for collecting the field data as well as for management of samples at Le Rheu (Elise Alix, Bernard Moulin, Solenn Guichard, Alina Tollenaere, Tiffany Bourlet), as well as Cécile Baron who handled the genotyping data, and Mathieu Rousseau-Gueutin who provided genomic support concerning the Darmor-bzh v10 genome version. The authors would also like to thank the “Domaine de la Motte” Experimental Unit (INRAE Bretagne Normandie, Domaine de la Motte, 35650 Le Rheu) for the provision of the experimental plots, the cultural interventions and the agri-environmental data recorded and used for this study. The authors are also grateful to the partners of the RAPSODYN project (Innolea, Limagrain Europe, Lidea, MAS seeds, Syngenta, RAGT, Terres Inovia) that provided the field data at the MET scale.

**Author contribution statement** AL and NN conceived and designed the analyses, planned the experiments, collected the data, contributed to the interpretation of the results and supervised the work. EC analyzed the data, interpreted the results and wrote the manuscript with inputs from all co-authors. CS contributed to the data collection and ran the diversity analyses. ML ran the genomic work and analyzed the exome capture data. All authors discussed the results, contributed, edited and validated the final manuscript.

**Funding** This research was supported by two national collaborative projects untitled GENERGY (ANR-07-GPLA-016) funded by the French National Research Agency (ANR) and RAPSODYN (ANR-11-BTBR-0004) funded by the program "Investments for the Future".

**Data availability** The datasets generated and analyzed in this study are available using the following link (https://entrepot.recherche.data.gouv.fr/dataset.xhtml?persistentId=doi:10.57745/0TWHMA).

## Declarations

**Conflict of interest** The authors declare no conflicts of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will

<figure><img src="image_8.png" /></figure>
need to obtain permission directly from the copyright holder. To view a
copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References

Angadi SV, Cutforth HW, Miller PR, Mcconkey BG, Entz MH, Brandt SA (2000) Response of three *Brassica* species to high temperature stress during reproductive growth. Can J Plant Sci 80:693-702. https://doi.org/10.4141/P99-152

Arifuzzaman M, Oladzadabbasabadi A, McLean P, Rahman M (2019) Shovelomics for phenotyping root architectural traits of rapeseed/canola (*Brassica napus* L.) and genome-wide association mapping. Mol Genet Gen 294:985-1000. https://doi.org/10.1007/s00438-019-01563-x

Astle W, Balding DJ (2009) Population structure and cryptic relatedness in genetic association studies. Stat Sci 24:451-471. https://doi.org/10.1214/09-STS307

Bates D, Mächler M, Bolker B, Walker S (2015) Fitting linear mixed-effects models using lme4. J Stat Softw 67:1-48. https://doi.org/10.18637/jss.v067.i01

Becker HC, Leon J (1988) Stability analysis in plant breeding. Plant Breed 101:1-23. https://doi.org/10.1111/j.1439-0523.1988.tb00261.x

Beillouin D, jeffroy MH, Gauffreteau A (2018) Characterization of spatial and temporal combinations of climatic factors affecting yields: an empirical model applied to the French barley belt. Agric Meteorol 262:402-411. https://doi.org/10.1016/j.agrformet.2018.07.029

Bell JL, Sloan LC, Snyder MA (2004) Regional changes in extreme climatic events: a future climate scenario. J Climate 17:81-87. https://doi.org/10.1175/1520-0442(2004)017%3c0081:RCIECE%3e2.0.CO;2

Bouchet AS, Laperche A, Bissuel-Belaygue C, Baron C, Morice J, Rousseau-Gueutin M, Dheu JE, George P, Pinochet X, Foubert T, Maes O, Dugué D, Guinot F, Nesi N (2016) Genetic basis of nitrogen use efficiency and yield stability across environments in winter rapeseed. BMC Genet 17:131. https://doi.org/10.1186/s12863-016-0432-z

Browning BL, Browning SR (2009) A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals. Am J Hum Genet 84:210-223. https://doi.org/10.1016/j.ajhg.2009.01.005

Chandler J, Corbesier L, Spielmann P, Dettendorfer J, Stahl D, Apel K, Melzer S (2005) Modulating flowering time and prevention of pod shatter in oilseed rape. Mol Breed 15:87-94. https://doi.org/10.1007/s11032-004-2735-4

Chaves SFS, Evangelista JSPC, Trindade RS, Dias LAS, Guimarães PE, Guimarães LJM, Alves RS, Bhering LL, Dias KOG (2023) Employing factor analytic tools for selecting high-performance and stable tropical maize hybrids. Crop Sci 63:1114-1125. https://doi.org/10.1002/csc2.20911

Chidzanga C, Mullan D, Roy S, Baumann U, Garcia M (2022) Nested association mapping-based GWAS for grain yield and related traits in wheat grown under diverse Australian environments. Theor Appl Genet 135:4437-4456. https://doi.org/10.1007/s00122-022-04230-9

Clarke WE, Higgins EE, Plieske J, Wieseke R, Sidebottom C, Khedikar Y, Batley J, Edwards D, Meng J, Li R, Lawley CT, Pauquet J, Laga B, Cheung W, Iniguez-Luy F, Dyrszka E, Rae S, Stich B, Snowdon RJ, Sharpe AG, Ganal MW, Parkin IA (2016) A high-density SNP genotyping array for *Brassica napus* and its ancestral diploid species based on optimised selection of single-locus markers in

the allotetraploid genome. Theor Appl Genet 129:1887-1899. https://doi.org/10.1007/s00122-016-2746-7

Cooper M, Tang T, Gho C, Hart T, hammer G, Messina C, (2020) Integrating genetic gain and gap analysis to predict improvements in crop productivity. Crop Sci 60:582-604. https://doi.org/10.1002/csc2.20109

Corlouer E, Gauffreteau A, Bouchet AS, Bissuel-Belaygue C, Nesi N, Laperche A (2019) Envirotypes based on seed yield limiting factors allow to tackle G × E interactions. Agronomy 9:798. https://doi.org/10.3390/agronomy9120798

Cowling WA, Castro-Urrea FA, Stefanova KT, Li L, Banks RG, Saradadevi R, Sass O, Kinghorn BP, Siddique KHM (2023) Optimal contribution selection improves the rate of genetic gain in grain yield and yield stability in spring canola in Australia and Canada. Plants 12:383. https://doi.org/10.3390/plants12020383

Danecek P, Auton A, Abecasis G, Albers AC, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, 1000 Genomes Project Analysis Group (2011) The variant call format and VCFtools. Bioinformatics 27:2156-2158. https://doi.org/10.1093/bioinformatics/btr330

Delourme R, Laperche A, Bouchet AS, Jubault M, Paillard S, Manzanares-Dauleux MJ, Nesi N (2018) Genes and quantitative trait loci mapping for major agronomic traits in *Brassica napus* L. In: Liu S, Snowdon R, Chalhoub B (eds) The *Brassica napus* genome compendium of plant genomes. Springer, Cham. https://doi.org/10.1007/978-3-319-43694-4_3

Deng C, Liu H, Yao Y, Guo S, Xiao L, Fu Z, Du D (2019) QTL analysis of frou yield related-traits for *Brassica napus* L. in multiple environments. Mol Breed 39:166. https://doi.org/10.1007/s11032-019-1067-3

Diepenbrock W (2000) Yield analysis of winter oilseed rape (*Brassica napus* L.): a review. Field Crop Res 67:35-42. https://doi.org/10.1016/S0378-4290(00)00082-4

Diouf I, Derivot L, Koussevitzky S, Carretero Y, Bitton F, Moreau L, Causse M (2020) Genetic basis of phenotypic plasticity and genotype x environment interactions in a multi-parental tomato population. J Exp Bot 1:5365-5376. https://doi.org/10.1093/jxp/eraa265

dos Santos Silva PP, e Sousa MB, de Oliveira EJ, Morgante CV, Silva de Oliveira CR, Vieira SL, Borel JC (2021) Genome-wide association study of drought tolerance in cassava. Euphytica 217:60. https://doi.org/10.1007/s10681-021-02800-4

El-Soda M, Malosetti M, Zwaan BJ, Koornneef M, Aarts MGM (2014) Genotype × environment interaction QTL mapping in plants: lessons from Arabidopsis. Trends Plant Sci 19:390-398. https://doi.org/10.1016/j.tplants.2014.01.001

Ferreira ME, Satagopan J, Yandell BS, Williams PH, Osborn TC (1995) Mapping loci controlling vernalization requirement and flowering time in *Brassica napus*. Theor Appl Genet 90:727-732. https://doi.org/10.1007/BF00222140

Finlay KW, Wilkinson GN (1963) The analysis of adaptation in a plant-breeding programme. AustJ Agric Res 14:742-754. https://doi.org/10.1071/AR9630742

Gajardo HA, Wittkop B, Soto-Cerda B, Higgins EE, Parkin IAP, Snowdon RJ, Frederico ML, Iniguez-Luy FL (2015) Association mapping of seed quality traits in *Brassica napus* L. using GWAS and candidate QTL approaches. Mol Breed 35:143. https://doi.org/10.1007/s11032-015-0340-3

Gao X, Starmer J, Martin ER (2008) A multiple testing correction method for genetic association studies using correlated single nucleotide polymorphisms. Genet Epidemiol 32:361-369. https://doi.org/10.1002/gepi.20310

Garin V, Malosetti M, van Eeuwijk F (2020) Multi-parent multi-environment QTL analysis: an illustration with the EU-NAM Flint population. Theor Appl Genet 133:2627-2638. https://doi.org/10.1007/s00122-020-03621-0
need to obtain permission directly from the copyright holder. To view a
copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References

Angadi SV, Cutforth HW, Miller PR, Mcconkey BG, Entz MH, Brandt SA (2000) Response of three *Brassica* species to high temperature stress during reproductive growth. Can J Plant Sci 80:693-702. https://doi.org/10.4141/P99-152

Arifuzzaman M, Oladzadabbasabadi A, McLean P, Rahman M (2019) Shovelomics for phenotyping root architectural traits of rapeseed/canola (*Brassica napus* L.) and genome-wide association mapping. Mol Genet Gen 294:985-1000. https://doi.org/10.1007/s00438-019-01563-x

Astle W, Balding DJ (2009) Population structure and cryptic relatedness in genetic association studies. Stat Sci 24:451-471. https://doi.org/10.1214/09-STS307

Bates D, Mächler M, Bolker B, Walker S (2015) Fitting linear mixed-effects models using lme4. J Stat Softw 67:1-48. https://doi.org/10.18637/jss.v067.i01

Becker HC, Leon J (1988) Stability analysis in plant breeding. Plant Breed 101:1-23. https://doi.org/10.1111/j.1439-0523.1988.tb00261.x

Beillouin D, jeffroy MH, Gauffreteau A (2018) Characterization of spatial and temporal combinations of climatic factors affecting yields: an empirical model applied to the French barley belt. Agric Meteorol 262:402-411. https://doi.org/10.1016/j.agrformet.2018.07.029

Bell JL, Sloan LC, Snyder MA (2004) Regional changes in extreme climatic events: a future climate scenario. J Climate 17:81-87. https://doi.org/10.1175/1520-0442(2004)017%3c0081:RCIECE%3e2.0.CO;2

Bouchet AS, Laperche A, Bissuel-Belaygue C, Baron C, Morice J, Rousseau-Gueutin M, Dheu JE, George P, Pinochet X, Foubert T, Maes O, Dugué D, Guinot F, Nesi N (2016) Genetic basis of nitrogen use efficiency and yield stability across environments in winter rapeseed. BMC Genet 17:131. https://doi.org/10.1186/s12863-016-0432-z

Browning BL, Browning SR (2009) A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals. Am J Hum Genet 84:210-223. https://doi.org/10.1016/j.ajhg.2009.01.005

Chandler J, Corbesier L, Spielmann P, Dettendorfer J, Stahl D, Apel K, Melzer S (2005) Modulating flowering time and prevention of pod shatter in oilseed rape. Mol Breed 15:87-94. https://doi.org/10.1007/s11032-004-2735-4

Chaves SFS, Evangelista JSPC, Trindade RS, Dias LAS, Guimarães PE, Guimarães LJM, Alves RS, Bhering LL, Dias KOG (2023) Employing factor analytic tools for selecting high-performance and stable tropical maize hybrids. Crop Sci 63:1114-1125. https://doi.org/10.1002/csc2.20911

Chidzanga C, Mullan D, Roy S, Baumann U, Garcia M (2022) Nested association mapping-based GWAS for grain yield and related traits in wheat grown under diverse Australian environments. Theor Appl Genet 135:4437-4456. https://doi.org/10.1007/s00122-022-04230-9

Clarke WE, Higgins EE, Plieske J, Wieseke R, Sidebottom C, Khedikar Y, Batley J, Edwards D, Meng J, Li R, Lawley CT, Pauquet J, Laga B, Cheung W, Iniguez-Luy F, Dyrszka E, Rae S, Stich B, Snowdon RJ, Sharpe AG, Ganal MW, Parkin IA (2016) A high-density SNP genotyping array for *Brassica napus* and its ancestral diploid species based on optimised selection of single-locus markers in

the allotetraploid genome. Theor Appl Genet 129:1887-1899. https://doi.org/10.1007/s00122-016-2746-7

Cooper M, Tang T, Gho C, Hart T, hammer G, Messina C, (2020) Integrating genetic gain and gap analysis to predict improvements in crop productivity. Crop Sci 60:582-604. https://doi.org/10.1002/csc2.20109

Corlouer E, Gauffreteau A, Bouchet AS, Bissuel-Belaygue C, Nesi N, Laperche A (2019) Envirotypes based on seed yield limiting factors allow to tackle G × E interactions. Agronomy 9:798. https://doi.org/10.3390/agronomy9120798

Cowling WA, Castro-Urrea FA, Stefanova KT, Li L, Banks RG, Saradadevi R, Sass O, Kinghorn BP, Siddique KHM (2023) Optimal contribution selection improves the rate of genetic gain in grain yield and yield stability in spring canola in Australia and Canada. Plants 12:383. https://doi.org/10.3390/plants12020383

Danecek P, Auton A, Abecasis G, Albers AC, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, 1000 Genomes Project Analysis Group (2011) The variant call format and VCFtools. Bioinformatics 27:2156-2158. https://doi.org/10.1093/bioinformatics/btr330

Delourme R, Laperche A, Bouchet AS, Jubault M, Paillard S, Manzanares-Dauleux MJ, Nesi N (2018) Genes and quantitative trait loci mapping for major agronomic traits in *Brassica napus* L. In: Liu S, Snowdon R, Chalhoub B (eds) The *Brassica napus* genome compendium of plant genomes. Springer, Cham. https://doi.org/10.1007/978-3-319-43694-4_3

Deng C, Liu H, Yao Y, Guo S, Xiao L, Fu Z, Du D (2019) QTL analysis of frou yield related-traits for *Brassica napus* L. in multiple environments. Mol Breed 39:166. https://doi.org/10.1007/s11032-019-1067-3

Diepenbrock W (2000) Yield analysis of winter oilseed rape (*Brassica napus* L.): a review. Field Crop Res 67:35-42. https://doi.org/10.1016/S0378-4290(00)00082-4

Diouf I, Derivot L, Koussevitzky S, Carretero Y, Bitton F, Moreau L, Causse M (2020) Genetic basis of phenotypic plasticity and genotype x environment interactions in a multi-parental tomato population. J Exp Bot 1:5365-5376. https://doi.org/10.1093/jxp/eraa265

dos Santos Silva PP, e Sousa MB, de Oliveira EJ, Morgante CV, Silva de Oliveira CR, Vieira SL, Borel JC (2021) Genome-wide association study of drought tolerance in cassava. Euphytica 217:60. https://doi.org/10.1007/s10681-021-02800-4

El-Soda M, Malosetti M, Zwaan BJ, Koornneef M, Aarts MGM (2014) Genotype × environment interaction QTL mapping in plants: lessons from Arabidopsis. Trends Plant Sci 19:390-398. https://doi.org/10.1016/j.tplants.2014.01.001

Ferreira ME, Satagopan J, Yandell BS, Williams PH, Osborn TC (1995) Mapping loci controlling vernalization requirement and flowering time in *Brassica napus*. Theor Appl Genet 90:727-732. https://doi.org/10.1007/BF00222140

Finlay KW, Wilkinson GN (1963) The analysis of adaptation in a plant-breeding programme. AustJ Agric Res 14:742-754. https://doi.org/10.1071/AR9630742

Gajardo HA, Wittkop B, Soto-Cerda B, Higgins EE, Parkin IAP, Snowdon RJ, Frederico ML, Iniguez-Luy FL (2015) Association mapping of seed quality traits in *Brassica napus* L. using GWAS and candidate QTL approaches. Mol Breed 35:143. https://doi.org/10.1007/s11032-015-0340-3

Gao X, Starmer J, Martin ER (2008) A multiple testing correction method for genetic association studies using correlated single nucleotide polymorphisms. Genet Epidemiol 32:361-369. https://doi.org/10.1002/gepi.20310

Garin V, Malosetti M, van Eeuwijk F (2020) Multi-parent multi-environment QTL analysis: an illustration with the EU-NAM Flint population. Theor Appl Genet 133:2627-2638. https://doi.org/10.1007/s00122-020-03621-0
Hassani M, Mahmoudi SB, Saremirad A, Taleghani D (2023) Genotype by environment and genotype by yield*trait interactions in sugar beet: analyzing yield stability and determining key traits association. Sci Rep 13:23111. https://doi.org/10.1038/s41598-023-51061-9

Jeuffroy MH, Valentin-Morison M, Champolivier L, Reau R (2006) Azote, rendement et qualité des graines : mise au point et utilisation du modèle Azodyn-colza pour améliorer les performances du colza vis-à-vis de l'azote. OCL 13:388-392. https://doi.org/10.1051/ocl.2006.0090

Josse J, Husson F (2016) missMDA: a package for handling missing values in multivariate data analysis. J Stat Softw. https://doi.org/10.18637/jss.v070.i01

Kuznetsova A, Brockhoff P, Christensen R (2017) lmerTest package: tests in linear mixed effects models. J Stat Softw. https://doi.org/10.18637/jss.v082.i13

Lamichhane JR, Debaek P, Steinberg C, You MP, Barbetti MJ, Aubertot JN (2018) Abiotic and biotic factors affecting crop seed germination and seedling emergence: a conceptual framework. Plant Soil 432:1-28. https://doi.org/10.1007/s11104-018-3780-9

Leveugle M, Duborjal H, West M, Pothier A, Clementi A, Nesi N, Dyrszka E, Pichon J, Riviere N (2015) Rapsodyn Whole Exome capture and genotyping in the polyploid *Brassica napus*, in: JOBIM. Clermont-Ferrand.

Li F, Chen B, Xu K, Gao G, Yan G, Qiao J, Li J, Li H, Li L, Xiao X, Zhang T, Nishio T, Wu X (2016) A genome-wide association study of plant height and primary branch number in rapeseed (*Brassica napus*). Plant Sci 242:169-177. https://doi.org/10.1016/j.plantsci.2015.05.012

Li S, Zhu Y, Varshney RK, Zhan J, Zheng X, Shi J, Wang X, Liu G, Wang H (2020) A systematic dissection of the mechanisms underlying the natural variation of silique number in rapeseed (*Brassica napus* L.) germplasm. Plant Biotech J 18:568-580. https://doi.org/10.1111/pbi.13224

Lippert C, Listgarten J, Liu Y, Kadie CM, Davidson RI, Heckerman D (2011) FaST linear mixed models for genome-wide association studies. Nat Methods 8:833-835. https://doi.org/10.1038/nmeth.1681

Lobell DB, Burke MB, Tebaldi C, Mastrandrea MD, Falcon WP, Naylor RL (2008) Prioritizing climate change adaptation needs for food security in 2030. Science 319:607-610. https://doi.org/10.1126/science.1152339

Lu K, Peng L, Zhang C, Lu J, Yang B, Xiao Z, Liang Y, Xu X, Qu C, Zhang K, Liu L, Zhu Q, Fu M, Yuan X, Li J (2017) Genome-wide association and transcriptome analyses reveal candidate genes underlying yield-determining traits in *Brassica napus*. Front Plant Sci 8:206. https://doi.org/10.3389/fpls.2017.00206

Malosetti M, Ribaut JM, Vargas M, Crossa J, van Eeuwijk FA (2008) A multi-trait multi-environment QTL mixed model with an application to drought and nitrogen stress trials in maize (*Zea mays* L.). Euphytica 161:241-257. https://doi.org/10.1007/s10681-007-9594-0

Malosetti M, Bustos-Korts D, Boer MP, van Eeuwijk FA (2016) Predicting responses in multiple environments: issues in relation to genotype × environment interactions. Crop Sci 56:2210-2222. https://doi.org/10.2135/cropsci2015.05.0311

Mangin B, Casadebaig P, Cadic E, Blanchet N, Boniface M-C, Carrère S, Gouzy J, Legrand L, Mayjonade B, Pouilly N, André T, Coque M, Piquemal J, Laporte M, Vincourt P, Muños S, Langlade NB (2017) Genetic control of plasticity of oil yield for combined abiotic stresses using a joint approach of crop modelling and genome-wide association. Plant, Cell Environ 40:2276-2291. https://doi.org/10.1111/pce.12961

Millet EJ, Welcker C, Kruijer W, Negro S, Coupel-Ledru A, Nicolas SD, Laborde J, Bauland C, Praud S, Ranc N, Presterl T, Tuberosa R, Bedo Z, Draye X, Usadel B, Charcosset A, Van Eeuwijk

F, Tardieu F (2016) Genome-wide analysis of yield in Europe: allelic effects vary with drought and heat scenarios. Plant Physiol 172:749-764. https://doi.org/10.1104/pp.16.00621

Moreau L, Gallais A, Charwsset A (2004) Use of trial clustering to study QTL x environment effects for grain and related traits in maize. Theor Appl Genet 110:92-105. https://doi.org/10.1007/s00122-004-1781-y

Morrison MJ (1993) Heat stress during reproduction in summer rape. Can J Bot 71:303-308. https://doi.org/10.1139/b93-031

Pal L, Sandhu SK, Bhatia D (2021) Genome-wide association study and identification of candidate genes for seed oil content in *Brassica napus*. Euphytica 217:66. https://doi.org/10.1007/s10681-021-02783-2

Parnaudeau V, Jeuffroy MH, Machet JM, Reau R, Bissuel C, Eveillard P, Ostergaard HS, Orson JA (2009) Methods for determining the nitrogen fertiliser requirements of some major arable crops in France. International Fertiliser Society, Cambridge 661:1-26

Quijada PA, Udall JA, Lambert B, Osborn TC (2006) Quantitative trait analysis of seed yield and other complex traits in hybrid spring rapeseed (*Brassica napus* L.): 1. Identification of genomic regions from winter germplasm. Theor Appl Genet 113:549-561. https://doi.org/10.1007/s00122-006-0323-1

Raboanatahiry N, Chao H, Dalin H, Pu S, Yan W, Yu L, Wang B, Li M (2018) QTL alignment for seed yield and yield related traits in *Brassica napus*. Front Plant Sci 9:1127. https://doi.org/10.3389/fpls.2018.01127

Raman H, Shamaya N, Pirathiban R, McVittie B, Raman R, Cullis B, Easton A (2023) Quantitative traits loci for genotype and genotype by environment interaction effects for seed yield plasticity to terminal water-deficit conditions in canola (*Brassica napus* L.) Plants 12:720. https://doi.org/10.3390/plants12040720

Rathke GW, Behrens T, Diepenbrock W (2006) Integrated nitrogen management strategies to improve seed yield, oil content and nitrogen efficiency of winter oilseed rape (*Brassica napus* L.): a review. Agr Ecosyst Environ 117:80-108. https://doi.org/10.1016/j.agee.2006.04.006

Ravier C, Meynard JM, Cohan JP, Jeuffroy MH (2017) Early nitrogen deficiencies favor high yield, grain protein content and N use efficiency in wheat. Eur J Agr 89:16-24. https://doi.org/10.1016/j.eja.2017.06.002

Rémy JC, Hébert J (1977) Le devenir des engrais azotés dans le sol. Compt-Ren L'acad D'agric Franç 63(1977):700-710

Resende RT, Piepho HP, Rosa GJM, Silva-Junior OB, e Silva FF, de Resende MDV, Grattapaglia D (2021) Enviromics in breeding: applications and perspectives on envirotypic-assisted selection. Theor Appl Genet 134:95-112. https://doi.org/10.1007/s00122-020-03684-z

Rincent R, Moreau L, Monod H, Kuhn E, Melchinger A, Malvar RA, Moreno-Gonzalez J, Nicolas S, Madur D, Combes V, Dumas F, Altmann T, Brunel D, Ouzunova M, Flament P, Dubreuil P, Charcosset A, Mary-Huard T (2014) Recovering power in association mapping panels with variable levels of linkage disequilibrium. Genetics 197:375-387. https://doi.org/10.1534/genetics.113.159731

Rousseau-Gueutin M, Belser C, Da Silva C, Richard G, Istace B, Cruaud C, Falentin C, Boideau F, Boutte J, Delourme R, Deniot G, Engelen S, Ferreira de Carvalho J, Lemainque A, Maillet L, Morice J, Wincker P, Denoeud F, Chevre AM, Aury JM (2020) Long-read assembly of the *Brassica napus* reference genome Darmor-bzh. GigaScience 9:1-16. https://doi.org/10.1093/gigascience/giaia137

Shi J, R, Qiu D, Jiang C, Long Y, Morgan C, Bancroft I, Zhao J, Meng J, (2009) Unraveling the complex trait of crop yield with quantitative trait loci mapping in *Brassica napus*. Genetics 182:851-861. https://doi.org/10.1534/genetics.109.101642
Hassani M, Mahmoudi SB, Saremirad A, Taleghani D (2023) Genotype by environment and genotype by yield*trait interactions in sugar beet: analyzing yield stability and determining key traits association. Sci Rep 13:23111. https://doi.org/10.1038/s41598-023-51061-9

Jeuffroy MH, Valentin-Morison M, Champolivier L, Reau R (2006) Azote, rendement et qualité des graines : mise au point et utilisation du modèle Azodyn-colza pour améliorer les performances du colza vis-à-vis de l'azote. OCL 13:388-392. https://doi.org/10.1051/ocl.2006.0090

Josse J, Husson F (2016) missMDA: a package for handling missing values in multivariate data analysis. J Stat Softw. https://doi.org/10.18637/jss.v070.i01

Kuznetsova A, Brockhoff P, Christensen R (2017) lmerTest package: tests in linear mixed effects models. J Stat Softw. https://doi.org/10.18637/jss.v082.i13

Lamichhane JR, Debaek P, Steinberg C, You MP, Barbetti MJ, Aubertot JN (2018) Abiotic and biotic factors affecting crop seed germination and seedling emergence: a conceptual framework. Plant Soil 432:1-28. https://doi.org/10.1007/s11104-018-3780-9

Leveugle M, Duborjal H, West M, Pothier A, Clementi A, Nesi N, Dyrszka E, Pichon J, Riviere N (2015) Rapsodyn Whole Exome capture and genotyping in the polyploid *Brassica napus*, in: JOBIM. Clermont-Ferrand.

Li F, Chen B, Xu K, Gao G, Yan G, Qiao J, Li J, Li H, Li L, Xiao X, Zhang T, Nishio T, Wu X (2016) A genome-wide association study of plant height and primary branch number in rapeseed (*Brassica napus*). Plant Sci 242:169-177. https://doi.org/10.1016/j.plantsci.2015.05.012

Li S, Zhu Y, Varshney RK, Zhan J, Zheng X, Shi J, Wang X, Liu G, Wang H (2020) A systematic dissection of the mechanisms underlying the natural variation of silique number in rapeseed (*Brassica napus* L.) germplasm. Plant Biotech J 18:568-580. https://doi.org/10.1111/pbi.13224

Lippert C, Listgarten J, Liu Y, Kadie CM, Davidson RI, Heckerman D (2011) FaST linear mixed models for genome-wide association studies. Nat Methods 8:833-835. https://doi.org/10.1038/nmeth.1681

Lobell DB, Burke MB, Tebaldi C, Mastrandrea MD, Falcon WP, Naylor RL (2008) Prioritizing climate change adaptation needs for food security in 2030. Science 319:607-610. https://doi.org/10.1126/science.1152339

Lu K, Peng L, Zhang C, Lu J, Yang B, Xiao Z, Liang Y, Xu X, Qu C, Zhang K, Liu L, Zhu Q, Fu M, Yuan X, Li J (2017) Genome-wide association and transcriptome analyses reveal candidate genes underlying yield-determining traits in *Brassica napus*. Front Plant Sci 8:206. https://doi.org/10.3389/fpls.2017.00206

Malosetti M, Ribaut JM, Vargas M, Crossa J, van Eeuwijk FA (2008) A multi-trait multi-environment QTL mixed model with an application to drought and nitrogen stress trials in maize (*Zea mays* L.). Euphytica 161:241-257. https://doi.org/10.1007/s10681-007-9594-0

Malosetti M, Bustos-Korts D, Boer MP, van Eeuwijk FA (2016) Predicting responses in multiple environments: issues in relation to genotype × environment interactions. Crop Sci 56:2210-2222. https://doi.org/10.2135/cropsci2015.05.0311

Mangin B, Casadebaig P, Cadic E, Blanchet N, Boniface M-C, Carrère S, Gouzy J, Legrand L, Mayjonade B, Pouilly N, André T, Coque M, Piquemal J, Laporte M, Vincourt P, Muños S, Langlade NB (2017) Genetic control of plasticity of oil yield for combined abiotic stresses using a joint approach of crop modelling and genome-wide association. Plant, Cell Environ 40:2276-2291. https://doi.org/10.1111/pce.12961

Millet EJ, Welcker C, Kruijer W, Negro S, Coupel-Ledru A, Nicolas SD, Laborde J, Bauland C, Praud S, Ranc N, Presterl T, Tuberosa R, Bedo Z, Draye X, Usadel B, Charcosset A, Van Eeuwijk

F, Tardieu F (2016) Genome-wide analysis of yield in Europe: allelic effects vary with drought and heat scenarios. Plant Physiol 172:749-764. https://doi.org/10.1104/pp.16.00621

Moreau L, Gallais A, Charwsset A (2004) Use of trial clustering to study QTL x environment effects for grain and related traits in maize. Theor Appl Genet 110:92-105. https://doi.org/10.1007/s00122-004-1781-y

Morrison MJ (1993) Heat stress during reproduction in summer rape. Can J Bot 71:303-308. https://doi.org/10.1139/b93-031

Pal L, Sandhu SK, Bhatia D (2021) Genome-wide association study and identification of candidate genes for seed oil content in *Brassica napus*. Euphytica 217:66. https://doi.org/10.1007/s10681-021-02783-2

Parnaudeau V, Jeuffroy MH, Machet JM, Reau R, Bissuel C, Eveillard P, Ostergaard HS, Orson JA (2009) Methods for determining the nitrogen fertiliser requirements of some major arable crops in France. International Fertiliser Society, Cambridge 661:1-26

Quijada PA, Udall JA, Lambert B, Osborn TC (2006) Quantitative trait analysis of seed yield and other complex traits in hybrid spring rapeseed (*Brassica napus* L.): 1. Identification of genomic regions from winter germplasm. Theor Appl Genet 113:549-561. https://doi.org/10.1007/s00122-006-0323-1

Raboanatahiry N, Chao H, Dalin H, Pu S, Yan W, Yu L, Wang B, Li M (2018) QTL alignment for seed yield and yield related traits in *Brassica napus*. Front Plant Sci 9:1127. https://doi.org/10.3389/fpls.2018.01127

Raman H, Shamaya N, Pirathiban R, McVittie B, Raman R, Cullis B, Easton A (2023) Quantitative traits loci for genotype and genotype by environment interaction effects for seed yield plasticity to terminal water-deficit conditions in canola (*Brassica napus* L.) Plants 12:720. https://doi.org/10.3390/plants12040720

Rathke GW, Behrens T, Diepenbrock W (2006) Integrated nitrogen management strategies to improve seed yield, oil content and nitrogen efficiency of winter oilseed rape (*Brassica napus* L.): a review. Agr Ecosyst Environ 117:80-108. https://doi.org/10.1016/j.agee.2006.04.006

Ravier C, Meynard JM, Cohan JP, Jeuffroy MH (2017) Early nitrogen deficiencies favor high yield, grain protein content and N use efficiency in wheat. Eur J Agr 89:16-24. https://doi.org/10.1016/j.eja.2017.06.002

Rémy JC, Hébert J (1977) Le devenir des engrais azotés dans le sol. Compt-Ren L'acad D'agric Franç 63(1977):700-710

Resende RT, Piepho HP, Rosa GJM, Silva-Junior OB, e Silva FF, de Resende MDV, Grattapaglia D (2021) Enviromics in breeding: applications and perspectives on envirotypic-assisted selection. Theor Appl Genet 134:95-112. https://doi.org/10.1007/s00122-020-03684-z

Rincent R, Moreau L, Monod H, Kuhn E, Melchinger A, Malvar RA, Moreno-Gonzalez J, Nicolas S, Madur D, Combes V, Dumas F, Altmann T, Brunel D, Ouzunova M, Flament P, Dubreuil P, Charcosset A, Mary-Huard T (2014) Recovering power in association mapping panels with variable levels of linkage disequilibrium. Genetics 197:375-387. https://doi.org/10.1534/genetics.113.159731

Rousseau-Gueutin M, Belser C, Da Silva C, Richard G, Istace B, Cruaud C, Falentin C, Boideau F, Boutte J, Delourme R, Deniot G, Engelen S, Ferreira de Carvalho J, Lemainque A, Maillet L, Morice J, Wincker P, Denoeud F, Chevre AM, Aury JM (2020) Long-read assembly of the *Brassica napus* reference genome Darmor-bzh. GigaScience 9:1-16. https://doi.org/10.1093/gigascience/giaia137

Shi J, R, Qiu D, Jiang C, Long Y, Morgan C, Bancroft I, Zhao J, Meng J, (2009) Unraveling the complex trait of crop yield with quantitative trait loci mapping in *Brassica napus*. Genetics 182:851-861. https://doi.org/10.1534/genetics.109.101642
Smith A, Norman A, Kuchel H, Cullis B (2021) Plant variety selection using interaction classes derived from factor analytic linear mixed models: models with independent variety effects. *Front Plant Sci* 12:737462. https://doi.org/10.3389/fpls.2021.737462

Snowdon RJ, Wittkop B, Chen TW, Stahl A (2021) Crop adaptation to climate change as a consequence of long-term breeding. *Theor Appl Genet* 134:1613–1623. https://doi.org/10.1007/s00122-020-03729-3

Sun F, Liu J, Hua W, Sun X, Wang X, Wang H (2016) Identification of stable QTLs for seed oil content by combined linkage and association mapping in *Brassica napus*. *Plant Sci* 252:388–399. https://doi.org/10.1016/j.plantsci.2016.09.001

Tang Y, Zhang G, Jiang X, Shen S, Guan M, Tang Y, Sun F, Hu R, Chen S, Zhao H, Li J, Lu K, Yin N, Qu C (2023) Genome-wide association study of glucosinolate metabolites (mGWAS) in *Brassica napus* L. *Plants* 12(3):639. https://doi.org/10.3390/plants12030639

Touzy G, Rincent R, Bogard M, Lafarge S, Dubreuil P, Mini A, Deswarte JC, Beauchêne K, Le Gouis J, Praud S (2019) Using environmental clustering to identify specific drought tolerance QTLs in bread wheat (*T. aestivum* L.) *Theor Appl Genet* 132:2859–2880. https://doi.org/10.1007/s00122-019-03393-2

Udall JA, Quijada PA, Lambert B, Osborn TC (2006) Quantitative trait analysis of seed yield and other complex traits in hybrid spring rapeseed (*Brassica napus* L.): 2. Identification of alleles from unadapted germplasm. *Theor Appl Genet* 113:597–609. https://doi.org/10.1007/s00122-006-0324-0

van Eeuwijk FA, Bink MCAM, Chenu K, Chapman SC (2010) Detection and use of QTL for complex traits in multiple environments. *Curr Opin Plant Biol* 13:193–205. https://doi.org/10.1016/j.pbi.2010.01.001

van Eeuwijk FA, Bustos-Korts DV, Malosetti M (2016) What should students in plant breeding know about the statistical aspects of genotype x environment interactions? *Crop Sci* 56:2119–2140. https://doi.org/10.2135/cropsci2015.06.0375

Wang B, Wu Z, Li Z, Zhang Q, Hu J, Xiao Y, Cai D, Wu J, King GJ, Li H, Liu K (2018) Dissection of the genetic architecture of three seed-quality traits and consequences for breeding in *Brassica napus*. *Plant Biotechnol J* 16:1336–1348. https://doi.org/10.1111/pbi.12873

Wang TC, Casadebaig P, Chen TW (2023) More than 1000 genotypes are required to derive robust relationships between yield, yield stability and physiological parameters: a computational study on wheat crop. *Theor Appl Genet* 136:34. https://doi.org/10.1007/s00122-023-04264-7

Weir BS, Cockerham CC (1984) Estimating F-statistics for the analysis of population structure. *Evolution* 38:1358–1370. https://doi.org/10.1111/j.1558-5646.1984.tb05657.x

Wezel A, Casagrande M, Celette F, Vian JF, Ferrer A, Peigné J (2014) Agroecological practices for sustainable agriculture. A review. *Agron Sustain Dev* 34:1–20. https://doi.org/10.1007/s13593-013-0180-7

Wricke G (1962) Evaluation method for recording ecological differences in field trials. *Z Pflanzenzücht* 47:92–96

Xavier A, Jarquin D, Howard R, Ramasubramanian V, Specht JE, Graef GL, Beavis WD, Diers BW, Song Q, Cregan PB, Neslon R, Mian R, Shannon JG, McHale L, Wang D, Schapaugh W, Lorenz AJ, Xu S, Muir WM, Rainey KM (2018) Genome-wide analysis of grain yield stability and environmental interactions in a multiparental soybean population. *G3 Genes Genom Genet* 8:519–529. https://doi.org/10.1534/g3.117.300300

Xie Y, Xu J, Tian G, Xie L, Xu B, Liu K, Zhang X (2020) Unraveling yield-related traits with QTL analysis and dissection of QTL x environment interaction using a high-density bin map in rapeseed (*Brassica napus* L.). *Euphytica* 216:171. https://doi.org/10.1007/s10681-020-02708-5

Young LW, Wilen RW, Bonham-smith PC (2004) High temperature stress of *Brassica napus* during flowering reduces micro- and megagametophyte fertility, induces fruit abortion, and disrupts seed production. *J Exp Bot* 55:485–495. https://doi.org/10.1093/jxb/erh038

Zandberg JD, Frenadez CT, Danilevicz MF, Thomas WJW, Edwards D, Batley J (2022) The global assessment of oilseed Brassica crop species yield, yield stability and the underlying genetics. *Plants* 11:2740. https://doi.org/10.3390/plants11202740

Zhao W, Wang X, Wang H, Tian J, Li B, Chen L, Chao H, Long Y, Xiang J, Gan J, Liang W, Li M (2016) Genome-Wide Identification of QTL for seed yield and yield-related traits and construction of a high-density consensus map for QTL comparison in *Brassica napus*. *Front Plant Sci* 7:17. https://doi.org/10.3389/fpls.2016.00017

Zheng M, Peng C, Liu H, Tang M, Yang H, Li X, Liu J, Sun X, Wang X, Xu J, Hua W, Wang H (2017) Genome-wide association study reveals candidate genes for control of plant height, branch initiation height and branch number in rapeseed (*Brassica napus* L.). *Front Plant Sci* 8:1246. https://doi.org/10.3389/fpls.2017.01246

Zou M, Shi T, Wang W, Ding G, Xu F, Shi L (2022) Genetic dissection of seed yield and yield-related traits in *Brassica napus* grown with constrasting nitrogen supplies. *Mol Breed* 42:15. https://doi.org/10.1007/s11032-022-01281-0

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Smith A, Norman A, Kuchel H, Cullis B (2021) Plant variety selection using interaction classes derived from factor analytic linear mixed models: models with independent variety effects. *Front Plant Sci* 12:737462. https://doi.org/10.3389/fpls.2021.737462

Snowdon RJ, Wittkop B, Chen TW, Stahl A (2021) Crop adaptation to climate change as a consequence of long-term breeding. *Theor Appl Genet* 134:1613–1623. https://doi.org/10.1007/s00122-020-03729-3

Sun F, Liu J, Hua W, Sun X, Wang X, Wang H (2016) Identification of stable QTLs for seed oil content by combined linkage and association mapping in *Brassica napus*. *Plant Sci* 252:388–399. https://doi.org/10.1016/j.plantsci.2016.09.001

Tang Y, Zhang G, Jiang X, Shen S, Guan M, Tang Y, Sun F, Hu R, Chen S, Zhao H, Li J, Lu K, Yin N, Qu C (2023) Genome-wide association study of glucosinolate metabolites (mGWAS) in *Brassica napus* L. *Plants* 12(3):639. https://doi.org/10.3390/plants12030639

Touzy G, Rincent R, Bogard M, Lafarge S, Dubreuil P, Mini A, Deswarte JC, Beauchêne K, Le Gouis J, Praud S (2019) Using environmental clustering to identify specific drought tolerance QTLs in bread wheat (*T. aestivum* L.) *Theor Appl Genet* 132:2859–2880. https://doi.org/10.1007/s00122-019-03393-2

Udall JA, Quijada PA, Lambert B, Osborn TC (2006) Quantitative trait analysis of seed yield and other complex traits in hybrid spring rapeseed (*Brassica napus* L.): 2. Identification of alleles from unadapted germplasm. *Theor Appl Genet* 113:597–609. https://doi.org/10.1007/s00122-006-0324-0

van Eeuwijk FA, Bink MCAM, Chenu K, Chapman SC (2010) Detection and use of QTL for complex traits in multiple environments. *Curr Opin Plant Biol* 13:193–205. https://doi.org/10.1016/j.pbi.2010.01.001

van Eeuwijk FA, Bustos-Korts DV, Malosetti M (2016) What should students in plant breeding know about the statistical aspects of genotype x environment interactions? *Crop Sci* 56:2119–2140. https://doi.org/10.2135/cropsci2015.06.0375

Wang B, Wu Z, Li Z, Zhang Q, Hu J, Xiao Y, Cai D, Wu J, King GJ, Li H, Liu K (2018) Dissection of the genetic architecture of three seed-quality traits and consequences for breeding in *Brassica napus*. *Plant Biotechnol J* 16:1336–1348. https://doi.org/10.1111/pbi.12873

Wang TC, Casadebaig P, Chen TW (2023) More than 1000 genotypes are required to derive robust relationships between yield, yield stability and physiological parameters: a computational study on wheat crop. *Theor Appl Genet* 136:34. https://doi.org/10.1007/s00122-023-04264-7

Weir BS, Cockerham CC (1984) Estimating F-statistics for the analysis of population structure. *Evolution* 38:1358–1370. https://doi.org/10.1111/j.1558-5646.1984.tb05657.x

Wezel A, Casagrande M, Celette F, Vian JF, Ferrer A, Peigné J (2014) Agroecological practices for sustainable agriculture. A review. *Agron Sustain Dev* 34:1–20. https://doi.org/10.1007/s13593-013-0180-7

Wricke G (1962) Evaluation method for recording ecological differences in field trials. *Z Pflanzenzücht* 47:92–96

Xavier A, Jarquin D, Howard R, Ramasubramanian V, Specht JE, Graef GL, Beavis WD, Diers BW, Song Q, Cregan PB, Neslon R, Mian R, Shannon JG, McHale L, Wang D, Schapaugh W, Lorenz AJ, Xu S, Muir WM, Rainey KM (2018) Genome-wide analysis of grain yield stability and environmental interactions in a multiparental soybean population. *G3 Genes Genom Genet* 8:519–529. https://doi.org/10.1534/g3.117.300300

Xie Y, Xu J, Tian G, Xie L, Xu B, Liu K, Zhang X (2020) Unraveling yield-related traits with QTL analysis and dissection of QTL x environment interaction using a high-density bin map in rapeseed (*Brassica napus* L.). *Euphytica* 216:171. https://doi.org/10.1007/s10681-020-02708-5

Young LW, Wilen RW, Bonham-smith PC (2004) High temperature stress of *Brassica napus* during flowering reduces micro- and megagametophyte fertility, induces fruit abortion, and disrupts seed production. *J Exp Bot* 55:485–495. https://doi.org/10.1093/jxb/erh038

Zandberg JD, Frenadez CT, Danilevicz MF, Thomas WJW, Edwards D, Batley J (2022) The global assessment of oilseed Brassica crop species yield, yield stability and the underlying genetics. *Plants* 11:2740. https://doi.org/10.3390/plants11202740

Zhao W, Wang X, Wang H, Tian J, Li B, Chen L, Chao H, Long Y, Xiang J, Gan J, Liang W, Li M (2016) Genome-Wide Identification of QTL for seed yield and yield-related traits and construction of a high-density consensus map for QTL comparison in *Brassica napus*. *Front Plant Sci* 7:17. https://doi.org/10.3389/fpls.2016.00017

Zheng M, Peng C, Liu H, Tang M, Yang H, Li X, Liu J, Sun X, Wang X, Xu J, Hua W, Wang H (2017) Genome-wide association study reveals candidate genes for control of plant height, branch initiation height and branch number in rapeseed (*Brassica napus* L.). *Front Plant Sci* 8:1246. https://doi.org/10.3389/fpls.2017.01246

Zou M, Shi T, Wang W, Ding G, Xu F, Shi L (2022) Genetic dissection of seed yield and yield-related traits in *Brassica napus* grown with constrasting nitrogen supplies. *Mol Breed* 42:15. https://doi.org/10.1007/s11032-022-01281-0

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.