Below we follow Maris et al. (2018)’s enlightening distinction between anticipatory and corroboratory predictions. Corroboratory predictions are hypothesis-derived predictions. Their goal is to be compared to observations for the purpose of understanding. They are necessary for corroborating or invalidating a hypothesis. In this paper, we analyze how quantitative genetic models of selection across two levels of diversity can provide operational anticipatory predictions. Anticipatory predictions consist in the application of extant knowledge. They can be used to achieve a transformative goal (Maris et al., 2018), i.e., the essence of plant breeding. Indeed, breeders usually have to act on a genetic system before knowing exactly how the system they are acting on will react (i.e., before knowing its response to selection). In this context, anticipatory predictions are expected to fuel a decision-making process by improving the intelligibility of the potential effects of breeders’ practices on complex, multispecies genetic systems. For this purpose, the three criteria mentioned in the previous section are tightly interlinked.

The breeder’s equation and its equivalent forms (Walsh and Lynch, 2018) predict the mean population change in a trait z from the product of a measure of genetic variation (such as heritability) by a measure of selection strength (such as a selection differential, S):

with the heritability h=2σG2σP2, i.e., the ratio of the (additive) genetic to total phenotypic variance of the measured trait. This simple formula is widely used by breeders and evolutionary biologists but under several equivalent forms (cf. Walsh and Lynch, 2018; Bijma, 2020). Breeders emphasize the selection intensity i¯ and the accuracy of selection h:

while the equivalent version by Lande (1979) is most commonly used in evolutionary biology:

with β the linear selection gradient, i.e., the slope of the linear regression of the relative fitness w on trait z (β = cov[z, w]/var[z]). This last version (eq. 3) can be extended to predict the change in the multivariate case through the concepts of G-matrix and multivariate selection gradient β (Lande, 1979):

or, under a developed matrix form for l traits:

with G_ii the genetic variance for the i^th trait and G_ij its genetic covariance with the j^th trait, and β_i the selection gradient for the i^th trait. By focusing on the case of artificial selection for plant breeding, we assume that selection is strong, and that β is controlled by the breeder, therefore known with high accuracy. Consequently below, the error due to the incidental influence of natural selection on potentially unmeasured traits (synthesis in Walsh and Lynch, 2018, Chapter 20) is treated as negligible.

In the following, depending on the approach described, we refer to each or all of the three equivalent formulations mentioned above (eqs. 1, 2, 3, or 4) to extend the effects of genetic evolution (in response to artificial selection) on a community-level variable. Below, this community-level variable is labeled c (refer to Box 1 for a definition and refer to Box 2 for biological examples). c potentially encompasses several components of community functions such as the total biomass yield of mixtures of perennial forage crops, or the yield of cereal grain when a cash crop and a companion forage legume are selected together, or any other ecological functions such as weed suppressing effects among grain legumes when a cereal species is intercropped for this purpose (refer to the example of faba bean in Box 2). As we will see below, c can be a community-level index including weighted performance components from several species and functions.

Our focus is on the links between the population and the community levels. Thus, we do not distinguish between the type of breeding scheme or methods, which are well summarized elsewhere (cf. Walsh and Lynch, 2018). We thus made our analytical framework as general as possible, so it can be used for different artificial (i.e., intentional) multispecies selection approaches based on randomized trials with available estimates of genotypic or phenotypic values for selection candidates to produce the next generation. It is therefore able to cover a wide a range of breeding strategies and tools depending on how the quantitative genetic parameters are estimated (e.g., mass, family-based, genomic selection). Our analytical framework can therefore be easily adapted to account for breeding cycle length, plant reproductive systems (self- and cross-pollinated crops), or any other parameters that affect breeding efficiency, for instance in simulation approaches (Bančič et al., 2021).

The most straightforward and intuitive approach to examine how genetic variation affects a variable c at the community level is to apply the standard quantitative genetic expression to the variable concerned:

This assumes that the value of c for the i^th community comprised of associated species (refer to definition in Box 1) will deviate from its mean μ, under the influence of the genotype g_Bi of a focal species B (Box 1) and under the influence of environmental effects e_ik (which might include genetic effects originating from other species not controlled for here, see below). This is the basic rationale in community genetics: genetically related individuals belonging to a focal species (also termed “foundation species” in the ecological literature) will favor similar patterns and similar values for the community-level variable (Whitham et al., 2003).

