Multidimensional (Co)Evolutionary Stability

The complexity of biotic and abiotic environmental conditions is such that the fitness of individuals is likely to depend on multiple traits. Using a synthetic framework of phenotypic evolution that draws from adaptive dynamics and quantitative genetics approaches, we explore how the number of traits under selection influences convergence stability and evolutionary stability in models for coevolution in multidimensional phenotype spaces. Our results allow us to identify three different effects of trait dimensionality on stability. First are (i) a “combinatorial effect”: without epistasis and genetic correlations, a higher number of trait dimensions offers more opportunities for equilibria to be unstable; and (ii) epistatic interactions, that is, fitness interactions between traits, which tend to destabilize evolutionary equilibria; this effect increases with the dimension of phenotype space. These first two effects influence both convergence stability and evolutionary stability, while (iii) genetic correlations (due, e.g., to pleiotropy or linkage disequilibrium) can affect only convergence stability. We illustrate the general prediction that increased dimensionality destabilizes evolutionary equilibria using examples drawn from well-studied classical models of frequency-dependent competition for resources, adaptation to a spatially heterogeneous environment, and antagonistic coevolution. In addition, our analyses show that increased dimensionality can favor diversification, for example, in the form of local adaptation, as well as evolutionary escape.


Introduction
In nature, the fitness of an individual will generally depend on local abiotic conditions, on interactions with other individuals from the same species, and on interactions with other individuals from different species. The complexity of both abiotic and biotic conditions, as well as the complexity of individual organisms themselves, make it very likely that many different individual traits will affect fitness (Blows 2007).
Even though models of quantitative trait evolution often consider the evolution of a single trait in a single species, seminal formulations of both evolutionary quantitative genetics (Lande 1979;Lande and Arnold 1983;Phillips and Arnold 1989) and adaptive dynamics (Dieckmann and Law 1996) are multidimensional, that is, describe simultaneous evolutionary dynamics of several traits. Moreover, these formulations can incorporate interactions between coevolving species, for instance in predator-prey systems (Dieckmann and Law 1996) or to study resource partitioning, character displacement, and ranges of coevolving species (Pacala and Roughgarden 1982;Taper and Case 1985;Case and Taper 2000). Recently, a few studies specifically looked at the effect of dimensionality in evolutionary and coevolutionary models and concluded that adding traits favors diversification (Doebeli and Ispolatov 2010, with an adaptive dynamics model) or the evolutionary escape of prey or hosts from an exploiter (Gilman et al. 2012, with a quantitative genetics approach). What is common to both these conclusions is that adding traits destabilizes evolutionary equilibria.
In this article, we present a general framework for analyzing the effect of the dimensionality of phenotype space on (co)evolutionary stability. Technically, this framework bridges the gap between quantitative genetics and adaptive dynamics-and thereby extends links already identified by Abrams et al. (1993) and Abrams (2001) to multiple dimensions. Using x * to denote combinations of traits at which fitness gradients vanish (i.e., x * is an evolutionary equilibrium), we investigate how dimensionality influences two types of stability, illustrated in figure 1: (i) convergence stability (Eshel and Motro 1981;Geritz et al. 1998), which refers to the question of whether x * is an attractor or a repellor for the evolutionary dynamics and relates to the direction of selection in the vicinity of x * (see fig. 1a, 1b), and (ii) evolutionary stability, which refers to the question of whether selection at x * is stabilizing or diversifying, that is, whether the trait variances tend to increase or decrease once the trait means are at x * (see fig. 1c, 1d). In the adaptive dynamics literature, evolutionary stability is defined as noninvadability (Geritz et al. 1998).   Figure 1: The different types of stability, illustrated when there are n p 2 traits under selection, focusing on one species. The initial population is plotted with gray points, and the postselection distribution is plotted with black points. The point x * is located at the intersection of the two axes. a and b correspond to the direction of selection near x * ; c and d correspond to whether selection is stabilizing at x * .
We will see that the conditions for noninvadability in an adaptive dynamics context match the conditions for selection to be stabilizing (as shown, e.g., by Sasaki and Dieckmann 2011;Débarre et al. 2013, with one trait). In addition, we will show that the impact of the number of traits under selection on (co)evolutionary stability can be dissected into three discrete effects that differently affect convergence stability and evolutionary stability. Finally, we illustrate the effects of dimensionality with specific examples of evolutionary dynamics due to frequencydependent competition and to victim-exploiter interactions.

