Optimal Contribution Selection Improves the Rate of Genetic Gain in Grain Yield and Yield Stability in Spring Canola in Australia and Canada

Crop breeding must achieve higher rates of genetic gain in grain yield (GY) and yield stability to meet future food demands in a changing climate. Optimal contributions selection (OCS) based on an index of key economic traits should increase the rate of genetic gain while minimising population inbreeding. Here we apply OCS in a global spring oilseed rape (canola) breeding program during three cycles of S0,1 family selection in 2016, 2018, and 2020, with several field trials per cycle in Australia and Canada. Economic weights in the index promoted high GY, seed oil, protein in meal, and Phoma stem canker (blackleg) disease resistance while maintaining plant height, flowering time, oleic acid, and seed size and decreasing glucosinolate content. After factor analytic modelling of the genotype-by-environment interaction for the additive effects, the linear rate of genetic gain in GY across cycles was 0.059 or 0.087 t ha−1 y−1 (2.9% or 4.3% y−1) based on genotype scores for the first factor (f1) expressed in trait units or average predicted breeding values across environments, respectively. Both GY and yield stability, defined as the root-mean-square deviation from the regression line associated with f1, were predicted to improve in the next cycle with a low achieved mean parental coancestry (0.087). These methods achieved rapid genetic gain in GY and other traits and are predicted to improve yield stability across global spring canola environments.


Introduction
The rate of gain in grain yield (GY) in the world's major grain crops ranges from 0.9 to 1.6% y −1 and should at least double to meet the expected global demand for grain crops in 2050 [1][2][3]. In some regions, GY is stagnating or declining due to climate change [2,4]. Plant breeders face a major challenge to accelerate genetic gain in GY during the forthcoming period of global population growth and climate change [1,5], and new methods of rapid breeding should involve wide-scale phenotyping, accurate selection, and international exchange of elite varieties [6].
One potential method of accelerating genetic gain in self-pollinating crops is the use of rapid cycles of early-generation recurrent selection [7]. Such methods were attempted previously and were successful for one or a few traits [8] but fell out of favour due to their lack of application to commercial crop breeding, where multiple traits were under selection. The aim of this study was to field-test stochastic models of early-generation recurrent selection in a self-pollinating crop [9,10], which showed great promise for high genetic gain in multiple economic traits. This field validation of genetic gain occurred within a commercial spring oilseed rape (canola) breeding program.
Optimal contributions selection (OCS) was first used in self-pollinating crops in stochastic models of the additive effects in rapid cycles of early generation selection [9,10,25]. OCS in MateSel [26] includes a crossing plan which optimises parental contributions and balances genetic gain for the index under a number of possible constraints or weightings on factors that affect genetic diversity, inbreeding rate and inbreeding level [27,28]. Importantly, MateSel permits the breeder to predict the genetic gain for each trait and achieved parental coancestry in the next cycle. This allows for timely adjustments in the breeding program for changes in future desired gains, such as adjustments to weightings on traits in the economic index or changes to selection targets.
In S0,1 family selection, the S0-derived S1 family provides an estimate of the breeding value of the S0 individual for selection purposes [20]. In stochastic models, rapid cycles of S0,1 family selection (two years per cycle) with OCS and optimised mating designs improved the rate of long-term genetic gain and reduced the rate of population inbreeding compared to truncation selection and random mating among selected parents [10]. For practical reasons, the S0,1 family seeds are grown in plots in field trials, and the agronomic and harvested grain traits recorded on these plots are used to predict the breeding value of the S0 individuals. Accurate PBV (r > 0.80) were generated on S0 individuals for GY, disease resistance, agronomic traits, and seed quality in stochastic models of S0,1 family selection based on a deep and highly interconnected pedigree [9,10].
In this field study, depending on seed availability, each S0,1 family was grown in replicated plots at multiple sites and in multiple regions; we grew concurrent field trials of the same S0,1 families in Australia and Canada ( Figure 1). We included data for GY and several major economic traits including flowering time, plant height, seed quality traits, and resistance to Phoma stem canker or blackleg disease caused by Leptosphaeria maculans. Resistance to blackleg disease is important for secure production of oilseed rape in Australia, Canada, and Europe [29,30].  Table 1.  Table 1.
Augmented S 0,1 family selection exploits complex relationships between genotypes within and between cycles [10] to generate accurate PBV in an individual model analysis with phenotypic and relationship data from all cycles of selection [17]. As a side bonus, the method also results in near-homozygous lines after two or three cycles that are ready for commercial evaluation [10]. Genomic relationship information can be added at any time and combined with pedigree relationship information in single-step genomic prediction [31]. New germplasm can be added to the program in any cycle. In some crops, it may be necessary to bulk seeds from single S 0 plants over 2 selfing generations (S 0,2 bulks) to generate sufficient seed for phenotyping [20], such as in common bean [32]. Table 1. (a) Cycle number, trial code, trial location and country of each trial site, the number of plots, ranges, rows, genotypes (including number tested in a single replicate) in each trial, and (b) connectivity of genotypes across trials in the study. Genotypes defined as S 0 were tested in field trials as S 0,1 bulks; likewise F 2 genotypes were tested as F 2,3 bulks, S 2 genotypes were tested as S 2,3 bulks, and so on. State abbreviations (Australia): Western Australia (WA), Victoria (VIC). Province abbreviation (Canada): Manitoba (MB).  In this study, we use MMM-FA analysis of multi-environment trials across cycles to assess the rate of genetic gain for GY and several economic traits inside an active global spring oilseed rape (canola) breeding program during three cycles of augmented S 0,1 family selection. We compare the rate of genetic gain in GY across cycles by two methods, OP and average PBV across environments. For the first time in crop breeding, we include yield stability across global environments based on RMSD in the economic index. This study also represents the first major attempt to integrate spring oilseed rape germplasm across the southern and northern hemispheres.

Terminology
We use the traditional definitions in crop breeding of 'F 1 ' as progeny of crosses between near-homozygous founder varieties, and 'S 0 ' as progeny of crosses between heterozygous parent plants. F 1 progeny are heterozygous and genetically uniform, and segregation occurs in the F 2 generation after self-pollination. S 0 progeny are non-inbred and therefore heterozygous and heterogeneous, that is, each S 0 progeny is a unique genotype and segregation occurs in the S 0 for multiple traits. Each plant in the ancestral pedigree is a unique genotype and has a relationship with every other genotype in the pedigree, whether the genotype is derived from crossing or selfing ( Figure 2). S 0,1 represents the selfed progeny from an S 0 plant [20]. may occur within a cycle. In this study, S0,1 family selection [20] is augmented by taking forwards S2 self progeny of parent plants in addition to S0 cross progeny for phenotyping in the next cycle (in field plots of S2,3 or S0,1 families), and this increases the number of collateral relatives and pedigree linkages across cycles of selection ( Figure 2).