The value of c is treated as a “community phenotype” or a community-level variable (Box 1): the phenotype of one or several associated species with a genetic variance σg2 and “community heritability” Hc2=σg2/(σg2+σe2) (Fritz and Price, 1988; Whitham et al., 2006). We use uppercase H to denote the broad sense heritability usually estimated in community genetics. Community genetics usually investigates the effect of the genetic variation in a focal plant species on its associated arthropod community. The experimental settings typically consist in a common garden with randomized replicated clones of the focal species and measurements taken at the community level such as arthropod abundance or species richness (reviews in Haloin and Strauss, 2008; Genung et al., 2011; Tack et al., 2012).

In such biological systems, the variable c only considers the IIGE on the arthropod community traits, excluding DGE on the plant species. Community heritability is therefore a measure of IIGE variance (which would not be true if c was, e.g., the total biomass of the whole plant-arthropod system). Thus, we now refer to IIGE heritability, HI⁢I⁢G⁢E2. The notion of IIGE heritability as presented above presents an analogy with population-level heritability without assuming that communities can be selected as a whole (Collins, 2003). It is thus a matter of some debate. Could HI⁢I⁢G⁢E2 be predictive of the response to selection of the community trait according to the breeder’s equation, i.e., Δ⁢c¯=HI⁢I⁢G⁢E2⁢Sc ? (with S_c a hypothetical selection differential on c). Community ecologists who quantified this parameter warned against such interpretation (Whitham et al., 2006; Genung et al., 2011). Whitham et al. (2006) asserted that this approach does not imply that communities have a fitness in the wild or evolve as populations do. Rather, they called for an interpretation of HI⁢I⁢G⁢E2 as an integrative measure of the cascading effects of genes of the focal species at the community level. This means that HI⁢I⁢G⁢E2 should be taken as a measure of an association to estimate how much of the variance in the community variable results from genetic variation in the focal species.

The notion of IIGE heritability was designed to investigate the strength of IIGE of host plants on other taxa. To what extent could it be useful for breeders interested in improving plant communities? The first practical advantage over a more fine-grained approach (see below) is that it is trait free: it does not investigate the effect of the phenotype of the focal species on the community variable (Figure 2A) and instead quantifies IIGE as a whole, thus limiting the need for trait phenotyping. For estimating IIGE heritability, a breeder would first have to choose a focal species and control for its genetic variation within a standard range of experimental crop communities. HIIGE2 would provide insights into the expected efficiency of selection of the focal species with respect to the improvement in the community variable.

The only example we found of HI⁢I⁢G⁢E2 estimates in a crop community in the literature was provided by Maamouri et al. (2015). These authors grew 46 contrasted lucerne genotypes (Medicago sativa) with a grass (Festuca arundinacea) grown as a companion crop. While they found strong average broad-sense heritability for the direct genetic effect on lucerne biomass (H² = 0.76, range across sampling rounds: [0.64 – 0.83]), the broad-sense heritability on the grass biomass was much lower (HI⁢I⁢G⁢E2 = 0.05 [0.00 – 0.17]). This suggests that selection on lucerne biomass would have at most minor consequences for the grass biomass.

Heritability estimates of IIGE or at the community level are exposed to the same flaws as heritability estimates at the population level (Hansen et al., 2011), including high dependency on the study context and assumptions and computation preferences (Wilson, 2008; Firmat et al., 2017). When the community variable is on an arithmetic scale such as yield, this can be partly improved by computing the coefficient of genetic variation for IIGE (i.e., the IIGE standard deviation standardized by the mean value of c, refer to Hansen et al., 2011). Such mean standardized estimates could help to perform meaningful comparisons across breeding populations and testing sites to adjust community-level selection strategies.

Although in nature, the value of c is not directly associated with a value of fitness (genetically based selection does not act at the community level), this might be the case in a breeder’s field if artificial selection is performed among isolated and genetically controlled communities. In this case, HI⁢I⁢G⁢E2 estimates might become interpretable within the breeders’ equation. However, DGE and IIGE can sometimes be selected as a whole. This is the situation addressed by the breeder’s version of the joint phenotype approach described in the following section.