Evolutionary Model
We consider a community of m coevolving species and assume that n i quantitative traits are under selection in species i. Each individual of species i is characterized by the values of its n i traits and is represented by a column vector , phenotype at the jth trait. We neglect features such as age structure or sex differences, that is, any within-species difference other than those due to differences in the traits under consideration. For simplicity, we also ignore the effects of environmental variation on traits and assume that we can identify phenotype and genotype. We are interested in the effect of selection on the distributions of traits in each species. In particular, we will describe how the means and variances of each trait in each species are affected by selection between two time steps. A time step typically corresponds to one generation (t p t ϩ 1) in models with discrete generations and to an infinitesimal time step (t p t ϩ dt) in continuous-time models. We will use flexible definitions that allow us to deal with the two kinds of models simultaneously.
We denote by f i the multivariate distribution of phenotypes in species i at time t (for notational convenience the time dependence is omitted). For each species i, the distribution f i is a function of n i variables, corresponding to the n i traits under selection in species i. We denote by w i (x i ) the relative fitness of individuals of species i with phenotype x i , given the current state of the community. This fitness measure is standardized such that the population average is equal to unity: w The fitness of an individual can depend directly on its own phenotype through adaptation to its abiotic environment as well as on the effects of interactions with other individuals; the strength of these interactions generally depends on the phenotypes of interacting individuals. Thus, in general, the fitness is a func- tion not only of x i , but also of the trait distributions f j for j p 1, ..., m, in the m coevolving species. We will derive explicit expressions for the fitness w i in "Illustrations" (eqq. [10], [13], [17]), but for the general arguments a specific expression of w i is not required. Given the various fitness functions w i , the change in the distribution of phenotypes due to selection, between two time steps is given by (1) To describe the evolutionary change in the phe-df (x )/dt i i notype distributions, we focus on summary variables such as means and variances. The mean of trait j in species i, , is given bȳ The covariance matrix G i contains terms G i,kl that are the covariances between traits k and l in species i:Ḡ We note that covariance matrices G i are real, symmetric (G i,kl p G i,lk ), and positive-semidefinite (all their eigenvalues are positive or null). We assume that the covariance matrices that we consider are positive-definite (all their eigenvalues are strictly positive); that is, we ignore the very particular case where some traits can be exactly expressed as linear combinations of the others (assuming that trait dimensionality has been reduced beforehand if this was the case; Lande and Arnold 1983). We further assume that this assumption remains correct even when covariance matrices evolve.
To obtain analytical results on the effect of selection, we assume that phenotypic variation in all m species is small, so that the difference between the trait value of an individual and the mean trait within its species is small, of the order d. This assumption, although limiting, is crucial to the derivation of our results (see "Discussion"). This allows us to Taylor-expand each fitness function w i around the mean traits in species i, so that a twice-x i differentiable fitness function can be approximated by a quadratic function: We note that in general, the arguments in both these functions are the mean trait in species i, , as well as thē x i distributions f j ( j p 1, ..., m). If these distributions are Gaussian with a fixed variance, they can be replaced by their means and covariances as arguments in the selection gradient and the Hessian. We also note that the selection gradient vector is commonly denoted by the Greek Dw (x ) i i letter b in quantitative genetics studies, while the Hessian matrix , also called matrix of quadratic and cor-2D w (x ) i i relational selection gradients (Blows et al. 2004), is denoted by g (Arnold 1992).
Because the Taylor expansion in equation (3a) is of second order, the expressions for the changes in each moment k of the distributions f i will also depend on the moments k ϩ 1 and k ϩ 2. Thus, changes in the means will depend on the covariance matrices and the third moments, and changes in the covariances will depend on third and fourth moments. Hence, our system of equations is not dynamically closed, as there are more variables than equations. We need to make moment closure assumptions, that is, assumptions about the shape of the trait distributions f i . As is classically done in quantitative genetics, we assume that the trait distributions f i are multivariate Gaussian. Gaussian distributions are symmetric: their third central moment is 0, and their fourth moments can be expressed as functions of covariances. Full expressions without the Gaussian approximation are presented in appendix A (apps. A-D available online), and simplified expressions are presented in the main text.