Figure 2.
Two cycles of augmented S0,1 family selection showing genotypes in the pedigree. In cycle 1, inbred founder parents from Australia (AU) were crossed with parents from Europe (EU) or Canada (CA), and their F1 progeny were intercrossed to generate S0 and F2 progeny (circles). Individual S0 and F2 progeny were selfed and the self-seed was sown in S0,1 and F2,3 plots in cycle 1 field trials (rectangles). Phenotypic data from S0,1 and F2,3 plots in cycle 1 field trials were used to calculate the predicted breeding value of S0 and F2 progeny. Crossing in cycle 2 was based on remnant S0,1 or F2,3 seed of superior S0 and F2 progeny. A similar process occurred in cycle 2, with phenotypic data from S0,1 and S2,3 plots used to assess the predicted breeding value of S0 and S2 progeny genotypes, respectively. From cycle 2 onwards, selection was based on an economic index composed of all traits, and mating designs were derived from optimal contributions selection. Two cycles of augmented S 0,1 family selection showing genotypes in the pedigree. In cycle 1, inbred founder parents from Australia (AU) were crossed with parents from Europe (EU) or Canada (CA), and their F 1 progeny were intercrossed to generate S 0 and F 2 progeny (circles). Individual S 0 and F 2 progeny were selfed and the self-seed was sown in S 0,1 and F 2,3 plots in cycle 1 field trials (rectangles). Phenotypic data from S 0,1 and F 2,3 plots in cycle 1 field trials were used to calculate the predicted breeding value of S 0 and F 2 progeny. Crossing in cycle 2 was based on remnant S 0,1 or F 2,3 seed of superior S 0 and F 2 progeny. A similar process occurred in cycle 2, with phenotypic data from S 0,1 and S 2,3 plots used to assess the predicted breeding value of S 0 and S 2 progeny genotypes, respectively. From cycle 2 onwards, selection was based on an economic index composed of all traits, and mating designs were derived from optimal contributions selection.
In this study, the term 'cycle' refers to a cycle of recurrent selection, that is, the generation interval (L) as used in the breeder's equation [20]. One cycle in this study takes two years, which includes the time taken to cross the parent plants, self the S 0 to generate S 0,1 families, phenotype the S 0,1 families in the field trials and the laboratory, and select parents for crossing to begin the new cycle ( Figure 2). One or more selfing generations may occur within a cycle. In this study, S 0,1 family selection [20] is augmented by taking forwards S 2 self progeny of parent plants in addition to S 0 cross progeny for phenotyping in the next cycle (in field plots of S 2,3 or S 0,1 families), and this increases the number of collateral relatives and pedigree linkages across cycles of selection ( Figure 2).

Founder Population, Crossing and Selfing to Begin Cycle 1
The founder population included 32 southern hemisphere (SH) elite spring canola lines from Australia (AU1 to AU32) and 32 northern hemisphere (NH) elite spring canola lines from Canada (CA1 to CA16) and Europe (EU1 to EU16) which were intercrossed in 32 pairwise combinations (SH × NH) in 2012. Crossing occurred in glasshouses at The University of Western Australia Field Station, Shenton Park, Western Australia. Their F 1 progeny were intercrossed from January to June 2013 in 16 combinations (F 1 × F 1 ). For example, the F 1 of AU2 × EU2 was intercrossed with the F 1 of CA5 × AU5 to generate S 0 progeny D1 ( Figure 2, Online Resource S1). All female parents were selected for triazine tolerance and therefore the breeding population was uniformly tolerant of triazine herbicides due to the cytoplasmic control of this trait [33]. This avoided any potential bias in the project stemming from segregation among progeny for herbicide tolerance.
In cycle 1, S 0 and F 2 progeny were grown in the field in Chile from October 2013 to February 2014 where selfing occurred inside pollination bags on single plants. S 0,1 and F 2,3 families were returned to Australia where they were sown in the cycle 1 field trial in May 2014 ( Figure 2). Recurrent selection cycles continued every two years and 70 new migrants (35 from SH and 35 from NH) were added during cycles 2, 3 and 4 as F 1 progeny of the crosses SH × NH (Online Resource S2).

Cycles of Augmented S 0,1 Family Selection
The cycle 1 field trial (trial code 2014AU1) was grown in Western Australia from May to November 2014 and included 668 S 0,1 families derived from F 1 × F 1 matings augmented with 577 F 2,3 families derived from F 2 seed harvested from each F 1 parent plant (Table 1). In addition, several S 1,2 and F 3,4 families were evaluated in 2014 (Table 1). Remnant seeds of S 0,1 , F 2,3 , S 1,2, and F 3,4 families were stored for potential future use as parents to begin cycle 2. Crossing to begin cycle 2 was among selections from cycle 1 progeny which performed well for GY and grain quality traits in trial 2014AU1, and 170 matings were made by the breeder to combine different pedigrees in male and female parents. Crossing decisions to begin cycles 3, 4, and 5 were based on optimal contributions selection (OCS) as described below in Section 2.8 'Optimal Contributions Selection and Crossing'.
The process of augmented S 0,1 family selection continued in two-year cycles with phenotyping of S 0,1 and S 2,3 families in multiple field trials in Australia and Canada in cycles 2 (2016), 3 (2018) and 4 (2020) ( Table 1, Online Resource S2). There were two field trials in cycle 2, three in cycle 3, and three in cycle 4. Trials in Canada were located at Sun Valley, Manitoba (trial codes 2018CA1 and 2020CA1), and in Australia were located in canola cropping zones near Perth, Western Australia (trial codes 2014AU1, 2016AU1, 2018AU1, and 2020AU1) and in northern Victoria (trial codes 2016AU2, 2018AU2 and 2020AU2). Sites identified as AU2 were 'disease nurseries' sown on straw of the previous year's canola crop to promote high levels of Phoma stem canker (blackleg) disease (Table 1a). Connectivity of genotypes, that is, the occurrence of the same genotypes in different trials, was high across trials within cycles but low across trials between cycles (Table 1b).

Pedigree Structure and Relationships
As required for pedigree analysis in ASReml-R v4 [34], each genotype was described in the pedigree file with the name of the genotype, the names of its male and female parents and its level of selfing ('fgen'). This was preceded in the file by rows which showed the same information for each parent. The pedigree included ancestors of founder varieties and control varieties. For most genotypes, fgen was the default value of zero because selfing was explicitly described in the pedigree file (the male and female parents were identical). For most inbred ancestral varieties fgen was set at 5 and for those derived from doubled haploids fgen was set at 10. The pedigree file was converted into an additive genetic relationship matrix and an inverse relationship matrix was formed for analysis in ASReml-R v4 [34]. The pedigree was highly interconnected within and across four cycles as a result of selfing and crossing in augmented S 0,1 family selection ( Figure 2, Online Resource S1).
Inbreeding coefficients (F) (Online Resource S3) and coefficients of coancestry ( f ) (Online Resource S4) were calculated from the additive genetic relationship matrix (Amatrix). F was calculated as the A-value on the diagonal of the A-matrix minus 1, and f was calculated as 1 2 the A-value of the 2 parent genotypes.