Summing the genetic contributions of several species in a community trait is an alternative approach to interfacing responses to selection between the population and the community-level. The simplest is to sum the genotypic values of each interacting species (Queller, 2014):

where g_Ai and g_Bj are the breeding values for individual i of species A and individual j of species B, respectively. According to this statistical expression, the joint phenotype c is defined by Queller and Strassmann (2018) as “a trait or outcome that can be affected, and potentially evolve, under the influence of two or more parties.” Both species are treated symmetrically and the joint phenotype results from the potential interaction between species, but the details of these interactions (values of interacting traits, complex feedbacks, etc.) are treated as a black box (Figure 2B). The evolutionary properties of joint phenotypes were investigated by Queller (2014) by applying the Robertson-Price equation (Frank, 2012). Queller (2014) showed that the change in c is predicted by the sum of the genetic variance for each species, respectively multiplied by their selective effects on c, i.e., the covariances σ(w_A, c) and σ(w_B, c). A non-zero covariance between fitness w and c indicates that natural selection acts on species’ traits that affect the community variable (Johnson et al., 2009). When the two covariances for each species are of opposite signs, a change in c results from an evolutionary conflict between the two species: an increase in fitness in one species is constrained by increased fitness in the other species.

If c is the amount of light captured within a plant community (Queller and Strassmann, 2018), an increase in size in species A will increase its competitive effect by negatively affecting the fitness of species B, and vice versa. The conflict can be resolved and c increased further if, for instance, one species evolves a strategy of shade tolerance, generating niche differentiation with respect to light. Recent studies in multispecies grasslands showed that niche differentiation and complementary resource use can evolve in a single generation (van Moorsel et al., 2018; Meilhac et al., 2020). As such, the concept of joint phenotype developed by Queller (2014) might provide a relevant theoretical framework to investigate the genetic determinants of such patterns. However, the analysis by Queller (2014) assumes regulation processes under natural selection in both species, which is of limited relevance for guiding artificial selection at the community level, which requires more specific modeling.

An analog version was coined for plant breeding many years ago. Wright (1985) proposed a selection model relying on the variance decomposition of a joint phenotype summing the “economic yield” of each species. This approach was inspired by factorial designs used to improve the parental population for breeding schemes for hybrid maize. The scale shift here relies on an analogy between the genetic combining ability of alleles (i.e., within a genome) from inbred lines and the ecological combining ability of genotypes (i.e., among distinct genomes) from different species.

Wright (1985) approach aimed at modeling genetic gain on the joint phenotype when selection is performed at the scale of the population of bispecific communities that differ in their genotypic composition. Fully dissecting the genetic variation in the joint phenotype across a population of genetically controlled artificial communities requires a factorial design including a large number of combinations of each genotype of species A sown with each genotype of species B. Accordingly, two sorts of variance components of c can be partitioned: (1) the average or general effect of a genotype of one species on the joint phenotypes c, i.e., the general mixing ability (GMA) including both DGE and IIGE of this genotype on c; (2) the specific effect resulting from genotypic interactions deviating from additivity and quantifying the specific mixing ability (SMA) of individual pairs of genotypes on c (for more details see: Annicchiarico et al., 2019; Sampoux et al., 2020). Put simply, the GMA component of variance includes both additive direct genetic effects and effects of additive ecological interactions on c (average IIGE of a genotype across the factorial), whereas the SMA component of variance quantifies the contribution of non-additive effects on c.

Today, 35 years after Wright’s publication, no prediction based on a full decomposition of the variance of c in a bispecific community has yet been made. This is likely due to the demanding experimental design required for such a decomposition. For two species and a minimum of 30 candidate genotypes each and three replicates per pair, the number of plots would be 3 × 30² = 2,700, corresponding to a huge single site experimental design beyond the reach of most breeding programs. With a third species, the required experiment becomes completely unrealistic (81,000 plots) (but refer to Haug et al., 2021).

Specifically targeting the GMA variance of interacting species (Hill, 1990) can substantially reduce experimental requirements. This involves parallel selection of the GMA of the two species, each with a specific (recurrent) selection process. For each of the two species, candidate genotypes of the selected species (let’s say species A) are sown in mixture with a mixture of representative genotypes from the other (unselected) species (let’s say species B), used in the mixture as a tester for genotypes of A. In parallel, the same procedure is performed for species B with a mixture of representative genotypes of species A used as testers (i.e., selection in parallel for GMA, refer to Figure 2 in Sampoux et al., 2020). According to Sampoux et al. (2020), when targeting GMA variance for c by selecting on, e.g., species A, the measurable component of the joint phenotype c is as follows:

