Haplotype sequences for both faba bean and triticale founder populations were simulated using Markovian Coalescent Simulator (MaCS) implemented in AlphaSimR, recreating the evolutionary process with multiple cycles of drift, mutation and selection (Chen et al., 2009). Notably, while our simulations were specifically conducted with faba bean and triticale, it’s important to emphasize that the framework we developed is adaptable and applicable to other crops as well. In total, 100 homozygous individuals of each crop were simulated to form the founder population. For computational efficiency, the genome of both component crops, triticale (an allo-polyploid) and faba bean (a diploid), has been simulated as a single diploid chromosome. For triticale, the physical length of the single chromosome was assumed to be 9.9e+08 basepair (bp), which is the average physical size of triticale’s wheat (Avni et al., 2017) and rye (Li et al., 2021) chromosomes. For the faba bean, the physical length of the single chromosome was assumed to be 1.1e+9 bp (Jayakodi et al., 2023).

In each scenario and for the two component crops, two correlated traits were simulated, representing the monocrop (MC) and intercrop (IC) performance of the breeding pool accessions. These traits were assumed to adhere to an additive genetic model that is controlled by 1000 biallelic QTL at polymorphic sites that are randomly distributed over the single chromosome and for which a reference allele effect is sampled from the standard normal distribution (Gaynor et al., 2021).

Genetic values (GVs) for both the MC and IC trait are computed through the summation of the QTL effects across all segregating sites. The loci influencing MC and IC trait performance are assumed to be identical but the different allele dosage effects are sampled in a way that fixes the genetic correlation between the two traits to one of three values from the set {0.3, 0.5, 0.9}. It should be clear that a positive correlation is required to allow for indirect selection for IC performance through evaluation of the MC response. The GV of an IC combination is obtained by averaging the IC GV of the faba bean and the triticale component as shown in Equation 1.

Where, g represents the GV of the IC combination, g_f denotes the IC GV of the faba bean component and g_t refers to the IC GV of the triticale component.

Phenotypes were simulated by adding a random environmental effect to the GVs of the MC and IC traits. These random environmental effects have been sampled from a normal distribution with zero mean and a variance σ e 2 r that matches the heritability of the trait and the number of replications r at a particular selection stage following Equation 2.

Where, h ² is the narrow sense heritability, σ a 2 the additive genetic variance of the founder population, σ e 2 the environmental variance.

This study focuses on two distinct MC breeding strategies for self-pollinating crops used as pure lines:

Figure 1 provides a schematic overview of the various breeding schemes that are examined in this study and have been partitioned according to the selection stage at which IC testing is introduced. In this breeding program, we assume there are three key selection stages: a preliminary yield trial (PYT), an advanced yield trial (AYT) and an elite yield trial (EYT). The PYT represents the initial screening stage, where a large number of lines are phenotypically evaluated without replication. AYT is the second stage, where promising lines from PYT are tested in two replications. Finally, the EYT represents the last phase, during which the top-performing lines from the AYT are assessed for potential commercial release. The presented breeding schemes have been specifically designed to have identical resource requirements for phenotyping in each of the selection stages. This implies that all schemes are constrained in the same way with respect to the number of candidate accessions that can be tested, enabling a fair comparison in terms of selection efficiency. Breeding schemes do differ in the number of initial crosses and the number of selected candidates at the end of each stage. The initial phase involves the creation of 14 initial crosses for both faba bean (Fb) and triticale (Tr). To facilitate these crosses, 14 parents are randomly sampled from the founder population of 100 individuals for each crop. In the EPR-based selection schemes, the 14 S0 plants are consecutively selfed until S4 assuming a visual, within family selection with low heritability (h² = 0.1). Subsequently, 8 percent of the individuals, 28 in total, are advanced to the preliminary yield trials. In the DH-based breeding schemes, the S0 individuals are assumed to produce fully homozygous inbred lines (DH0) without selection. In the MC scenario, 392 Fb and 392 Tr S4/DH0 plants enter PYT hosting a total of 784 single replicate plots on a single location. In the scenario where IC testing is introduced at PYT, 28 fb lines are combined with 28 tr lines, giving rise to 784 unique Fb x Tr line combinations.

Following the PYT, AYT allows to test 144 entries in two replicates, involving either 12 Fb x 12 Tr line combinations in the IC scenario or 72 Fb and 72 Tr lines in the MC setting. In the following EYT, there is room to test 36 IC combinations from 6 Fb and 6 Tr lines in four replicates. The lines that make it to the AYT stage are recycled as parents for making new initial crosses in the subsequent breeding cycle. In each scenario, a total of 20 breeding cycles are simulated to monitor genetic progress, derived from the IC GVs during EYT stage, involving 36 IC combinations.