Field Trials and Phenotyping
S 0 and S 2 genotypes were evaluated as S 0,1 or S 2,3 families in field plots ( Figure 2) with one or two replicates per site and at multiple trial sites per cycle depending on seed availability. Phenotypes of the S 0,1 or S 2,3 families were used to predict the breeding value of the S 0 and S 2 individuals for selection purposes [20], as described below in Section 2.  [35].
Partially replicated field trials [36] were designed in DiGGer (available from http: //nswdpibiom.org/austatgen/software/, accessed 9 January 2023) with the default spatial model at sites in Australia and Canada (Table 1). Each site was grown under natural rainfall with date of sowing and agronomic treatments similar to commercial canola production in the region. Each field plot was 6-m long and six rows wide (1.8 m centre-to-centre) and sown with 3 g seed per plot. Plots were cut back to 4 m prior to harvest. Control cultivars and S 0,1 and S 2,3 families with sufficient seed were replicated at each site ( Table 1). The proportion of plots at a site composed of single-replicate genotypes varied from 37.7% to 83.9% over the four cycles (Table 1).
During the growing season at some sites, plots were scored for days to 50% flowering (DTF) and plant height at maturity (PlHt, cm). At disease nursery sites (AU2), plots were also scored for blackleg resistance at maturity. The blackleg (BL) disease resistance score was based on the number of plants that were visibly dead or lodging as a result of blackleg disease in the plot at pod filling stage, and BL scores ranged from 1 (very susceptible, no plants survived blackleg disease) to 9 (very resistant, all plants survived) ( Figure 3). Two assessors scored each plot and the average of their scores was recorded as the BL score on the plot.
Field trial plots were harvested by small plot harvester and the weight of harvested seed per plot was converted to grain yield per hectare (GY, t ha −1 ). Harvested grain from each plot was stored in a dry environment and one hundred seeds were randomly extracted from harvest bags to measure 100 seed weight (SW100, g).
Samples of harvested seed from each plot (5 g) were assessed by near-infrared radiation spectroscopy for moisture content (%), seed oil (Oil, % of seed, adjusted to 6% moisture), protein in meal (ProM, % of meal, adjusted to 10% moisture), glucosinolates (GSL, µmole g −1 ), and oleic acid (OL, % of total fatty acids) after pre-calibration against standard seed samples with known levels of moisture, Oil, ProM, GLS, and OL [37].

Preliminary Single Site Analysis
Univariate single-site analyses were conducted in ASReml-R v4 [34] which produces residual maximum likelihood (REML) estimates of the variance parameters [43] and best linear unbiased predictions (BLUP) [44] of the random effects. Univariate analyses were general linear mixed models which followed the structural specifications presented in Gilmour et al. [45] (Online Resource S5).
The spatial variation of each trait within each trial was assessed following the mixed model approach described by Gilmour et al. [46] and Stefanova et al. [47]. Trends along rows and columns were assessed in an autoregressive residual model (AR1 × AR1) and significant linear row and column trends and/or random row and column effects were fitted to the model following Stefanova and Buirchell [21] and Beeck et al. [38]. The effects of local spatial trends were assessed using a plot of residuals, the sample variogram, row and column faces of the empirical variogram, and REML likelihood ratio tests (LRT) [47].
'Genotype' was considered a random genetic effect since large numbers of progeny genotypes were randomly assigned to field trials, there was no prior selection by the breeder, and the goal was to rank progeny genotypes for relative performance. Random effects were assessed for significance by the Z-test of variance components and fixed effects were assessed by the significance of the Wald statistic.
Finally, the additive genetic relationship matrix was added to the model to estimate the additive and nonadditive genetic effects (Online Resource S5). The subsequent improvement in the model was assessed by REML LRT, Akaike information criterion (AIC) [48], and Bayesian information criterion (BIC) [49]. Possible outliers were assessed and removed.

Base Univariate Model Across-Sites
The base model univariate analysis across-sites was a MMM and included the additive genetic relationship matrix and significant fixed and random terms from the single site analysis but ignored covariances between sites, following Stefanova and Buirchell [21] and Beeck et al. [38]. 'Site' was considered a fixed effect. The main effect of 'Genotype' was excluded and 'Genotype × Site' was included as a random effect since the goal was to rank genotypes for relative performance at each site (Online Resource S5). Significant spatial terms from the single-site analyses were transferred to the base model. Variance components were estimated for additive and nonadditive genetic effects at each site and narrow-sense heritability for each trait at each site was calculated as the ratio of the additive variance component divided by the sum of additive, nonadditive, and error variance components. The average accuracy of the additive effects (r i ) was calculated for the i-th individual for each trait at each site, following Gilmour et al. [45] for the animal model.

Factor Analytic Modelling of Each Trait across Sites
Across-sites analysis of each trait proceeded with an MMM-FA model (Online Resource S5) following the approach in Stefanova and Buirchell [21] and Beeck et al. [38] where 'Site' was considered a fixed effect and 'Genotype × Site' a random effect. MMM-FA included the additive genetic relationship matrix and significant fixed and random spatial terms from the base model. Both additive and non-additive genetic effects were assessed in the MMM-FA when significant [38]. Analysis proceeded sequentially in ASReml-R from the base model to first-order factor FA(1), second-order factor FA(2), and so on until the optimum MMM-FA model was chosen based on a non-significant or trivial improvement in the REML LRT, AIC and BIC tests or percentage variance accounted for in the next higher order model.
The PBV for an individual was calculated by averaging the predicted variety means across environments based on the optimum MMM-FA model [21].
Genetic correlations of additive effects of genotypes across sites were obtained from the optimum MMM-FA model.

Genotype Overall Performance and Stability
The average PBV for an individual across environments was a measure of genotype performance across environments based on the optimum MMM-FA model, following the approach for variety predicted means across environments [21]. A second measure of genotype performance across environments was overall performance (OP) which was based on f 1 scores in the optimum MMM-FA model and expressed in trait units. OP was deemed to be a generalised main effect of variety performance when scale differed across environments [22].
We compared two measures of stability of genotype performance across environments: scores for the f 2 from the optimum MMM-FA model [21] and the root-mean-square deviation (RMSD) from the regression line associated with f 1 in the optimum MMM-FA model [22].

Economic Index
The PBV for each trait from the optimum MMM-FA model for each genotype was weighted by an economic value and summed across traits to form an economic index for each genotype in the ancestral pedigree. An example of a hypothetical genotype is provided in Table 2.
Firstly, the total GY of each genotype was calculated as the sum of the PBV GY of the genotype and population mean GY, expressed in t ha −1 . The total GY of the genotype was multiplied by the average contemporary market value of harvested grain in US$ t −1 to find the contribution of the GY of the genotype to the economic index in US$ ha −1 ( Table 2).
The economic weight for Oil was based on the bonification rate for oil in Australia where there is a 1.5% grain price premium (or deduction) for each 1% seed oil above (or below) the base seed oil content of 42% [50]. The economic weight for +1 unit PBV Oil is converted into US$ t −1 (US$8.25 t −1 ) based on grain price US$550 t −1 . The contribution of Oil to the economic index in US$ ha −1 was calculated as the PBV Oil of the genotype × economic weight for +1% Oil (US$8.25 t −1 ) × total GY of the genotype (2.343 t ha −1 ) ( Table 2). Table 2. Calculation of economic index based on economic weights for traits to achieve desired genetic gains in a hypothetical genotype in cycle 4. The assumed grain price is US$550 t −1 . The economic weight of a trait for +1 unit predicted breeding value (PBV) is shown as % grain price and in US$ t −1 . The contribution of a trait to the economic index in US$ ha −1 is calculated by multiplying the economic weight for +1 PBV unit (US$ t −1 ) × PBV of the individual in units of trait × total GY of individual (t ha −1 ). The economic weights on ProM, DTF, PlHt, GSL, OL, SW100, and BL in the economic index were subjectively assigned and modified to achieve the desired outcomes in the next cycle of the optimised mating design. Positive weights were assigned for traits where genetic gains were desired upwards (ProM, BL, SW100), and negative weights where genetic gains downwards were desired (DTF, PlHt, GSL), or zero weights where no change was desired (OL) ( Table 2).

Optimal Contributions Selection and Crossing
The inter-crossing of founders to begin cycle 1 is described above in Section 2.2 'Founder Population, Crossing and Selfing to Begin Cycle 1'. Crossing to begin cycle 2 was among selections from cycle 1 progeny which performed well for GY and grain quality traits in trial 2014AU1, and 170 matings were made by the breeder with an emphasis to combine different pedigrees in male and female parents (Online Resource S2).
Crossing designs for cycles 3, 4, and 5 were based on OCS in the implementation platform 'MateSel' for the construction of an optimised mating design [26]. The selection and mate allocation method used previously [9,10,28] involves a function whose key components relate to genetic gain and genetic diversity. The optimisation of this function results in OCS. The practical implementation of this method is based on an evolutionary algorithm, with weightings and constraints easily invoked to ensure practical relevance, precise control of the response to each trait in the economic index, and other requirements of progressive breeding programs. MateSel dictates which individuals to select and the actual mating allocations and/or selfings to be made, based on the breeder's identification of candidates for crossing and the maximum permissible number of crosses for each candidate.
MateSel was instructed to design 250 matings (to initiate cycles 3 and 4) and 150 matings (to initiate cycle 5) (Online Resource S2) from candidates based on a conservative balance strategy of target 45 degrees or higher, where 0 degrees gives full emphasis to short term genetic gain in index and 90 degrees gives full emphasis to minimizing parental coancestry [26]. Each genotype was limited to a maximum of 30 matings which was normally the maximum number of remnant self seeds. Remnant self seeds of selected genotypes were sown in the glasshouse as parent plants for crossing, based on the mating design derived from MateSel.

Assessment of Genetic Gain
Two types of genetic gain were assessed in this study -genetic gain in the breeding population over cycles, and genetic gain in historical cultivars included as controls in the trials over year of release.
The rate of genetic gain in the breeding population was estimated from the linear regression over cycles of PBV or OP of candidate genotypes in cycles 2, 3, and 4 from the optimum MMM-FA analysis which included all data from cycles 2, 3, and 4.
The rate of genetic gain in Australian historical cultivars was assessed by regressing the PBV or OP of historical cultivars based on the optimum MMM-FA analysis over their year of release. The Australian historical cultivars are listed in Section 2.5 'Field Trials and Phenotyping' (above).

Environmental Trends over Cycles and Countries
On average, 0.6% of plots at each trial site were invalid for various reasons and were excluded from the analysis of GY (Online Resource S6).
The predicted site mean GY from the MMM-FA base model varied in Australia and Canada from site to site and cycle to cycle with no significant environmental trend in GY over cycles ( trend, the average PlHt was higher in Australia (129.2 cm) than in Canada (103.1 cm) (Online Resources S6, S7). Average Oil was 46.9% in Australia and 40.4% in Canada, in contrast to average ProM which was 39.8% in Australia and 43.6% in Canada. The average GSL was 11.3 μmole g -1 , the average OL was 61.7%, the average SW100 was 0.325 g and the average BL score was 5.1 (Online Resources S6, S7). There were no significant environmental trends over cycles in site mean values of any trait except a very small upward trend in SW100 (Online Resources S6, S7).  Table 2). The slope of the linear regression of the predicted site mean GY across cycles was not significant (NS). For a description of trial codes, see Table 1.

Connectivity of Genotypes within and across Cycles
The connectivity of genotypes between trials within cycles was high (each S0,1 bulk was tested at multiple sites within cycles) but low across cycles (between 1 and 9 control cultivars were common across cycles) (Table 1b). However, the connectivity of the pedigree across cycles was high as a result of the use of both cross-and-self-sibs, that is, S2 selfprogeny (evaluated as S2,3 families) and S0 cross progeny (evaluated as S0,1 families), and this increased the accuracy of additive genetic values in the base model (Online Resource S6). A small portion of the pedigree shows this high level of connectivity between cycles  Table 2). The slope of the linear regression of the predicted site mean GY across cycles was not significant (NS). For a description of trial codes, see Table 1.
The average DTF was 104.1 days in Australia and 44.3 days in Canada, which reflects the different responses of spring canola to the winter-spring growing season in Australia and the summer growing season in Canada (Online Resources S6, S7). Following this trend, the average PlHt was higher in Australia (129.2 cm) than in Canada (103.1 cm) (Online Resources S6, S7). Average Oil was 46.9% in Australia and 40.4% in Canada, in contrast to average ProM which was 39.8% in Australia and 43.6% in Canada. The average GSL was 11.3 µmole g −1 , the average OL was 61.7%, the average SW100 was 0.325 g and the average BL score was 5.1 (Online Resources S6, S7). There were no significant environmental trends over cycles in site mean values of any trait except a very small upward trend in SW100 (Online Resources S6, S7).