with μA and μB the mean value of species A and B contributing to the community value c. Here, the controlled genetic variance of c in the design is only due to species A, with vAi and aAi the direct genetic (DGE) and the indirect genetic effects (i.e., IIGE, a for “associated effect” in breeding terms) of species A, respectively. Together, these terms model the GMA of the ith genotype of species A. From this expression, improving c by selecting on species A only gives the following expected response for the joint phenotype:

i¯c is the selection intensity on the value of c and σ_c the total phenotypic variance for c. Under the procedure of selection in parallel for GMA described above, the analogous expression can be derived for species B (with A as a tester), also contributing to genetic gain for c. The value between brackets is the variance of the joint phenotype caused by genetic variation only in species A on which selection is performed. This variance expression is the sum of two components: the direct response to selection on c for the contribution of species A (σvA2+σvA,aA) and the correlated response on the value of species B (σvA,aA+σaA2) [for a detailed argument including a release from the assumption of the equal species weights, refer to Sampoux et al. (2020)]. Both components include the covariance between the DGE and the IIGE of the selected species. Interestingly, eq. 8 parallels the results obtained by Griffing (1967) for the response to group selection at the intraspecific level, when assuming unrelated interacting individuals.

The main advantage of this joint phenotype approach is that it emphasizes the constraining role of this covariance term for community-scale genetic improvement. If σ_vA,aA < 0, improving the contribution of species A (DGE) to the mixture while not accounting for its IIGE on B (a_A) is not possible without weakening the contribution to the mixture of species B (Wright, 1985). This typically happens when selection for yield in A increases its competitive effect, thus weakening the contribution of B. Breeders that has to cope with an evolutionary conflict (sensu Queller, 2014) and manage a trade-off between species vs. community performance. Performing two selection processes in parallel (one for each species), as proposed by Sampoux et al. (2020) and described above, could be a promising way to deal with such an evolutionary conflict.

By allowing integrative decomposition of the genetic interactions underlying the variation in a community variable, the joint phenotype approach provides quantitative genetic expression of the among-species conflicts the plant breeder will have to deal with. However, this approach enables prediction across biological scales (from species quantitative traits to c) through estimates of covariances (i.e., σ_vA,aA) between traits measured in two different species (as the IIGE of species A affects, e.g., the yield of species B). The nature of this parameter exposes quantitative predictions to strong limitations of both biological and methodological origins, as we will see below (refer to section headed “Common limitations and links between the three approaches”).

In contrast to the two previous approaches, this approach relies on the functional trait assumption: a set of functional traits, and not genotypes, affects the community variable (Figure 2C). Investigating how a set of traits measured in a single species affects community functioning requires a multivariate approach linking the evolution of genetically correlated functional traits to a community-level phenotype. The multivariate selection model proposed by Johnson et al. (2009) relies on the notion of a “community-trait gradient”: a multiple regression of a single community trait c on a set of l evolving traits for a focal species. Variation in the community variable is therefore modeled as follows:

where α_j is the partial regression coefficients for the j^th trait and g_ij its genetic value for the i^th genotype at trait j. The estimated vector of coefficients makes it possible to project the multivariate evolutionary change in a focal species on the ecological variable. This is a way of modeling the IIGE of a complex, multivariate phenotype on a predefined community variable. This “matrix projection model” therefore builds on Lande’s (1979) formulation of the breeder’s equation to extend it to a community variable (Johnson et al., 2009):

or, for a more synthetic notation:

The second and last terms are the elements of the standard multivariate breeder’s equation, respectively, the G-matrix, and the linear gradient of selection (Lande and Arnold, 1983), i.e., the vector of partial regression coefficients of each j trait on fitness, whose product gives the expected response to selection in the vector of traits. The first term is the vector of partial regression coefficients (estimated independently in eq. 9), allowing the projection the genetic changes on the community variable (the subscript T denotes the transpose of this vector). This projection of evolutionary change corresponds to the sum of the evolutionary changes for the l traits, weighted by the αj coefficients (refer to Appendix 1, for a derivation of this expression using the Robertson theorem). In ecological terms, these coefficients represent the level to which each j trait is “functional” (refer to Litrico and Violle, 2015), i.e., interact with the species of the targeted community and affect the community-level variable.