Changes in Mean Traits
Directional Change. To analyze directional evolution, we assume that the distributions f j , (j p 1, ... ,m), are Gaussian; hence, the fitness functions w i are functions of the trait x i , the mean trait values , and covariancēx , … , x 1 m matrices G 1 , ..., G m . Using the Taylor expansion (3a), the change in the vector of mean traits in species i is The derivation of equation (4) is presented in appendix A. This equation is the coevolutionary equivalent of traditional quantitative genetics equations (Lande 1979;Barton and Turelli 1987;Phillips and Arnold 1989;Abrams et al. 1993 Convergence Stability. The convergence stability of an equilibrium x * of the evolutionary dynamics (4) is given by the eigenvalues of the Jacobian matrix J, whose expression is given below in equations (5). We already assumed that phenotypic distributions are Gaussian; we now briefly assume that their covariance matrices are constant. This is a safe assumption because changes in the mean traits (eq. [4]) and changes in the covariance matrices (eq. [6] below) occur at different timescales. Hence, we can consider here that fitnesses w i are functions of the mean traits in all species (j p 1, ..., m) and we can treat x j covariance matrices G j as parameters-this assumption will be relaxed in the section "Changes in Trait Covariances." With this, the Jacobian matrix J is made up of blocks of size n i by n h , J ih : Each block J ih can be written as a product: where G i is the genetic covariance matrix of species i (defined in eq. [2b]), and C ih is a matrix that describes the effect of a change in the mean traits of species h on the selection gradient acting on species i: The equilibrium x * is convergent stable if all eigenvalues of the Jacobian J have negative real parts (see fig. 1a). Once x * is reached, mean traits in all species are at equilibrium, but this does not necessarily imply that the distributions of traits are frozen: higher moments can still change, and in particular, the covariance matrices G i may still evolve (Steppan et al. 2002). In the next section, we look at the effect of selection on trait (co)variances in each species.

Changes in Trait Covariances
Using the Taylor expansion (3a) around an equilibrium x * and the assumption that distributions are multivariate Gaussian, we find that the changes in trait covariances in species i due to selection can be written as The derivation of equation (6) is presented in appendix A. A similar expression, when only m p 1 species is evolving, can be found in, for example, Phillips and Arnold (1989, their eq. [2], and references therein; note that since we are at the equilibrium x * , the selection gradient b p Dw i is equal to 0). We say that x * is evolutionarily stable if selection reduces trait variances in all species (i.e., selection is stabilizing). In contrast, we say that the equilibrium x * is evolutionarily unstable if there is at least one (composite) direction in phenotype space along which the genetic variance increases (i.e., selection is diversifying, which means that a wider range of phenotypes will coexist). In terms of equation (6), this means that x * is evolutionarily stable if all eigenvalues of are negative; otherwise, x * is evo- The matrix D 2 w i appearing in equation (6) is the Hessian matrix of second derivatives of w i , evaluated at x * (see eq. [3b]); by definition, it is symmetric (D 2 w i,kl p D 2 w i,lk ). As noted previously, covariance matrices G i are also symmetric, and all their eigenvalues are real and positive. As a result, the signs of the eigenvalues of are 2 G 7 D w 7 G i i i the same as the signs of the eigenvalues of D 2 w i (Coppel 2009, pp. 238-239). Therefore, evolutionary stability of x * is completely determined by the Hessian matrices D 2 w i : x * is evolutionarily stable if and only if all eigenvalues of D 2 w i are negative (see fig. 1c). Similar conclusions are reached using an adaptive dynamics framework (Leimar 2009), a framework in which evolutionary stability refers to noninvadability.

Dimensionality and Stability
An evolutionary equilibrium x * is a point in the space of all traits of all species at which the fitness gradients of all species vanish. As explained above, two different kinds of stability need to be considered for such equilibria. On the one hand, x * is convergent stable if it is a (local) attractor for evolutionary dynamics (4), that is, if the real parts of all eigenvalues of the Jacobian (5b) are negative. Denoting by l(A) the eigenvalue with the largest real part of a square matrix A and by Re(z) the real part of a complex number z, this condition can be rewritten as On the other hand, x * is evolutionarily stable if trait variances decrease once the population means have reached the equilibrium x * , that is, if all eigenvalues of the Hessians D 2 w i are negative: 2 Gi, l(D w ) ! 0 (7b) i (note that Hessians are symmetric real matrices and that all eigenvalues of symmetric real matrices are real). The main goal of this article is to determine the effects of the dimensionality of phenotype space on these two types of stability. Thus, we are interested in determining the effect of dimensionality on the occurrence of eigenvalues with positive real part in the Jacobian matrix J (eq. [5]) and the Hessian matrices D 2 w i (eq. [3c]). Note: (Yes) means that there is no effect in symmetric competition models.
To tease apart different effects of dimensionality, we make the distinction of whether different traits are independent or not. Nonindependence can be due to epistatic interactions among traits or to trait correlations. Trait correlations correspond to off-diagonal elements in the covariance matrices G i ; in the absence of correlations in species i, G i is a diagonal matrix. Epistatic interactions occur when the effect of one phenotypic component on fitness depends on other phenotypic components. Formally, epistasis can be defined through the second derivative of the fitness function w i : for two traits k and l ( k, epistatic interactions occur if the second derivative 2 D w p i,kl is nonzero. In the absence of epistatic in- teractions, D 2 w i is a diagonal matrix. Keeping this definition in mind, we can identify three classes of effects of dimensionality on stability. These effects are described below and summarized in table 1.
Combinatorial Effect. The first effect of dimensionality occurs even when all traits are independent, that is, when there are no epistatic interactions and no correlations among traits; we call it "combinatorial effect." The combinatorial effect can most easily be seen in the context of evolutionary stability. When all traits are independent, the Hessian matrices D 2 w i are diagonal: their eigenvalues are their diagonal elements. Stability conditions can therefore be formulated separately in each phenotypic dimension. The more traits there are, the more likely it becomes that at least one diagonal element in D 2 w i will be positive, which destabilizes the entire system of multidimensional phenotypes in the direction of the corresponding trait. Thus, evolutionary instability becomes more likely as the dimension of phenotype space increases. A similar effect occurs for convergence stability. Even when all traits are independent, that is, in the absence of epistatic interactions, the Jacobian matrix J determining convergence stability is in general not diagonal, and destabilization, when it occurs, is generally along composite phenotypic directions instead of single traits. However, the same reasoning still applies. Traits being independent, the more dimensions there are, the more eigenvalues of J there will be, and the higher the chance that at least one of them is positive.
Epistasis. The second effect of dimensionality is due to epistatic interactions. Such interactions lead to correlational selection (Schluter and Nychka 1994;Brodie et al. 1995;Blows and Brooks 2003) and correspond to Arnold's (1992) third horizon of selective constraints. Assuming that a high-dimensional system without epistasis is stable (either convergent stable or evolutionarily stable), we can ask how introducing epistatic interactions into such a system affects stability. This effect can be most easily seen by considering evolutionary stability, that is, the Hessian matrices D 2 w i . In these matrices, epistatic interactions correspond to nonzero off-diagonal elements. Let us consider a "full" Hessian matrix D 2 w i and a corresponding diagonal matrix , whose diagonal elements are the same as 2 (0) D w i the diagonal elements of D 2 w i but whose off-diagonal elements are 0. A linear algebra theorem (the eigenvalue interlacing theorem; Cauchy 1891) tells us that the largest eigenvalue of D 2 w i (which is a symmetric matrix) is greater than (or equal to) its largest diagonal element: i i This means that starting from a system that is evolutionarily stable in the absence of epistasis (i.e., , adding epistasis can lead to a situation where . This mechanism is gen- eral: epistatic interactions increase the likelihood of evolutionary instability. We will illustrate this in our examples, and we will see that the magnitude of the destabilizing effect of epistasis increases with the dimension of phenotype space. Convergence stability is also affected by epistasis, since epistatic interactions also affect the Jacobian matrix. In general, the J matrix is too complex for any general statement to be made. However, as we will see in "Illustrations," in some particular and classical models it can be shown that epistatic interactions tend to render an equilibrium x * convergence unstable (i.e., turn x * into a repeller) and that this effect is again larger in higher dimension.
Trait Correlations. The third effect of dimensionality is due to correlations among traits, that is, to nonzero offdiagonal elements in the covariance matrices G i (see fig.  2c). Such correlations can, for instance, be due to pleiotropy and correspond to Arnold's (1992) second horizon of genetic constraints. As mentioned above, covariance matrices affect only the Jacobian matrix J at an equilibrium x * and not the Hessian matrices. Hence, covariance matrices can affect only convergence stability of x * but not evolutionary stability.
In some situations, correlations do not have an effect on convergence stability. For example, if a single species is evolving, the Jacobian matrix is of the form J p (Co)Evolutionary Stability 163  , where G is the covariance matrix and C is the G 7 C Jacobian of the selection gradient. If the matrix C is symmetric, which as is the case for symmetric competition models (see "Illustrations" and Matessi and Schneider 2009;Doebeli and Ispolatov 2013), then multiplying it by G has no effect on the sign of the eigenvalues, because G is positive definite. Therefore, genetic covariances have no effect on convergence stability if C is symmetric. In particular, if a symmetric C matrix has only negative eigenvalues, then the equilibrium x * is convergent stable regardless of correlations between phenotypic components. This situation was termed "strong convergence stability" by Leimar (2009).
When C is not symmetric, however, or if there are more than one coevolving species, trait correlations can affect convergence stability. For a single species and an asymmetric C, this effect is actually weak (see fig. D1, available online). When multiple species coevolve, it is not possible to draw general conclusions about the magnitude and direction of this effect, except for some special cases where C is simple enough (as in Gilman et al. 2012). In the example of coevolution in victim-exploiter systems given in the next section, genetic covariances tend to destabilize convergent stable equilibria.

Illustrations
We illustrate the effects of the dimensionality of phenotype space on convergence stability and evolutionary stability using three examples: a model of competition for resources, a model of adaptation to two different habitats under soft selection (the Levene model; Levene 1953;Christiansen 1975;Kisdi 2001;Débarre and Gandon 2011), and a model of victim-exploiter coevolution. These models correspond to classical models in evolutionary ecology that were either published before with only one dimension or published in less general forms (i.e., with specific functions, often Gaussian, or with more restrictive assumptions). The first two examples can be grouped into a broader class of models, symmetric competition models, which we now describe.

Symmetric Competition Models
Consider a single species (m p 1) in which competition for resources determines the evolutionary dynamics of an n-dimensional phenotype (for notational convenience we drop the species subscript i in the following). The ndimensional phenotype is assumed to determine both the strength of competition and resource availability. The strength of competition between individuals with phenotypes x p {x 1 , ..., x n } T and y p {y 1 , ..., y n } T is described by a function a(x, y), called the competition kernel. For symmetric competition, it is assumed that the strength of competition increases with phenotypic similarity of competing individuals, and hence that the competition kernel a(x, y) is maximal along the diagonal x p y (we also assume a(x, x) p 1 for all x). The amount of resources available to individuals of phenotype x is described by the carrying capacity function K(x).
If N(x, t) denotes the density distribution of the phenotype x at time t, then for a given competition kernel a(x, y) and a given carrying capacity function K(x), the dynamics of this density distribution is given by the following equation: proportion of simulations where x * is stable when there are epistatic interactions given that x * is stable in their absence. Trait correlations do not influence stability in this class of models and are therefore not shown. The x * equilibrium is always convergent stable in the model with a unimodal distribution of resources K, as a direct consequence of the shape of K, so the gray dots are not displayed in a and b. There are 10 6 replicates for each dot of each figure.
We know from equation (8) that adding diagonal elements to the Hessian matrix, that is, adding epistatic interactions among traits while keeping the same diagonal elements (a jj and k jj terms), increases the range of eigenvalues. This means that x * can be evolutionarily stable if the traits are independent but may become evolutionarily unstable if epistatic interactions are added. To numerically investigate how frequently epistasis destabilizes x * , we draw random negative definite HK and Ha matrices ("full" matrices, with off-diagonal elements corresponding to epistatic interactions), and derive from them equivalent diagonal matrices HK (0) and Ha (0) , which have the same diagonal elements and correspond to the cases where all traits are independent (no epistasis); and as previously, we set K(0) p 1. Figure 3b represents the proportion of simulations where x * is stable with epistatic interactions given that it is stable in the absence of epistatic interactions. The decrease in this proportion as the number of phenotypic dimensions increases represents the destabilizing effect of epistasis, controlling for the combinatorial effect. This gen-eralizes the results of Doebeli and Ispolatov (2010) to any well-behaved (twice-differentiable) competition and carrying capacities functions, and shows the links between their adaptive dynamics approach and the moment-based quantitative genetics approach used here (Abrams et al. 1993;Abrams 2001).
Levene Models. In its simplest form, the Levene model (Levene 1953) describes the evolution of a population of individuals divided into two demes in proportions c and (1 Ϫ c) (with 0 ! c ! 1) with different selective conditions and under soft selection (Christiansen 1975). Generations are discrete and nonoverlapping, and all the individuals are pooled and redistributed at the end of each generation. There are no explicit density dynamics, but selection is frequency dependent. We again assume a single evolving species (m p 1) and drop the species index for convenience. The n-dimensional phenotype x determines fecundity in the two habitats, denoted by q h (x), h p 1, 2. With these assumptions, the postselection distribution of phenotypes is where the term in brackets, w(x), is the relative fitness of phenotype x, and where is the mean fecundity in habitat h (h p 1, 2). It was shown previously that Roughgarden's competition model (Roughgarden 1972(Roughgarden , 1979Doebeli and Dieckmann 2000) and the Levene model (Levene 1953) can be seen as two ends of a continuum (Débarre 2012). Here we go one step further and show that the Levene model is in fact equivalent to a symmetric competition model of the form (9), with suitable competition kernel a and carrying capacity K. The derivation of K and a is explained in more detail in appendix B. Here we note that the corresponding carrying capacity function is c 1Ϫc 1 2 We refer to equation (B8b) for the expression of the corresponding competition kernel. The multidimensional Levene model is hence a symmetric competition model. It follows that the covariance matrix G has no effect on convergence stability and that the convergent stable equilibrium points x * are exactly the local maxima of the function K(x) given by equation (14). Note, however, that contrary to the previous example with unimodal carrying capacity functions, an equilibrium x * is not necessarily a local maximum of K (it could be a saddle point or a local minimum) and hence is not necessarily convergent stable. Whether derived directly using w(x) given in equation (13), or via the equivalent K and a functions and the expressions (11) and (12), we find that the matrix C controlling convergence stability is where is the gradient of the fecundity function q 2 in * Dq 2 habitat 2, evaluated at x * , and where are the corre-* Hq h sponding Hessian matrices, also evaluated x * . Similarly, the Hessian matrix D 2 w determining evolutionary stability is To illustrate the effects of dimensionality on convergence and evolutionary stability in the Levene model, we make the further symmetry assumptions that both habitats are in equal frequencies (c p 1/2) and that the fecundity functions are symmetric (q 2 (x) p q 1 (Ϫx)); this ensures that the point x * p 0 is an equilibrium. Epistatic interactions are then present if and only if Ѩ 2 q 1 /Ѩx k Ѩx l ( 0 for at least one pair k ( l. The consequences of the combinatorial effect and of epistatic interactions on the destabilization of x * , estimated from our simulations, are illustrated in figure 3c and 3d, respectively. Random diagonal Hq 1 and Hq 2 matrices are drawn to assess the combinatorial effect, and figure 3c represents the proportion of simulations where x * is evolutionarily stable (black dots) or convergent stable (gray dots). As previously, we then control for the combinatorial effect when assessing the magnitude of the effect of epistasis on stability; figure  3d represents the proportion of simulations where x * is evolutionary (black) or convergent (gray) stable in the presence of epistatic interactions, given that it is stable in their absence. We see that both the combinatorial effect and epistasis tend to generate both convergence instability and evolutionary instability, and the likelihood of destabilization increases with increasing dimension of phenotype space.
We note that x * is a generalist strategy; when x * is not evolutionarily stable, the population diversifies and specializes to the local habitats. Hence, because increasing the number of dimensions makes selection more diversifying at x * , it also increases the potential for local adaptation.

Victim-Exploiter Coevolution
In our last example, we consider the coevolution of two different species, a victim and an exploiter, so that the community contains m p 2 species, labelled 1 (victim) and 2 (exploiter). In each species i, n i traits are either under stabilizing selection and/or involved in interspecific interactions. We use a classical Lotka-Volterra formulation of the model, as presented, for instance, in Doebeli (2011). The variables n 1 (x 1 ) and n 2 (x 2 ) represent the densities of individuals of species 1 and 2 with traits x 1 and x 2 , at time t. We assume that species 1 grows logistically in the absence of species 2, and that the carrying capacity, K(x 1 ), depends on the phenotype x 1 (first term in eq. [16a]). Interactions between the two species affect the fitness of individuals of both species, with functions b i describing the effects of interactions (second term in eq. [16a] for the effect on the victims and first term in eq. [16b] for the effect on the exploiters). In the absence of interactions, the individuals of species 2 (exploiter) die at a constant rate d (second term in eq. [16a]). This model is expressed in continuous time, and interactions are density dependent: ͵ 1 1 2 2 2 2 1 1 ( ) Ѩf (x ) f (y ) 1 Contrary to our previous examples, this C matrix is in general not symmetric even when there are no epistasic interactions among traits in the K and b i functions. As a result, the signs of the eigenvalues of the Jacobian J given by equation (5b) depend both on C and on the covariance matrices G 1 and G 2 . In particular, correlations among traits (i.e., nonzero off-diagonal elements in the G i matrices) will play a role in whether x * is convergent stable. To investigate the magnitude of this effect numerically, we randomly draw G i matrices (positive definite matrices whose elements are drawn from uniform distributions), and compare the stability of J with full G matrices to the stability of J with equivalent diagonal G (0) matrices-all the off-diagonal elements being set equal to 0; then, we assess the effect of trait covariances while controlling for the combinatorial effect. Figure 4 confirms that the three effects of dimensionality that we identified are at work (see gray dots). All effects (combinatorial effect, epistasis, trait correlations) affect convergence stability and increase with the number of dimensions. In particular, and as already noted by Gilman et al. (2012) with a more limited model, trait correlations affect convergence stability ( fig. 4c): an equilibrium that is an attractor (CS) in the absence of trait correlations may become a repellor (non-CS) with trait correlations, leading to the evolutionary escape of the species that benefits from being different; it is the victim (species 1) with this particular interaction function, although, as emphasized previously, this could as well be the exploiter (species 2) with another type of interaction function.
We now consider evolutionary stability. The Hessian matrices of fitness for each species are By assumption in this example, is negative definite * Hb 2 (this reflects the fact that exploiters benefit from matching the victims), and hence so is D 2 w 2 : selection on the exploiters' traits is stabilizing when the victims and exploiters match.
The situation for evolutionary stability in the victim as described by D 2 w 1 is similar to the situation in one-species competition models. In the absence of epistatic interactions between the different traits, D 2 w 1 is a diagonal matrix, and selection is diversifying if at least one of the diagonal elements is positive (combinatorial effect; see black points in fig. 4a). Even if all these diagonal elements are negative, epistatic interactions between traits, obtained by adding  nondiagonal elements to D 2 w 1 , expand the range of the eigenvalues and increase the chance that at least one of the eigenvalues is positive ( fig. 4b). Again, the strength of this effect increases with the dimensionality of phenotype space, so that evolutionary instability, and hence diversification, becomes more likely in higher dimensions. Such diversification in species 1 may in turn favor diversification in species 2 (Doebeli and Dieckmann 2000;Doebeli 2011). We note that Gilman et al. (2012) focused on convergence stability and did not investigate the effect of the dimensionality of phenotype space on evolutionary stability; they also did not consider the effect of epistatic interactions.

Three Effects of Dimensionality
In this article, we used a moment-based approach to investigate the effect of dimensionality of phenotype space on the coevolutionary stability of a community of m species when selection acts on quantitative traits. Each species i is characterized by a multivariate distribution f i of traits, and we focus on the effect of frequency-dependent selection on the first two moments of these distributions, namely, their means and covariance matrices G i . Co-x i evolutionary equilibria x * are points in the space of all traits of all species at which fitness gradients of all species vanish. Two types of stability must be considered. On the one hand, an equilibrium x * is convergent stable (Eshel and Motro 1981) if it is a local attractor for the (co)evolutionary dynamics. On the other hand, an equilibrium x * is evolutionarily stable (Geritz et al. 1998) if selection is stabilizing, and hence phenotypic variances are decreasing, once the population means are at the equilibrium point. Figure 1 illustrates these two types of equilibria. Our analysis shows that, generally, increased dimensionality of phenotype space destabilizes equilibria in both senses: x * is more likely to be a repeller for the evolutionary dynamics, and selection at x * is more likely to be disruptive as the number of phenotypic dimensions increases. We specifically pinpointed which features of dimensionality are responsible for this and summarize these effects in table 1.
Dimensionality affects stability even when all traits are independent (i.e., when there are no epistatic interactions among them and when they are uncorrelated). We call this first effect of dimensionality a "combinatorial effect": increasing the number of traits increases the chance of instability in at least one phenotypic direction.
Epistatic interactions correspond to nonadditive interactions between different traits and result in the addition of off-diagonal terms in the Hessian matrices of the fitness functions (i.e., there are epistatic interactions between two different traits k and l in species i when Ѩ 2 w/Ѩx i,k Ѩx i,l ( 0, where w is the fitness function). Essentially, adding epistasis terms offers new (composite) directions in phenotype space along which equilibria can be unstable and hence increases the likelihood of instability.
The effect of correlations among traits is more limited. First, trait correlations never affect evolutionary stability. Second, their effect on convergence stability is limited. When the Jacobian of the selection gradient (i.e., the matrix C defined in eq. [5c]) is symmetric, which, for example, is the defining feature of symmetric competition models (Doebeli and Ispolatov 2013), trait correlations have no effect on convergence stability. Convergence stability with symmetric C matrices has therefore been called "strong convergence stability" by Leimar (2009). When only one species is evolving, and when C is not symmetric, the destabilization due to trait correlations is rather limited (see fig. D1). It is only when there are multiple species coevolving, as in the example of victim-exploiter coevo-(Co)Evolutionary Stability 169 lution presented in the previous section, that trait correlations can turn equilibria that are attractors in the absence of correlations into repellers.

The Incomplete Tale of Vectors and Matrices
Technically speaking, it is interesting to note that our moment-based approach yields equations that are similar to those used in evolutionary quantitative genetics (Lande 1979;Lande and Arnold 1983;Phillips and Arnold 1989, for m p 1 species) and that the conditions we derived for convergence stability and evolutionary stability are the same as the ones obtained in adaptive dynamics models (Dieckmann and Law 1996;Geritz et al. 1998;Leimar 2009). Indeed, we focused on the effect of selection on trait distributions, an approach that is the same for clonal or sexual organisms. These links between the quantitative genetics and adaptive dynamics frameworks had already been discussed in one-dimensional, single-species cases (see, e.g., Abrams et al. 1993;Abrams 2001;Sasaki and Dieckmann 2011;Débarre et al. 2013). Our results show how to extend them to evolution of multiple species in multiple phenotypic dimensions.
Two critical assumptions allowed us to derive our results. The first is the assumption of small phenotypic variation in each species, which allowed us to write Taylor expansions of the fitness functions and to neglect high orders (see eq. [3a]). We could therefore express individual fitnesses as a function of fitness of mean traits in each species. But even with this assumption, the expressions for the changes of particular moments of the phenotype distributions depend in general on the two next-higher moments, leading to dynamic insufficiency. A momentclosure approximation-here the assumption that all distributions are multivariate Gaussians-is needed to circumvent this problem. In a model with m coevolving species, we thus arrive at a system of 2m equations: m equations for the change in mean traits (see eq. [4]) and m equations for the change in covariance matrices (see eq. [6]).
It is important to keep in mind that these expressions rely on our moment-closure approximation and that higher moments of the trait distributions would matter should the assumption of normality not hold, for example, when spatially heterogeneous selection generates asymmetric (skewed) trait distributions (Yeaman and Guillaume 2009;Débarre et al. 2013). Third moments, corresponding to asymmetries in trait distributions, do for instance appear in the general expression for the change in mean traits (see eq. [A8]; Rice 2004b; Brodie and McGlothlin 2007). Information about third moments cannot be represented in matrices but requires three-dimensional arrays. Similarly, kurtoses (represented in four-dimensional arrays) appear in the general expression for the change in trait variances (see eq. [A9d]). Hence, selection gradient vectors, covariance matrices, and matrices of correlational selection, while extremely helpful, are in general not sufficient to describe the effect of selection on trait distributions (Rice 2004a). It remains to be seen whether models including higher-order effects would corroborate our results. Building on Blows and Walsh's (2009) analogy ("spherical cows grazing in flatland"), we can say that not only are cows not spherical, but they are not just ellipsoidal either: both their asymmetry (third moments) and tails (fourth moments) potentially matter too.

Developmental Evolution: The Origin of Epistatic Interactions and Trait Correlations
Our model takes ecological interactions within and between species into account. In particular, selection can be both density and frequency dependent. We did not, however, specify how epistatic interactions or covariances among traits could emerge-what corresponds to Arnold's (1992) fourth horizon of developmental and functional constraints. Rice (2002Rice ( , 2004b integrated the effect of complex developmental interactions in a population genetics model. Extending this framework to coevolutionary interactions may be particularly challenging, but it constitutes an interesting avenue for future research.

Multidimensional Santa Rosalia
Our results show that increased dimensionality makes selection more disruptive. A key motivation for our work, as argued in the introduction, was that selection in nature is likely to affect more than just one trait at a time. Yet, our study is limited to the effect of selection on trait distributions, and it does not include the effect of sexual reproduction, nor does it allow us to investigate the evolution of reproductive isolation or the effect of genetic constraints on speciation (Felsenstein 1981). In other words, while we can conclude that increased dimensionality makes it more likely that selection will be disruptive in at least one (composite) direction, it remains to be investigated whether this would result in speciation. Nevertheless, the consideration of trait dimensionality may contribute to explaining the maintenance of phenotypic variation in nature.