Connectivity of Genotypes within and across Cycles
The connectivity of genotypes between trials within cycles was high (each S 0,1 bulk was tested at multiple sites within cycles) but low across cycles (between 1 and 9 control cultivars were common across cycles) (Table 1b). However, the connectivity of the pedigree across cycles was high as a result of the use of both cross-and-self-sibs, that is, S 2 selfprogeny (evaluated as S 2,3 families) and S 0 cross progeny (evaluated as S 0,1 families), and this increased the accuracy of additive genetic values in the base model (Online Resource S6). A small portion of the pedigree shows this high level of connectivity between cycles as a result of selfing and crossing of parent plants (Online Resource S1). The connectivity of genotypes between cycle 1 and the subsequent cycles was very low (Table 1b), and hence this cycle was excluded from further analysis.

Inbreeding Coefficients and Coefficients of Coancestry
The inbreeding coefficient (F) of individuals in the pedigree increased as the level of selfing increased and the coefficient of coancestry ( f ) between genotypes also increased as selfing increased in parent genotypes. For example, the F-values of individuals A1 and A2 in the pedigree diagram (Online Resource S1) were 0.0 and 0.5, respectively, as expected for F 1 and F 2 progeny of elite homozygous founder lines AU1 × EU1 (Online Resource S3). Founder lines AU1 and EU1 are assumed to be unrelated because the breeding programs in the southern hemisphere and northern hemisphere were mostly genetically isolated prior to 2000 [51]. Cycle 1 genotype C1 is an S 0 progeny of the cross between two F 1 genotypes A1 × B with F-value 0.019 (Online Resource S3), which reflects a small level of inbreeding between the two Australian founder lines AU1 and AU3 or between the European founder lines EU1 and EU3 prior to their use as founder parents. Self-progeny of C1 increased in F-value as the level of selfing increased from 0.509 in S 1 genotype C2, 0.755 in S 2 genotype C3, through to 0.969 in S 5 genotype C6 (Online Resources S1, S3).
The F-values in S 0 cross progeny varied greatly between crosses depending on the level of common ancestors in the parents, but independently of the F-values in the parent genotypes. Cycle 4 S 0 genotype Q (progeny of cross J4 × O3) had an F-value 0.326 due to several ancestors in common between parents J4 and O3 prior to cycle 4, but cycle 4 S 0 genotype P (progeny of cross C6 × O2) had a much lower F-value 0.077 as a result of very few common ancestors between parents C6 and O2, despite the fact that one parent C6 had a very high F-value of 0.969 (Online Resources S1, S3).
The coefficient of coancestry ( f ) among full-sib progeny increased as the selfing level in parent genotypes increased. For example, the f -value of cycle 1 full-sib progeny G1 and H1 was 0.262 (Online Resource S4), which is close to the expected value of 0.25 because the parent plants were F 1 genotypes and not inbred and not related. The f -value is slightly elevated due to a small level of inbreeding between the two Australian founder lines AU5 and AU6 or between the Canadian founder lines CA5 and CA6 prior to their use as founder parents (Online Resources S1, S4). However, the f -value of full-sibs J1 and K in Cycle 2 was higher than expected at 0.389 due to partial inbreeding in their S 1 parents E2 and H2 (Online Resource S4). J1 is a cross progeny and E3 is a self-progeny of their common parent, E2, and J1 and E3 are cross-and-self-sibs with an expected f -value equivalent to full-sibs. The actual f -value between J1 and E3 was 0.388 due to partial inbreeding in their common S 1 parent E2 (Online Resources S2, S4).
Coefficients of coancestry between offspring of different parents in the population varied greatly. The f -value between cycle 4 S 0 genotypes P and Q (0.262) was relatively high due to the presence of several common ancestors (especially O1 and J1), but the f -value between cycle 4 S 0 genotypes J1 and I (0.080) was low since they shared no recent common ancestors (Online Resources S2, S4).
In summary, coefficients of coancestry between genotypes were often higher than expected due to selfing in parent genotypes, and this was helpful to increase the accuracy of PBV. Likewise, the presence of cross-and-self-sibs in addition to full-sibs and many other collateral relatives helped to increase the accuracy of additive genetic values (Online Resource S6) and increased pedigree connections across cycles. High coefficients of coancestry induced by selfing within cycles did not necessarily contribute to population inbreeding because optimized crossing designs favoured matings between unrelated genotypes with OCS in MateSel [26].