Johnson et al. (2009) originally used this model to investigate the effect of functional trait variation in a wild species Oenothera biennis (Onagraceae) on the associated arthropod community. The model described above helped these authors formalize a set of conditions for the evolution of IIGE on the focal plant species: the trait j should be genetically variable (G_jj > 0), selection should affect the trait (β_j ≠ 0), and the trait should cause variation in the community variable (α_j ≠ 0).

Johnson et al. (2009)’s approach is strictly analogous to a linear selection index: it weights the sum of change in each dimension of a multivariate phenotypes to project it on a single variable, i.e., the selection index. The only – but nevertheless significant – difference is that the coefficients of the index are estimated statistically (eq. 9, refer to Appendix 1). Selection indexes are widely used by breeders and their properties are well known (Hazel et al., 1994). According to the properties of the linear selection index (Lin and Allaire, 1977; Nordskog, 1978) and using Johnson et al. (2009)’s notation, the heritability of the community variable described by k underlying traits is given by:

with P the phenotypic variance–covariance matrix among traits. Thus, the response to selection in equation (10b) can be expressed in the selection intensity form, which is more familiar to breeders:

with i¯c the univariate intensity of selection applied to the community variable.

From a practical point of view, this approach includes two independent steps. The first step aims to estimate how the phenotypic trait values of the focal species affect the community variable (i.e., estimating α). The second step aims to quantify the genetic architecture of the ecologically relevant traits previously identified (i.e., estimating G).

This approach has two potential advantages. First, by explicitly modeling ecological functions in the first step, it avoids the potentially problematic use of covariances between IIGE and DGE, which play a central role in joint phenotype approaches, as this covariance is suspected of strong variability among experimental contexts (discussed in the following section). Second, the main practical advantage of this trait-centered approach is that it makes it possible to discard traits that have no influence on the target community variable c and to focus on the traits that have the main effects. From a multiple traits-community variable regression (eq. 9), a standard model selection procedure using information theoretical approaches (review in Grueber et al., 2011) could be used in the first step to discard traits that do not deserve to be taken into account for the improvement of c. This would make it possible to considerably reduce the size of the sample required for further quantitative genetic predictions in the following prediction step (eq. 10). Estimates of α could rely on biological and experimental knowledge from functional (agro)ecology (Garnier and Navas, 2012). The analysis of α would provide an innovative way to design for the focal species what we call “community-ideotypes,” i.e., ideotypes that not only target the focal species’ own performance but also the expected performance generated through IIGE by such an ideotype at the scale of the community. However, one should keep in mind that a trait j with no community-level effect (α_j ≈ 0) may still be useful to define the selection index due to its genetic correlation with other ecologically relevant traits.

Another key advantage of this approach is that it allows breeders to map the consequences of genetic constraints (i.e., genetic covariances among the focal species’ traits) in terms of response at the community scale. Let us assume a simple situation where the genetic covariance among two traits is positive (G₁₂ > 0) and both α₁ and α₂ for these traits are also positive. Then, selecting for trait 1 alone is expected to improve the contribution of both traits to the community performance (i.e., as Δ⁢c¯=α1⁢G11⁢β1+α2⁢G12⁢β1). But if α₁ and α₂ are of opposite signs (α₁ > 0, denoting, e.g., a direct genetic effect on total yield and α₂ < 0, denoting, e.g., a negative competitive effect on a companion legume that supplies nitrogen), selecting for trait 1 alone is expected to reduce the value of trait 2 due to the positive genetic covariance. In other words, this could enable breeders to easily objectivize the consequences of their usual (i.e., species-centered) work at the community scale, by generating anticipatory predictions at this level. For instance, by extending the standard breeding indices based on a series of traits measured in a pure stand to predict performance in a mixed stand (e.g., Annicchiarico, 2003), this approach would make it possible to predict the consequences of this species-specific performance at the community scale. This includes the quantitative genetic exploration of the promising notion of “biological interaction function” coined by Haug et al. (2021). To sum up, Haug et al. (2021) suggest correlating both direct and indirect genetic effects with measured traits that generate different types of species interactions in the mixture, leading to the identification of traits that produce favorable biological interactions at the community scale. To identify suitable cultivars for cereal–legume mixtures, Kammoun et al. (2021) argue for a simple statistical approach that predicts reliable mixtures based on pure stand performances for both components combined with a description of their interaction function. We believe that such a basic set of variables would be a reliable starting point for designing selection criteria for this type of simple crop mixture.