The realized genetic progress of each breeding scheme is quantified by calculating the average IC GVs of the 36 IC combinations that are tested at the EYT stage for each of the 20 consecutive breeding cycles. All reported estimates are averaged over 100 independent iterations of the simulation routine. Figures 2 , 3 specifies the IC GVs for three levels of trait heritability (0.3, 0.5, 0.9) and three levels of genetic correlation between MC and IC trait performance (R = 0.3, R = 0.5, R = 0.9). Phenotypic evaluation of IC performance is either introduced in the PYT, AYT or EYT selection stages, as described in Figure 1 .

In all simulated scenarios, GVs increase over the initial breeding cycles and converge asymptotically to a maximum value. Convergence is assumed when the difference between the average GVs of the IC trait between two consecutive breeding cycles becomes smaller than 0.005 ( Supplementary Material Table 1 ).

Genetic variances in IC trait diminishes consistently across all scenarios, showing a gradual decline as breeding cycles advances ( Supplementary Material 2.2 ).

This asymptotic behavior of GV is readily explained by the gradual depletion of genetic variance over breeding cycles, as these schemes do not implement a variance preserving strategy like those advocated by Allier et al., 2020 and Vanavermaete et al., 2021. For this particular study, however, both the initial rate of genetic progress and the asymptotic limit are informative criteria to compare the efficiency of the examined selection schemes under different assumptions of trait heritability and genetic correlation between MC and IC trait performance.

IC breeding presents challenges due to the complicated interactions and involvement of various species (Thierfelder et al., 2012; Brooker et al., 2015; Bourke et al., 2021; Huss et al., 2022; Weih et al., 2022b). To facilitate this process, a simplified and customizable breeding strategy is required, one that can be adapted for a wide range of crops and is straightforward to implement. Numerous publications underscore the critical need for an adaptable breeding strategy to enable effective IC breeding (Annicchiarico et al., 2019; Bourke et al., 2021; Weih et al., 2022a; Moore et al., 2023). Our research aimed to determine the most suitable selection stage for transitioning from MC breeding to IC breeding. All examined breeding schemes assume phenotypic selection with a fixed number of available field plots at each selection stage, maintaining a constant phenotyping workload during both MC and IC settings. The presented approach facilitates the integration of IC trait selection into MC breeding programs without significant increases in resource allocation, suggesting a methodologically sound and economically viable pathway for enhancing crop yields and sustainability through IC breeding.

The relationship between MC and IC performance, specifically in terms of trait correlations, has been a consistent focus, with research over the decades affirming a existence of variable correlation among crops in MC and IC setting (Atwood and Garber, 1942; Dijkstra and De Vos, 1972; Gomez and Gomez, 1983; Caradus et al., 1989; de Oliveira Zimmermann, 1997; Holland and Brummer, 1999; Santalla et al., 2001; Gebeyehu et al., 2006; Maamouri et al., 2017; Homulle et al., 2021; MacLaren et al., 2023). Gomez and Gomez (1983) reviewed 30 studies comparing yields of MC and IC systems for cereals, legumes, and sweet potatoes. For cereals, correlations in yield variations across 28 varieties ranged from 0.35 to 0.90. Legume trials with 19 varieties, showed a wider correlation range from -0.36 to 0.91, these findings underline the critical role of trait correlations in IC performance. In this study, the influence of the genetic correlations between MC and IC performance is assessed in DH and EPR breeding schemes for a range of trait heritabilities, replicating real-world IC traits.

In scenarios characterized by a high genetic correlation (R = 0.9), our findings consistently advocate for the initiation of IC phenotyping at the EYT stage for traits spanning the full heritability range. This conclusion is readily explained by the observation that selection for MC performance in the PYT and AYT stages allows to test a much larger set of candidates while the high genetic correlation between MC and IC traits guarantees an indirect selection for IC performance. Furthermore, delaying IC phenotyping efforts until the EYT stage implies a considerable simplification of the phenotypic evaluation activities, making it easier to integrate IC breeding into existing MC programs. Results from previous studies also suggest that the calculation of a correlation coefficient between MC and IC can give an indication whether separate breeding efforts are necessary, if correlation is high and significant, there probably is no need for a separate breeding effort for IC improvement (Gomez and Gomez, 1983; Francis and Smith, 1985).

For traits that show a moderate genetic correlation (R = 0.5), our results indicate that IC phenotyping should be introduced at the AYT stage or later for traits of low to medium heritability. The findings of previous experiments suggest that there exists a moderate correlation between the MC performance of plant cultivars and their performance when grown in IC environment. This observation was consistently noted across different studies, particularly when pure stand performances were documented (Atwood and Garber, 1942; Caradus et al., 1989). Conversely, for traits with high heritability, IC phenotyping should be initiated at the EYT stage to maximize genetic progress. High heritability implies that the variation in these traits is largely due to genetic differences rather than environmental factors. A strong genetic signal makes it possible to delay intensive IC phenotyping efforts until the EYT stage, where the focus shifts to refining and selecting the top-performing lines. By this stage, breeding programs have already filtered out lines based on preliminary assessments, allowing resources to be concentrated on the most promising IC combinations. Low and medium heritability traits require more IC phenotyping efforts to identify the superior combinations.