Optimum MMM-FA Models
For GY, the optimum MMM-FA model was FA (2), where the percentage genetic variance accounted for (%VAF) was 68.3% and there was a trivial increase in %VAF in the FA(3) model (Table 3). FA(2) was the optimum model for most traits except SW100 and DTF where the optimum model was FA(1) ( Table 3). The average PBV for an individual across environments was generated from the optimum MMM-FA model for each trait at the end of cycle 4. Table 3. Evaluation of MMM-FA models of the genotype × environment interaction effects for each trait, and number of estimated variance components in each model. The optimum model a was selected on the basis that it increased percentage genetic variance accounted for (%VAF), and significantly increased REML log-likelihood (RLL) and decreased Akaike information criterion (AIC) and Bayesian information criterion (BIC). The base model included information from the additive genetic relationship matrix and significant fixed and random terms. Most of the %VAF occurred in FA(1) with minor increases in FA(2), and for most traits, the loadings for FA(1) were all positive (Online Resource S8). In this case, OP is a "generalised main effect which allows for heterogeneity of scale between environments" [22], and OP is a valid measure of the overall genetic value of a variety for most traits. The exception is SW100, where there were negative and positive loadings across sites in FA(1) as a result of crossover genotype × environment interaction (Online Resource S8).

Genetic Correlations across Trials
There were moderate to high genetic correlations of additive effects for GY across trials in Australia and Canada within and between cycles ( Figure 5) and high correlation coefficients across trials for all other traits (Online Resource S9). That is, trials in Canada and Australia identified similar high-performing and low-performing genotypes for all traits.  Table 2). For description of trial codes, see Table 1.

Pairwise Correlations of PBV across Traits
Cluster analysis of pair-wise correlations of PBV across traits from the optimum MMM-FA model revealed two main groups of traits ( Figure 6). GY, BL, Oil, PlHt, and DTF clustered together as a group of positively correlated traits, and this group tended to be negatively correlated with a second group of traits including GSL, ProM, OL, and SW100. DTF was strongly positively correlated with PlHt and late flowering lines tended to have smaller seeds ( Figure 6). GY, Oil, PlHt, and DTF were all positively correlated with BL, that is, BL-resistant genotypes tended to be high-yielding, later flowering, and taller with higher seed oil. GSL was negatively correlated with both Oil and GY ( Figure  6).  Table 2). For description of trial codes, see Table 1. Genetic correlations of GY between the disease nursery sites 2016AU2 and 2020AU2 and other sites in Australia or Canada were low to moderate, whereas genetic correlations of GY between non-disease sites were much higher ( Figure 5, Online Resource S9). This may result from some BL-susceptible genotypes performing relatively poorly for GY at disease nursery sites but relatively better at non-disease nursery sites, and vice versa.

Pairwise Correlations of PBV across Traits
Cluster analysis of pair-wise correlations of PBV across traits from the optimum MMM-FA model revealed two main groups of traits ( Figure 6). GY, BL, Oil, PlHt, and DTF clustered together as a group of positively correlated traits, and this group tended to be negatively correlated with a second group of traits including GSL, ProM, OL, and SW100. DTF was strongly positively correlated with PlHt and late flowering lines tended to have smaller seeds ( Figure 6). GY, Oil, PlHt, and DTF were all positively correlated with BL, that is, BL-resistant genotypes tended to be high-yielding, later flowering, and taller with higher seed oil. GSL was negatively correlated with both Oil and GY ( Figure 6). and 100 seed weight (SW100). Correlation coefficients are shown as values below the diagonal and as circles above the diagonal, where the relative size and colour of circles represents the size and sign of the correlation coefficients as indicated by the colour code on the right-hand side.

PBV, Overall Performance and Yield Stability for GY
PBV GY and OP GY from the optimum MMM-FA model were closely correlated across genotypes (r = 0.957), which confirms that both assess average performance of genotypes although the range of values was wider for PBV GY than for OP GY (Figure 7). RMSD values for GY (t ha -1 ) were perfectly correlated to the absolute value of f2 scores as expected from the FA(2) MMM-FA model for GY [21,22], although this may not necessarily be the case for FA(3) and higher-level models. RMSD is preferred over f2 scores since it is scaled in trait units and is applicable to FA models at all levels [22]. The traits correspond to grain yield (GY), days to 50% flower (DTF), plant height (PlHt), seed oil (Oil), protein in meal (ProM), glucosinolates (GSL), oleic acid (OL), phoma stem canker (blackleg) disease score (BL) and 100 seed weight (SW100). Correlation coefficients are shown as values below the diagonal and as circles above the diagonal, where the relative size and colour of circles represents the size and sign of the correlation coefficients as indicated by the colour code on the right-hand side.