This approach is based on the predictions of c from genetic (co)variation in a single species. It therefore assumes that one species needs to receive more attention than any other species from the breeder. We observed that this is the most frequently reported case in the plant breeding literature (for recent examples: Annicchiarico et al., 2021; Ergon and Bakken, 2022; Moutier et al., 2022).

The three approaches we identified are distinct in several conceptual, statistical, and practical dimensions (Table 1). These approaches rely on distinct strategies for modeling causal pathways between genetic variation and community-level variation (Figure 2). However, a general model linking all three approaches is possible. The community-trait genetic gradient approach could theoretically be extended to more than one species. The evolution of c affected by two species A and B represented by l traits each (kept equal for the sake of simplicity) can also be modeled as the evolution of a linear index:

This equation simply sums the ecologically weighted effects of each evolving trait for each species. But note that changes in traits in species A and B are not independent, but are linked through genetic covariance across each pair of species traits, caused by genetic interactions among species (IIGE) mediated by interaction traits (sensu Litrico and Violle, 2015). The sum of the trait’s specific α for each species corresponds to their respective weights regarding their contribution to c. The generalization of equation (13) for m interacting species with l traits each is:

with the first and second sigma summing effects over species and over traits, respectively. Interestingly, obvious links exist between this general expression and the integrative functional parameters used to describe the trait-based diversity of functional types in a plant community (refer to Violle et al., 2007; Garnier and Navas, 2012). The issue, however, is that full predictions of selection response for such an index require a complex set of interacting parameters involving both species-specific G-matrices and reciprocal cross-species genetic covariance matrices (CSG), describing the interspecific, among trait genetic relations between DGE and IIGE between pairs of species (Figure 3). To our knowledge, the only estimate for such a combined G-CSG-matrix, i.e., including both the direct and indirect genetic effects for a focal species, was provided by Riedel et al. (2018). However, these authors did not report estimate uncertainties and did not compute the effect on an integrative, index-like, community-level variable.

A CSG-matrix combines the covariance between the DGE of a given species A with its IIGE (originating from A) on the partner species’ traits. It thus involves the association between traits measured in two different species (e.g., a biomass trait measured in a forage crop species A [DGE], that through IIGE, affects a grain yield trait measured in an associated cash crop species B), i.e., in totally unrelated individuals. This type of covariances is therefore distinct from standard genetic covariances among traits measured in the same series of individuals as well as for DGE and IGE measured in interacting individuals of the same species (Griffing, 1967). These genetic covariances among species’ traits are complex parameters that include both genetic and ecological (i.e., species interactions) causal pathways. They are pivotal parameters in the joint phenotype approach (Wright, 1985; Sampoux et al., 2020) and in some versions of multitrait approaches to community genetics (Riedel et al., 2018) as they make it possible to link population genetics to the community scale.

If the capacity of the general equation cited above (eq. 14) to provide reliable predictions should be considered with great caution, it is nevertheless useful for clarifying the statistical and conceptual links between the three approaches described above. With l measured traits in a single species (m = 1) eq. (14) reduced to Johnson et al.’s. (2009) community-trait gradient approach (eqs. 10a, 10b, and 12). With two species (m = 2), each with a single measured trait (l = 1) and no estimated ecological gradient (α_A = α_B = 1), this becomes Wright’s (1985) style joint phenotype approach (eq. 8). If the latter assumption is further reduced to one species, then the last remaining parameter is a variance component analogous to IIGE H² approaches. This is also equivalent to a joint phenotype approach with only the IIGE effect measured (i.e., in eq. 6, only σaA2 is assumed to be non-zero), underlining the fact that community H² approaches are a special case of the joint phenotype, with neglected direct effects. Breeding based on the community H² approach would then be possible if the community is taken as the unit of selection (i.e., the breeder assigns a fitness value – whether or not the genotype is selected – to the community-level performance of a genotype).