The results of this study illustrate that for traits characterized by a low genetic correlation (R = 0.3), it is imperative to initiate IC phenotyping at the PYT stage, irrespective of trait heritabilities. or breeding method (i.e. DH or EPR). Given the low genetic correlation between MC and IC genotype performance, it becomes clear that selecting for MC performance does not allow to drive improvement in the IC trait. This is particularly challenging for breeding programs because genotypes that perform well in MC conditions might not exhibit the same level of performance or adaptability in intercropping setups. This discrepancy necessitates early initiation of IC phenotyping.

In this situation, postponing IC phenotyping until the EYT stage limits the opportunity to accurately select candidate genotypes for traits that are beneficial in intercropping systems early in the breeding process. Since genetic progress for IC traits is slow due to the low correlation with MC traits, waiting until the later stages of breeding to focus on IC performance risks overlooking lines that could have shown promising IC traits earlier. This delay can hinder the ability to achieve the maximum genetic potential for these traits, affecting the overall efficiency and effectiveness of breeding program.

In our study, we considered specific correlation level between MC and IC traits, however, it is essential to recognize that real world scenarios introduce additional complexities. For instance, a study conducted in Colombia examined the performance of climbing bean cultivars across two consecutive seasons, revealing a high correlation between yields in both systems during the first season (R = 0.82), while the correlation was lower in the subsequent season (R = 0.41). Similarly, yields of bush bean cultivars exhibited a high correlation in the initial season (R = 0.88) and a lower correlation in the following season (R = 0.51). Additionally, mungbean data from the International Rice Research Institute (IRRI) in the Philippines depicted similar results (Francis and Sanders, 1978; Francis and Smith, 1985).

Intercrop breeding for perennials faces unique challenges, notably due to their longer lifecycles and complex root systems, demanding innovative but patient strategies for effective selection and optimization. A recent review by Moore et al., 2022 highlights distinct objectives across various intercropping systems, focusing on optimizing productivity within specific forage mixtures, breeders aim to enhance the overall productivity of grass and legume mixes harvested together. Perennial grain-forage systems prioritize grain yield, with forage as a secondary aim. Strategies vary across systems, emphasizing temporal and spatial niche differentiation to reduce competition and enhance complementarity. These findings underscore the dynamic nature of intercropping systems and the importance of considering variability across seasons and cultivars in optimizing IC breeding schemes.

The suggested adaptations to existing MC breeding schemes should be considered as a first, low-threshold step towards dedicated IC breeding strategies. The next step likely involves the use of tester varieties to reduce the number of IC combinations that are evaluated in the field. In this setup, all selection candidates of one component crop are intercropped with a single accession of the other component crop. In hybrid breeding programs, this single accession is referred to as the tester line which is specifically chosen for its superior combining ability when crossed with accessions of the complementary heterotic group. This scheme translates to a selection for general mixing ability in an IC setting which is expected to increase IC selection efficiency.

While discussing the implications of our findings, it is critical to acknowledge the characteristic limitations of stochastic simulations. While these routines are instrumental for the theoretical optimization of breeding schemes, they inevitably simplify the complicated genetic interactions and environmental variables that characterize real-world scenarios.

Our simulations, predicated on various assumptions, may not encapsulate the entire spectrum of gene-environment interaction effects, particularly in the dynamic setting of intercropping. Although phenotype-based selection is employed as a realistic stand-in for assessing genetic potential, real-life phenotypes can be influenced by genotype-by-environment interactions which have not been considered in this study. Furthermore, most breeders are confronted with multiple-trait objectives that cover a range of heritabilities and as such obfuscate the optimal stage for intercrop (IC) testing. However, the commonly used multi-trait selection indexes allow to aggregate multiple trait objectives in a univariate pseudo-trait that is generally cursed with a low heritability but nonetheless adheres to the breeder’s equation and is therefore covered by the presented simulation scenarios. A similar logic can be followed for other commonly used multi-trait selection strategies such as independent culling or tandem selection that also represent specific aggregations on the involved traits.

It should also be noted that the cost and time per breeding cycle can differ substantially between selection strategies such as EPR and DH. This study, however, assumes that both approaches incur identical costs prior to field testing. However, the cost associated with DH breeding largely depends on the specificities of the applied DH protocol and the cost of EPR is a function of the number of selfing generations and the applied protocol for phenotypic evaluation prior to field testing. A more in-depth comparison of the EPR and DH selection strategies would require detailed cost estimates of the involved pre-testing steps for the crop combination under study.

Real-world breeding programs are also shaped by a confluence of factors beyond genetic data, including ecological considerations, economic viability and the decision-makers’ expertise, which are beyond the scope of our simulation models. Consequently, while our simulation-based approach yields valuable insights for prospective IC breeders, these insights must be corroborated through empirical research to ensure their practical relevance and accuracy in actual breeding programs.