PBV, Overall Performance and Yield Stability for GY
PBV GY and OP GY from the optimum MMM-FA model were closely correlated across genotypes (r = 0.957), which confirms that both assess average performance of genotypes although the range of values was wider for PBV GY than for OP GY (Figure 7). RMSD values for GY (t ha −1 ) were perfectly correlated to the absolute value of f 2 scores as expected from the FA(2) MMM-FA model for GY [21,22], although this may not necessarily be the case for FA(3) and higher-level models. RMSD is preferred over f 2 scores since it is scaled in trait units and is applicable to FA models at all levels [22].

Genetic Gain across Cycles
There was a very high rate of genetic gain in GY in the breeding population of 0.087 t ha −1 y −1 or 4.3% y −1 (expressed as a percentage of the population mean) over cycles 2, 3, and 4 from 2016 to 2020 as assessed by PBV for candidates for crossing in each cycle (Figure 8a, Table 4). Genetic gain assessed by OP resulted in a lower slope of 0.059 t ha −1 y −1 or 2.9% y −1 over cycles 2, 3, and 4 from 2016 to 2020, mostly due to lower cycle 3 and cycle 4 means in OP (Figure 8b) compared with PBV (Figure 8a).  Table 2). *** p < 0.001.

Genetic Gain across Cycles
There was a very high rate of genetic gain in GY in the breeding population of 0.087 t ha -1 y -1 or 4.3% y -1 (expressed as a percentage of the population mean) over cycles 2, 3, and 4 from 2016 to 2020 as assessed by PBV for candidates for crossing in each cycle (Figure 8a, Table 4). Genetic gain assessed by OP resulted in a lower slope of 0.059 t ha -1 y -1 or 2.9% y -1 over cycles 2, 3, and 4 from 2016 to 2020, mostly due to lower cycle 3 and cycle 4 means in OP (Figure 8b) compared with PBV (Figure 8a).
The significant genetic gain in GY over cycles (Figure 8) occurred in the absence of any environmental trend in site mean GY across cycles ( Figure 4). Therefore, the estimates of genetic gain in GY over cycles ( Figure 8) were not influenced by environmental trends and are true genetic trends.
There was a very high rate of genetic gain in the economic index of 7.26% y -1 (US$69.96 ha -1 y -1 ) as assessed by PBV and 5.05% y -1 (US$47.71 ha -1 y -1 ) as assessed by OP over cycles 2, 3 and 4 from 2016 to 2020 (Table 4, Online Resource S10).  Table 2). *** p < 0.001.   Table 4. Genetic gain in predicted breeding values (PBV) and overall performance (OP) across cycles 2, 3 and 4 for each trait and the economic index. PBV and OP were obtained from the optimum MMM-FA model of the genotype×environment interaction effects (see Table 2). The linear regression coefficient of candidate PBV and OP across cycles 2 to 4 is shown in units cycle −1 , and annual genetic gain from cycles 2 to 4 as units y −1 and as % of the population mean y −1 for each trait. Italicized values of coefficient and intercept in grey were not significant at p = 0.05. a Trait abbreviations: GY = grain yield, DTF = days to 50% flowering, PlHt = plant height, Oil = seed oil at 6% moisture, ProM = protein in meal at 10% moisture, GSL = glucosinolates, OL = oleic acid, BL = Phoma (blackleg) disease score from 1 (very susceptible) to 9 (very resistant), and SW100 = 100 seed weight. b Genetic gain in PBV and OP per year is calculated by dividing the linear regression coefficient across cycles by 2, and is expressed as units y −1 and as % of population mean y −1 . The significant genetic gain in GY over cycles (Figure 8) occurred in the absence of any environmental trend in site mean GY across cycles (Figure 4). Therefore, the estimates of genetic gain in GY over cycles ( Figure 8) were not influenced by environmental trends and are true genetic trends.
There was a very high rate of genetic gain in the economic index of 7.26% y −1 (US$69.96 ha −1 y −1 ) as assessed by PBV and 5.05% y −1 (US$47.71 ha −1 y −1 ) as assessed by OP over cycles 2, 3 and 4 from 2016 to 2020 (Table 4, Online Resource S10).
Significant genetic gain (or loss) in PBV and OP was found for all traits under selection in the economic index except PlHt (Table 4, Online Resource S10). There were no environmental trends in site means across cycles for any trait except for a very small environmental trend upwards in SW100 across cycles (Online Resource S7) which was opposed to the small genetic trend downwards in SW100 across cycles (Table 4). Therefore, the estimates of genetic trend in PBV over cycles for all traits are true genetic trends.
The economic weight in the selection index for increasing ProM was double that for increasing Oil (Table 2) in order to achieve the desired positive gains in both ProM (+0.167% y −1 ) and Oil (0.276% y −1 ) when assessed by PBV. The genetic gain in ProM assessed by OP (+0.187% y −1 ) was superior to the assessment by PBV, but the genetic gain in Oil by OP (0.249% y −1 ) was inferior to the assessment by PBV ( Table 4).
The average BL score in the population began at a low level in cycle 2 but increased rapidly by +8.3% y −1 as assessed by PBV, and +7.0% y −1 as assessed by OP (Table 4, Online Resource S10).

Genetic Gain in Historical Cultivars
The rate of genetic gain in GY in historical cultivars was 0.026 t ha −1 y −1 (1.3% y −1 ) over their year of release from 1996 to 2015 as assessed by PBV (Figure 9a), and 0.023 t ha −1 y −1 (1.1% y −1 ) as assessed by OP (Figure 9b). Significant genetic gain (or loss) in PBV and OP was found for all traits under selection in the economic index except PlHt (Table 4, Online Resource S10). There were no environmental trends in site means across cycles for any trait except for a very small environmental trend upwards in SW100 across cycles (Online Resource S7) which was opposed to the small genetic trend downwards in SW100 across cycles (Table 4). Therefore, the estimates of genetic trend in PBV over cycles for all traits are true genetic trends.
The economic weight in the selection index for increasing ProM was double that for increasing Oil (Table 2) in order to achieve the desired positive gains in both ProM (+0.167% y -1 ) and Oil (0.276% y -1 ) when assessed by PBV. The genetic gain in ProM assessed by OP (+0.187% y -1 ) was superior to the assessment by PBV, but the genetic gain in Oil by OP (0.249% y -1 ) was inferior to the assessment by PBV ( Table 4).
The average BL score in the population began at a low level in cycle 2 but increased rapidly by +8.3% y -1 as assessed by PBV, and +7.0% y -1 as assessed by OP (Table 4, Online Resource S10).

Genetic Gain in Historical Cultivars
The rate of genetic gain in GY in historical cultivars was 0.026 t ha -1 y -1 (1.3% y -1 ) over their year of release from 1996 to 2015 as assessed by PBV (Figure 9a), and 0.023 t ha -1 y -1 (1.1% y -1 ) as assessed by OP (Figure 9b).
Neither Oil nor ProM changed significantly over the year of release in historical cultivars, although the trend was towards higher Oil and lower ProM (Online Resource S10). There were no significant trends in BL, PlHt, and DTF over the year of release, but large variation between cultivars (Online Resource S10).
In contrast, the breeding population showed significant genetic gains in PBV over cycles for most traits, although in the case of GY and BL, the breeding population began below the average of historical cultivars in cycle 2 and exceeded historical cultivars by cycle 4 (Figures 8 and 9, Online Resource S10).

Impact of Genetic Correlations between Traits on Genetic Gain
The strong positive genetic correlations between BL, GY, DTF and PlHt resulted in the population becoming taller and later flowering in cycles 2 and 3 ( Figure 6, Online Neither Oil nor ProM changed significantly over the year of release in historical cultivars, although the trend was towards higher Oil and lower ProM (Online Resource S10). There were no significant trends in BL, PlHt, and DTF over the year of release, but large variation between cultivars (Online Resource S10).
In contrast, the breeding population showed significant genetic gains in PBV over cycles for most traits, although in the case of GY and BL, the breeding population began below the average of historical cultivars in cycle 2 and exceeded historical cultivars by cycle 4 (Figures 8 and 9, Online Resource S10).

Impact of Genetic Correlations between Traits on Genetic Gain
The strong positive genetic correlations between BL, GY, DTF and PlHt resulted in the population becoming taller and later flowering in cycles 2 and 3 ( Figure 6, Online Resource S10). This group of traits was negatively correlated with ProM, and therefore the economic weight for ProM was double that for Oil (Table 2) in order to achieve desired genetic gains in ProM and Oil and all other traits. It was also necessary to introduce a negative economic weight against PlHt and DTF in cycles 3 and 4 ( Table 2), and as a result, these began to fall in cycle 4 (Online Resource S10). The negative correlation between GY and GSL (Figure 6), coupled with the positive genetic gain in GY and Oil, contributed to a significant decrease in GSL over cycles (Online Resource S10).

Predictions from OCS in Cycle 5 with and without RMSD GY in the Economic Index
Two scenarios were assessed for the selection of parents and optimised mating design in MateSel at the end of cycle 4: an economic index excluding RMSD GY, and an economic index including a negative economic weight on RMSD GY ( Table 5). The goal of including RMSD GY in the index was to improve yield stability (that is, to decrease RMSD GY) in the population and to evaluate the impact RMSD GY on the genetic gain in GY and other traits.    In cycles 3 and 4, RMSD GY ranged widely from close to zero (stable across environments) to 0.40 t ha −1 (unstable across environments) with no signs of decreasing across cycles, while PBV GY in the population increased ( Figure 10). With a negative economic weight on RMSD GY in the index and an optimum mating design generated in OCS, the population mean RMSD GY was predicted to decrease from 0.160 t ha −1 in cycle 4 to 0.109 t ha −1 in cycle 5, equivalent to 32.2% reduction or 16.1% y −1 ( Table 5). As a result, the candidates selected for crossing in cycle 4 in MateSel tended to have lower RMSD GY (higher yield stability across environments) and higher PBV GY than non-selected candidates ( Figure 10).
The predicted increase in PBV GY from cycle 4 to cycle 5 was 0.164 t ha −1 or 4.1% y −1 relative to the population mean without RMSD GY in the index, and 0.170 t ha −1 or 4.2% y −1 with RMSD GY in the index (Table 5). That is, the inclusion of RMSD GY in the index had no negative impact on the genetic gain in GY or other traits except for DTF and PlHt, which had a slower response to selection for earlier and shorter types when RMSD GY was included in the index (Table 5). Interestingly, even without RMSD GY in the index, RMSD GY was predicted to decrease by 19.8% in cycle 5, which suggests that high GY is associated with a more stable GY when selection is based on OCS with optimum mating designs (Table 5).
There was a small increase in achieved parental coancestry in cycle 5 from 0.085 without RMSD GY to 0.087 when RMSD GY was included in the index (Table 5). This is a very low level of achieved parental coancestry in cycle 5. We conclude that the inclusion of RMSD GY in the index did not limit potential genetic gain for GY and other key traits in cycle 5.
The population mean PBV of both DTF and PlHt was predicted to continue to decline in cycle 5 (Table 5) in response to negative economic weights in the index (Table 2). Mean PBV BL in the population was predicted to increase only marginally in cycle 5 (Table 5, Online Resource S10).

Discussion
A very high rate of genetic gain in GY of 0.087 t ha −1 y −1 (4.3% y −1 ) was achieved inside a spring oilseed rape (canola) breeding program based on field testing in Australia and Canada over three cycles of S 0,1 family selection from 2016 to 2020, as assessed by average predicted breeding values (PBV) for GY in the optimum MMM-FA model (Figure 8a). The rate of gain assessed by overall performance (OP) based on f 1 scores in the optimum MMM-FA model and expressed in trait units was lower at 0.059 t ha −1 y −1 (2.9% y −1 ) (Figure 8b). This very high rate of genetic gain in GY in spring canola was the result of index selection of multiple economic traits with optimized mating designs based on OCS, and was predicted to continue with low rates of achieved parental coancestry in cycle 5 (Table 5). This high rate of genetic gain has the potential to continue into the foreseeable future as anticipated in stochastic models based on this breeding approach [9,10].
Both measures of genetic gain, PBV and OP, are valid and interesting outcomes of MMM-FA models, which help to explain the environmental impact on the assessment of genetic gain. Both PBV and OP support the conclusion that rapid genetic gain in spring oilseed rape was achieved across two disparate global regions (summer growing season in Canada; winter-spring growing season in Australia). All the evidence, including high genetic correlations of additive effects across trials for most traits, suggests that the genetic gains in this population were realised in both Australia and Canada.
The reasons for the lower rate of genetic gain assessed by OP vs. PBV in this study are complex but may reflect the relatively low loadings to FA(1) for GY in the optimum MMM-FA model allocated to disease nursery trials 2018AU2 and 2020AU2 in cycles 3 and 4 (Online Resource S8). At these disease nursery sites, BL-resistant candidates with relatively high GY may contribute more to PBV than to OP as a result of the low loadings to FA(1) for GY at these sites. Consequently, PBV would be higher than OP for the candidates in cycles 3 and 4 ( Figure 8).
The genetic gain was achieved in several agronomic and grain quality traits in these trials, including increases in Oil, ProM, and BL disease resistance and a decrease in GSL (Table 4). These estimates of genetic gain were not biased by environmental trends over years [19,20].
In historical cultivars in the same trials, the estimate of genetic gain over years of release was considerably less than inside the breeding population. Genetic gain in GY in Australian historical triazine tolerant cultivars was 0.026 t ha −1 y −1 (1.3% y −1 ) from 1996 to 2015 as assessed by PBV (Figure 9a), and 0.023 t ha −1 y −1 (1.1% y −1 ) as assessed by OP (Figure 9b). These values are close to the genetic gain in historical cultivars of lupins, which were 0.020 t ha −1 y −1 over 31 years assessed by PBV and 0.017 t ha −1 y −1 assessed by OP [21]. Similarly, the average achieved a genetic gain in GY in the world's major self-pollinating grain crops is approximately 1% y −1 [1].
High rates of genetic gain for GY and other traits inside this actively evolving crop breeding program in Australia and Canada verify the predictions of stochastic models of augmented S 0,1 family selection with OCS based on an economic index [9,10]. Also, for the first time in crop breeding, we incorporate a parameter related to yield stability (RMSD) [22] in the economic index in cycle 4, and RMSD GY was predicted to decrease (that is, GY stability to increase) in cycle 5 without negatively impacting the rate of genetic gain in GY (Table 5). This result, and the high genetic correlations (>0.80) in GY and all other traits across sites in Canada and Australia ( Figure 5, Online Resource S9), support the conclusion that yield stability across both local and global environments can be improved with the inclusion of RMSD in the economic index.
In addition to high rates of genetic gain in GY, rapid increases were achieved in Phoma stem canker (blackleg) disease resistance (+8.3% y −1 ), seed oil (+0.62% y −1 ), and protein in meals (+0.41% y −1 ), while seed glucosinolate content decreased (−3.4% y −1 ) (Table 4). However, genetic correlations between traits were both helpful and detrimental, and changes were made to the weightings of traits in the economic index across cycles as more knowledge became available to achieve desired gains. OCS in MateSel was important to evaluate the best weighting on traits in the index to achieve the desired predicted genetic gains in the next cycle, for example, to increase the weighting against tall and late flowering genotypes in the index for the selection of parents for cycle 5 ( Table 2). This resulted in a predicted decrease in the population mean DTF and PlHt in cycle 5 ( Table 5, Online Resource S10).
These high rates of genetic gain in GY and other traits in the economic index occurred with an achieved parental coancestry in cycle 5 parents of 0.087 with RMSD GY in the index (Table 5). This was lower than predicted in stochastic models of augmented S 0,1 family selection with OCS based on an economic index, where achieved parental coancestry was >0.20 with OCS after five cycles [10]. In this study, 35 new migrants (F 1 's of SH × NH crosses) were added in cycles 1, 2, 3, and 4, which helped to reduce achieved parental coancestry in cycle 5 (Online Resource S2). The low achieved parental coancestry in cycle 5 bodes well for future genetic gain in this breeding program, which should continue at similarly high rates for 20 or 30 cycles as predicted in stochastic models [9,10].
The breeding method outlined here follows the principles of BRIO crop breeding, based on accurate breeding values, rapid cycles, index selection, and OCS [52], and adherence to these principles was important to achieve rapid and sustainable improvements in several economic traits in this study. OCS or similar optimised mating systems are important to conserve genetic diversity while optimising long-term genetic gain [53,54]. BRIO principles are flexible to include new technology as it arises. For example, we applied these principles to the breeding of a common bean with genomic relationship information [32]. Genomic relationship information (G-BLUP) may be combined with additive relationship information (A-BLUP) in single-step H-BLUP analysis [31] which should further improve the accuracy of the PBV and help to accelerate genetic gain. Another attractive feature of BRIO principles is that selection goals can be adapted to changing circumstances over time, just as we adapted selection goals to react to undesirable PlHt and DTF in cycles 2 and 3.
In the future, we may include heat stress tolerance in the index to pre-empt the negative impact of global warming on GY [10].
The size and cost of this augmented S 0,1 family selection breeding program were relatively small-a total of 900 crosses were undertaken and 14,424 small field plots were sown and harvested in field plot trials at 9 locations in Australia and Canada from 2014 to 2020 ( Table 2). The cost of the field plots and grain quality testing is the same as in traditional plant breeding trials. The analyses we applied to augmented S 0,1 family selection, including MMM-FA analysis and OCS in MateSel, were carried out on a laptop computer. OCS in MateSel was used to improve GY, cooking time, and seed Fe and Zn content in a common bean breeding program in East Africa [32] and to improve different traits in Desi and Kabuli chickpea varieties [55]. Another positive feature of augmented S 0,1 family selection is that it delivers inbred lines for advanced testing as shown in the pedigree diagram (Online Resource S1) and thereby combines the population improvement phase with the product development phase of plant breeding [56].
Interestingly, we achieved only a small increase in the accuracy of additive genetic values in inbred lines over non-inbred genotypes of approximately 0.02 (Online Resource S6), which supports the use of non-inbred parents to reduce breeding cycle time [5]. Our results show that selfing to purity is not essential before selection and use as a parent and may slow the rate of genetic gain by extending the generation interval (L). One reason for the relatively high accuracy of additive genetic values in non-inbred lines in this study is that we carried forward both self and cross progeny into the next cycle, and this increased the number of collateral relatives within cycles and increased the connectivity of genetic relationships across cycles. Cross-and-self-sibs have a coefficient of coancestry (f ) equivalent to full-sibs which helps to increase the accuracy of PBV.
MMM-FA provides a means to assess yield stability (RMSD GY) and overall performance (OP) of genotypes in the presence of genotype by environment interactions [21,22]. OP is a useful concept when a diverse breeding population is evaluated over multiple sites and cycles, and genotype by environment interactions are relatively low compared to the generalised main effect of variety performance [22]. PBV GY values of genotypes in the optimised MMM-FA model were closely correlated to OP GY values in this study (r = 0.957) (Figure 7), and both PBV and OP are valid measures of genetic gain. RMSD GY in this study reflects global yield stability across multiple sites in Australia and Canada. RMSD GY was perfectly associated with f 2 scores for GY in the FA(2) model in this study but this may not be the case in higher-order FA models.
Future selection for adaption to global spring rapeseed environments may require economic weights on traits specific to regions (for example, there may be specific weights for important spring canola target regions such as Canada, Europe, China, India, and/or Australia), and this may be accommodated through 'multiple end-uses' in MateSel [26]. The population would be managed with multiple versions of the economic selection index for each target region and a global index applicable to all regions, with specific crossing programs targeted for each region.
The integration of multiple traits into an economic index follows the logic of Hazel and Lush [57] that selection on an "index of net desirability is much more efficient than selection for one trait at a time". However, index selection must be accompanied by OCS or a similar optimised mating system to optimize long-term genetic gain. OCS based on an economic index was essential in this study to achieve concurrent genetic improvements in GY, Oil, ProM, and BL while decreasing GSL and maintaining OL, SW100, DTF, and PlHt with low rates of population inbreeding. The breeder can change economic weightings in the economic index to achieve desired gains over cycles as more data become available or breeding goals change; for example, here we increased the index weightings against DTF and PlHt in cycle 4 (Table 3) which greatly reduced the predicted population mean of these two traits in cycle 5 (Table 5). We also increased the index weighting of ProM to double the economic value of Oil (Table 3) in order to achieve a genetic gain in both traits (Online Resource S10). This is clearly a subjective decision of the breeder since there are no current