Multidimensional divergent selection, local adaptation, and speciation

Divergent selection applied to one or more traits drives local adaptation and may lead to ecological speciation. Divergent selection on many traits might be termed “multidimensional” divergent selection. There is a commonly held view that multidimensional divergent selection is likely to promote local adaptation and speciation to a greater extent than unidimensional divergent selection. We disentangle the core concepts underlying dimensionality as a property of the environment, phenotypes, and genome. In particular, we identify a need to separate the overall strength of selection and the number of loci affected from dimensionality per se, and to distinguish divergence dimensionality from dimensionality of stabilizing selection. We then critically scrutinize this commonly held view that multidimensional selection promotes speciation, re‐examining the evidence base from theory, experiments, and nature. We conclude that the evidence base is currently weak and generally suffers from confounding of possible causal effects. Finally, we propose several mechanisms by which multidimensional divergent selection and related processes might influence divergence, both as a driver and as a barrier.

Populations adapt in response to natural selection to optimize fitness in their native environment. Any number of environmental pressures might drive adaptation, either through stabilizing, directional, or divergent selection. Interactions between the environment and the traits of an organism determine fitness; typically, many selection pressures and traits contribute to fitness, some operating independently, whereas others interact (Débarre et al. 2014). Divergent selection generates local adaptation, potentially leading to speciation, by driving the divergence of populations toward distinct fitness optima in a heterogeneous environment (Kawecki and Ebert 2004;Rundle and Nosil 2005;Nosil 2012). For any pair of locally adapted populations, many shared environmental variables are likely to generate stabilizing selection in both habitats (Gavrilets 1997;Langerhans and Riesch 2013), but there is one axis in niche space that separates them (Yamaguchi and Otto 2020). This axis of separation between habitats might be dominated by a single environmental variable, or might impose selection on a single trait, in which case divergent selection could be described as "unidimensional." Alternatively, one might identify differentiation in multiple environmental variables or traits, generating "multidimensional" divergent selection (Rice and Hostert 1993; sometimes called "multifarious selection," e.g., Feder and Nosil 2010).
There is a broad consensus that the higher the dimensionality of divergent selection, the more likely it is to generate local adaptation and speciation. This stems from classic reviews (Nosil et al. 2009b; Rice and Hostert 1993;Nosil and Harmon 2009) and is often repeated (Smadja and Butlin 2011;Butlin et al. 2012;Langerhans and Riesch 2013;Ravinet et al. 2017). However, nearly three decades on from the original proposal, there remains little clear evidence supporting this hypothesis. Furthermore, the multiple mechanisms by which increased dimensionality might influence local adaptation and speciation have not been fully distinguished, either theoretically or in empirical studies. Here, we argue that the proposed effects of increased dimensionality of divergent selection might instead be attributed to increased overall strength of divergent selection, increased dimensionality of stabilizing selection, increased number of loci under selection, or to other possible correlates of dimensionality. We highlight the need for theoretical and empirical work to test the impact of dimensionality on local adaptation and speciation in ways that help to understand the mechanisms of action.

Defining Dimensionality
Although we are mainly concerned about response to divergent selection, it helps to begin by considering the dimensionality of selection in a single habitat. This requires an understanding of the mapping between environmental, phenotypic, and genetic variation.

OVERALL DIMENSIONALITY
First, there is the dimensionality of the environment. The habitat occupied by a population can be described by measuring many environmental variables. For each environmental variable, there is a range within which the population can survive and reproduce and this defines a hypervolume that describes the population's niche (Hutchinson 1957). Because the environmental variables are likely to be correlated, the effective dimensionality of this volume is lower than the number of measurable variables. The leading eigenvectors of the matrix of environmental variables define a set of orthogonal environmental axes and the complexity or dimensionality of the environment can be described by the number of these axes required (de Paula-Souza and Diniz-Filho 2020).
Phenotypic variation can be described in a similar way. A large number of possible phenotypic traits can be measured but correlations among traits mean that the dimensionality of phenotypic variation is likely to be much lower (Kirkpatrick and Meyer 2004;McGuigan et al. 2005). Therefore, phenotypic variation can be described by a smaller number of orthogonal phenotypic axes. The dimensionality of phenotypic variation depends in part on the underlying pattern of genetic variation. This can be described in quantitative genetic terms by the genetic variancecovariance matrix (G-matrix; Guillaume 2011), which also has a limited number of orthogonal axes that may differ from the major axes of the phenotypic matrix. Alternatively, genetic variation can be described at the level of individual genetic variants and their patterns of linkage disequilibrium. The forms of the genetic and phenotypic variation will depend on the history of selection on the population, and so on the environment (Yamaguchi and Otto 2020).
The dimensionality of selection depends on the interaction between environmental and genetic variation (Kirkpatrick and Meyer 2004). This is implicit in Hutchinson's niche definition because the environmental variables that matter are those that impose limits on the region that the population can occupy. There are environmental variables that have little or no effect on these limits and, similarly, there are genes and phenotypes whose variability does not influence fitness within the given habitat. Tenaillon (2014, p. 194) defines complexity as "a quantitative measure that reflects the number of variationally quasi-independent traits an organism is exposing to the action of natural selection in a given environment." This definition arises from consideration of Fisher's Geometric Model of selection toward an environmental optimum (Fisher 1930). Tenaillon emphasizes that its value is labile, as environments, phenotypes, and genotypes evolve, and depends on the time scale under consideration, typically being lower for shorter durations. It can be thought of as the number of orthogonal phenotypic or genetic axes on which there is effective stabilizing selection. The number of environmental axes that impose appreciable stabilizing selection may be similar but the mapping between the two sets of axes may not be simple.

DIVERGENCE DIMENSIONALITY
To consider local adaptation, it is necessary to extend this thinking to two habitats, each of which can be represented by Fisher's Geometric Model with a single optimum that differs between habitats. Therefore, there is stabilizing selection around each optimum, and divergent selection between habitats. Making the simplifying assumption that the same set of environmental or phenotypic axes underlies selection in both environments, these optima are two points in the same multidimensional space. Clearly they can be connected by a single axis, which we will refer to as the axis of divergent selection. In this sense, divergent selection is always unidimensional (unless there are more than two habitats under consideration). However, the axis of divergent selection might be aligned with a single axis in the environmental or phenotypic space, or even with a single underlying environmental variable or phenotypic trait (Fig. 1). This might be considered unidimensional divergent selection in contrast to cases where the divergent selection axis implies selection on multiple phenotypic traits or axes in response to multiple environmental variables or axes. Unfortunately, the literature on the role of multidimensional selection in local adaptation and speciation rarely makes these distinctions (see below). MacPherson et al. (2015) provide an example to illustrate their model of local adaptation in a metapopulation. Femur length and head width are phenotypes in the cricket, Gryllus firmus (Bégin and Roff 2001), that are genetically correlated (G-matrix, covariance is positive). MacPherson et al. also envisage some environmental selection on head width and on femur length. We might imagine, for the sake of argument, that selection is due to two environmental variables: a resource variable and a predation risk variable that impose selection on the two phenotypes, respectively. Rather than varying independently across demes, there is a positive correlation between these environmental variables is an optimum (dark blue dot). These optima are separated along one environmental axis (x 2 ) with stabilizing selection along a different axis (x 1 ) creating unidimensional divergent selection between multidimensional niches. We might envisage that an ancestral population evolved from the origin to fill these niches via the trajectories shown by the red dashed arrows. Transitioning to the middle-left plot (B), divergent selection is now applied over two environmental axes (x 1 and x 2 ). There is still a single axis in multidimensional environment space separating the optima, but as two axes are now divergent, there is the potential for stronger divergent selection as optima become more distantly separated. Multidimensional divergent selection may have correlated responses. The example in the bottom-left plot (C) shows two orthogonal traits have now diverged: t 1 in response to selection on x 1 , and t 2 in response to selection on x 2 . By chance, t 1 is a multiple-effect trait and greater divergence in t 1 produces assortative mating. This additional barrier would not have arisen without multidimensional divergent selection. An additional barrier mechanism is shown in the bottom-right plot (D). Multidimensional divergent selection has produced regions of divergence, as measured by FST, around multiple locally adaptive loci (LA loci). Linkage between locally adaptive loci and DMI loci produces correlated divergence in the DMI loci, generating an additional barrier to gene flow. Furthermore, returning to the initial two-dimensional niche representation, overall dimensionality might increase from two-dimensional (A) to threedimensional (D) via the introduction of a third environmental axis (x 3 ). This additional axis might provide either stabilizing or directional selection, but regardless will increase transgressive incompatibilities as hybrid fitness deviates from the optimum along the additional axis.
EVOLUTION SEPTEMBER 2021 2 1 6 9 (described by an E-matrix). The strengths of the environmental or genetic correlations influence the effective dimensionality: if they are strong, dimensionality is close to one; if weak, it is close to two. Because the genetic and environmental correlations need not be the same, dimensionality depends on the viewpoint. The dimensionality of selection depends on the interaction between genetic and environmental variation: for example, if head width and femur length were perfectly genetically correlated, the dimensionality of selection would be one, regardless of the correlation between environmental variables. MacPherson et al. (2015) show that local adaptation increases with dimensionality in their model (see below) but the increase is greater if the G and Ematrices are aligned, illustrating this interaction. This example can be extended to show the distinction between divergence dimensionality and the overall dimensionality of the selective environment. Suppose that the crickets also vary in color and that matching to the background color of the environment influences fitness, but that this background color does not vary among patches. The dimensionalities of the environment and of the phenotype are increased but the divergence dimensionality is not. Finally, the crickets might vary in bristle pattern, which has no influence on fitness in this environment, and the habitats might vary in sward height, over a range that has no impact on cricket survival and reproduction. These variables would increase phenotypic or environmental dimensionality but not the dimensionality of selection. A "roadmap" of the perspectives of overall versus divergence dimensionality and their correlated effects is provided in Figure 1.

Multidimensionality and Extrinsic Isolation
Why should multidimensional selection increase local adaptation and the chances of progress toward speciation, compared with unidimensional selection? If a full picture of the effects of dimensionality per se is to be achieved, there are several potentially covarying effects that must be addressed. Here, we distinguish three of these possible effects: (1) The intensity of divergent selection increases with divergence dimensionality; (2) Other components of selection increase with overall dimensionality; and (3) The genetic dimensionality increases with divergence dimensionality.

STRONGER DIVERGENT SELECTION
Considering divergent selection between two distinct habitats, overall selection for local adaptation can be conceptualized as the distance between environmental optima in multidimensional space (defined by orthogonal phenotypic axes). For simplicity, the axes can be scaled so that fitness in each habitat follows a Gaussian decline, equally in all directions from the habitat op-timum, and we can assume equal fitness for well-adapted phenotypes in the two habitats. In reality, this scaling might not be possible across two different habitats but this complexity does not influence our arguments here. The fitness of an individual depends on its phenotype and the habitat in which it is found. Here, we will focus on the fitness of individuals that are well adapted to one habitat when they are either in their home environment or the alternative environment, and on the fitness of hybrids in the habitat in which they have higher fitness (cf. Thompson et al. 2019, for example). The Euclidean distance between environmental optima then determines both the fitness of a migrant between habitats and the reduction in fitness of a hybrid whose phenotype is at the mid-point between optima (Fig. 2).
In this framework, there is no necessary relationship between the intensity of divergent selection and either the dimensionality of divergence or the dimensionality of the trait space. With increasing dimensionality of divergent selection, that is, where the two optima differ on a greater number of orthogonal phenotypic axes under selection, there are two possible extreme modeling assumptions. The first is that overall selection increases with the number of selection pressures, as each selection pressure contributes an additional fitness reduction for migrants or hybrids. Alternatively, overall selection might remain constant but be spread over more axes, implying weaker selection per-axis (Nosil et al. 2009b). In the first scenario, selection is "additive" across dimensions, whereas in the second it is "diluted." We will use these terms to refer to the two modes of multidimensional selection while recognizing that there is a continuum of intermediate possibilities (Fig. 2). One might expect to find cases throughout this continuum in nature and the empirical challenge is to measure both divergence dimensionality and the overall intensity of selection if their effects are to be separated. Attempting to predict a priori how stressors will interact is a significant challenge, mired by ecological and temporal complexity, although significant strides are being made (Galic et al. 2018;Birk et al. 2020;Orr et al. 2020).

INCREASED OVERALL DIMENSIONALITY
Hybrids between populations adapted to two distinct habitats have reduced fitness because of their intermediate phenotypes, which fall between two adaptive optima. However, they also have reduced fitness due to segregation of alleles that influence other phenotypic axes, that is, because their phenotypes fall away from the line directly connecting the two optima resulting in further fitness reduction (Fig. 2). This fitness reduction is known as "transgressive incompatibility" (Chevin et al. 2014). It increases as overall dimensionality, not divergence dimensionality, increases (Chevin et al. 2014;Thompson et al. 2019). It is also more dependent on the history of the populations than on their current separation in phenotypic space, even being experienced by populations Concepts in adaptation to a multidimensional landscape. Panels depict overlays of two-dimensional adaptive landscapes in two different environments (separated by thick black lines). Environments have distinct fitness optima (intersection of dashed black lines), which may vary along one or both phenotypic axes. Selection over the two environments may therefore be either stabilizing (single optimum) or divergent (two optima) for each phenotypic axis. Contours represent maximum fitness of a given phenotypic combination over both environments. In all cases, we assume that phenotypes are scaled such that fitness surface is Gaussian with equal variances in both phenotypes and no covariance between phenotypes. The adaptive landscapes of panels A and B are identical, and only phenotype 1 is under divergent selection with two adaptive optima; phenotype 2 is under stabilizing selection. Comparison of this landscape with panels C and D shows the two modes of multidimensional selection. Here, both phenotypes are under divergent selection. Under additive multidimensionality (comparison of A/B with C), the per-axis divergence is held constant, producing a greater Euclidean distance

between fitness optima and deeper fitness valley (contours), lowering the fitness of migrants and hybrids, respectively. Under diluted multidimensionality (comparison of A/B with D), overall selection remains constant: the distance between multidimensional optima is equal to the unidimensional case. The divergence of individual phenotypes (distance between dashed lines) is smaller than under unidimensional divergent selection (A and B) or additive two-dimensional selection (C). In all panels, adaptation proceeds from an ancestral point to adaptive optima via mutations (yellow arrows). Each mutation has pleiotropic effects on both phenotypic axes. Local adaptation
can proceed via few large-effect mutations (C) or via many small-effect mutations (D), irrespective of mode or dimensionality. With more mutations, there is an increased probability of producing a constitutive incompatibility. Theory states that at higher overall environmental dimensionality, more mutations are required for adaptation to a local optimum (Chevin et al. 2014). Panels A and B also depict the effect of transgressive incompatibilities. In panel B, the mutational trajectory is less closely aligned with the axis between optima in phenotype space than in panel A. These off-axis mutational effects produce "segregation variance" in hybrid offspring, as phenotypes (red dots) vary widely along axes that are orthogonal to the discriminant axis separating optima (e.g., axis 2 in panels A and B; Chevin et al. 2014; Thompson et al. 2019). Off-axis variance in hybrid phenotypes is predicted to decrease with the alignment seen in panel A and to increase with the overall environmental dimensionality.
in identical environments. This is because the set of loci at which two populations differs depends on the evolutionary trajectory by which they have reached their current state. Transgressive incompatibilities can increase the barrier to gene flow between populations and this might increase local adaptation as well as making speciation more likely. Therefore, although the effect is not dependent on divergence dimensionality, it can produce a positive relationship between overall dimensionality and local adaptation. It is possible to distinguish this fitness cost experimentally, for example, by measuring the fitness of hybrids between populations independently adapted to similar environments as well as those adapted to distinct environments (e.g., Johansen-Morris and Latta 2006; Van Der Sluijs et al. 2008), or how hybrid fitness varies by phenotype and environment (Arnegard et al. 2014).

INCREASED GENETIC DIMENSIONALITY
There is potentially a correlation between the divergence dimensionality and the number of loci under divergent selection. In turn, the number of loci, or other aspects of genetic dimensionality determined by covariance among loci, might influence the potential for local adaptation and speciation. The number of loci under selection is a measure of genomic dimensionality, but the extent to which each locus represents an independent dimension is modulated by pleiotropy, because an allele can influence multiple traits, and by genetic architecture, because nearby loci do not evolve independently due to linkage disequilibrium. Asexual organisms have fewer dimensions of genomic variability as loci are locked together, whereas unlinked loci in sexual organisms are more orthogonal as they can be inherited independently. The independence of loci is also impacted by epistasis; alleles that interact in their effects on fitness will tend to be inherited together. Decomposition of genomic variation via principal component analysis is a familiar concept when analyzing genomic divergence (for recent examples, see Hu et al. 2019;Morales et al. 2019;Tusso et al. 2021). The genetic variance-covariance matrix (the G-matrix) describes standing genetic variation and can also be decomposed into a smaller number of orthogonal axes. There is evidence that widespread pleiotropy makes the dimensionality of the G-matrix lower than the number of measurable traits, but also that pleiotropy causes fitness effects of unmeasured traits to influence the selection observed (e.g., Sztepanacz and Blows 2017a,b). The dimensionality of the G-matrix needs to be considered in the context of the fitness landscape to determine the dimensionality that is relevant here: there will be genomic axes of variation that are neutral, as well as those that are relevant to selection but not to local adaptation. The G-matrix is both a product of mutation and (multidimensional) selection (e.g., Matuszewski et al. 2014) and a determinant of the short-term response to selection (MacPherson et al. 2015).
In the long term, measures of genomic dimensionality are dynamic because genomic architecture can evolve. Where multiple selection pressures are correlated and gene flow is present, high levels of recombination may be costly as haplotypes containing adaptive alleles for different selection pressures are disrupted (Felsenstein 1981;Kirkpatrick and Barton 2006). As the cost of recombination between co-adapted alleles increases, genetic architectures that reduce recombination and thereby lower genomic dimensionality are likely to be favored (Yeaman 2013). The extreme of this would be the formation of supergenes where many previously independent loci underlying a set of coevolving traits become tightly associated (Thompson and Jiggins 2014). This has the effect of rewriting the genetic variance-covariance matrix and redefining the dimensionality of orthogonal traits (Svensson et al. 2021). Alternatively, it may be beneficial to increase genomic dimensionality and break associations between loci, for instance if adaptation to a newly available multidimensional niche requires separation of two associated phenotypes into independent traits.
The genomic dimensionality relevant to local adaptation is likely to increase with divergence dimensionality, but this relationship is not necessarily strong. A broad genomic response, involving many loci (polygenic), could be termed "multidimensional," whereas selection on a single locus might be termed "unidimensional" (Kinsler et al. 2020), regardless of the number of selection pressures to which they respond. On average, one might expect that multidimensional divergent selection elicits multidimensional genomic responses (Nosil et al. 2009a), but this is not guaranteed in every case. A single large-effect locus might pleiotropically affect adaptation of multiple phenotypes to multiple environmental axes, as with the cricket "body size" ( = head width + femur length) example (MacPherson et al. 2015). In contrast, many loci might contribute to divergence on a single phenotypic axis. If the divergence axis is multidimensional in a space defined by orthogonal axes of genetic variation, then the genetic basis of divergence must be more complex than for unidimensional selection in this space.
The extent to which the mapping of loci to traits influences local adaptation has been explored in studies of restricted pleiotropy (Chevin et al. 2010;Le Nagard et al. 2011;Kinsler et al. 2020;Yamaguchi and Otto 2020). MacPherson et al. (2015) show that local adaptation increases more strongly with dimensionality if genetic and environmental axes of variation are correlated. Strong selection on a few loci might overcome gene flow more readily but divergence due to many loci of small effect is possible and may, ultimately, result in a stronger barrier to gene flow (Flaxman et al. 2014;Nosil et al. 2017). Therefore, an impact of dimensionality on local adaptation could be mediated by its effect on genetic or genomic complexity but the extent of this contribution remains largely an open theoretical and empirical question.

AN ALTERNATIVE VIEW ON EXPERIMENTS
Experimental speciation (experimental evolution of diverging populations) is a direct way to test these proposed effects on local adaptation and speciation, and yet no study has explicitly varied the dimensionality of divergent selection. The view that experimental studies support the role of multidimensional selection in promoting speciation dates back to a classic review of speciation experiments (Rice and Hostert 1993), and is repeated in a later review of niche dimensionality (Nosil and Harmon 2009). However, we argue that there was not, nor is there currently, any strong experimental evidence to support this conclusion.
Experimental speciation studies test aspects of ecological speciation by varying the environmental conditions for adaptation (Fry 2009). Therefore, they are well-suited to examine how environmental dimensionality shapes the speciation process. This is in contrast to artificial selection studies in which the experimenter determines fitness based on traits (e.g., Koopman 1950), and so might examine trait dimensionality. Unfortunately, the overwhelming majority of experimental speciation studies have selected along only one axis (although, as established above, selection along one environmental axis might produce divergence in multiple traits). Just five (Kilias et al. 1980;Rice 1985;Rice andSalt 1988, 1990;Rundle 2003) out of 59 studies reviewed by Nosil and Harmon (2009) used multidimensional selection. That these studies were associated with higher levels of reproductive isolation has been taken as evidence that multidimensional selection promotes speciation (Rice and Hostert 1993;Nosil and Harmon 2009).
However, all five of these multidimensional studies used Drosophila species and three involved the same experimental setup, Rice's "habitat maze" (Rice 1985;Rice andSalt 1988, 1990). This design deliberately selects for multiple-effect traits (traits that are under divergent selection and impact other components of reproductive isolation-sometimes referred to as "magic traits"; Smadja and Butlin 2011) because habitat choice is experimentally tied to mate choice. Note that here, in the context of multiple-effect traits, we use the term "traits" in the sense of individual phenotypes rather than orthogonal phenotypic axes. In these experiments, multiple-effect traits are significantly more likely to have driven reproductive isolation than multidimensional selection (Fry 2009;White et al. 2020). This conclusion is reinforced by a recent experimental speciation study of the parasitic feather louse, Columbicola columbae, which produced rapid reproductive isolation via unidimensional divergent selection (via host body size) on a multiple-effect trait (louse body size; Villa et al., 2019). Since the Nosil and Harmon review, few studies have imposed divergent selection (Sharon et al. 2010;Castillo et al. 2015;Markov et al. 2016;Bush et al. 2019) and hence there remains little experimental evidence for the role of multidimensional divergent selection in local adaptation and speciation (White et al. 2020).

EVIDENCE FROM NATURE
In natural populations, several studies cite divergence in multiple phenotypes (Nosil and Sandoval 2008;Gompert et al. 2013;Egea-serrano et al. 2014;Stankowski et al. 2015;Stuart et al. 2017;Aguirre-Liguori et al. 2019). This is to be expected: any two habitats will often generate multiple different demands on an organism, although the extent to which these represent orthogonal axes of selection is unclear. For example, the marine to freshwater transition in sticklebacks changes the osmoregulatory environment, food availability, and predation pressure (Hendry et al. 2013), whereas the coastal to inland transition in monkeyflowers alters water availability, season length, and competition (Lowry and Willis 2010). In Littorina winkles, environmental axes of selection can be combined in different ways, but each axis imposes selection on multiple traits (Morales et al. 2019). There are examples where single traits dominate divergence, such as coloration in Timema stick insects (Sandoval and Crespi 2008) or beach mice (Steiner et al. 2007), or heavy metal tolerance in plants (Singh et al. 2016). There are good examples of studies in which divergence in environments, traits, and genomes has been parsed. For instance, a study of Phrynocephalus lizards identified divergent clusters along the first principal component axis from nine environmental variables, and along the first principal component axis of 11 phenotype measurements. These major axes of environment and phenotype were significantly correlated, indicating divergence along one major environmental/trait axis (Hu et al. 2019). However, systematic comparisons among studies to identify associations between dimensionality and patterns of divergence are rare.
Perhaps the strongest comparison across studies of variable dimensionality is a meta-analysis in which the dimensionality of environmental divergence (one to four traits) was estimated consistently across 35 plant reciprocal transplant studies (MacPherson et al. 2015). The degree of local adaptation was significantly correlated with the dimensionality of divergent selection, which accounted for 20% of the variance. Only 4% of the variance in local adaptation was explained by the overall extent of divergence in environmental measurements, using the same set of sites (Hereford 2009). However, this should not be interpreted as a separation between dimensionality and total selection: MacPherson et al. (2015) acknowledge that the main driver of the dimensionality effect that they observe is likely to be an increase in total selection with more environmental dimensions, with the response to selection probably aided by alignment of environmental and genetic axes. A study of Scutellaria plants similarly concluded that speciation-with-gene-flow and species coexistence are facilitated by multidimensional selection, on the basis that no single environmental factor among 71 variables assessed could account for niche separation. Rather, multidimensional selection partitioning niches along several axes was more consistent with the observed niche distribution among species (Huang et al. 2017).
Quantifying dimensionality of present-day environmental divergence is a difficult issue but it is also necessary to reconstruct it at the time of divergence (Orsini et al. 2012;Pfrender 2012). Although there are good examples of quantification for the present, the past is rarely considered (Öhlund et al. 2020). This is problematic because theory shows that where multiple species coexist at evolutionary equilibrium, a single key change (unidimensional divergence) can destabilize evolutionary equilibria and so cause co-evolutionary ripple effects that lead to multidimensional divergence in response to biotic selection pressures, conflating cause and effect (Gilman et al. 2012;Chevin et al. 2014;Débarre et al. 2014;Yamaguchi and Otto 2020). A study of present-day conditions could thus implicate multidimensional selection, whereas this was not the primary cause of divergence (Öhlund et al. 2020). Furthermore, it is important to consider that present-day divergence on one environmental/trait axis might be the end point of many possible adaptive trajectories with implications for genomic parallelism and transgressive incompatibilities (Chevin et al. 2014;Thompson et al. 2019).

THEORY AND SIMULATION
Upon first impression, most modeling/simulation studies addressing this question indicate that multidimensional selection promotes speciation. However, these have not yet made important distinctions between effects of the number of traits, the overall strength of selection, and the genomic basis of adaptation. For instance, a model of a mosaic metapopulation with complex spatial structure and environmental heterogeneity showed that local adaptation is strongly correlated with dimensionality, and that its impact increases with migration (MacPherson et al. 2015). However, the effects of multidimensionality cannot be separated from overall selection as there was no test of the diluted alternative. Because the model was not genetically explicit, it could not test the effects of number of loci. Furthermore, models generally make the simplification that selection and trait dimensionalities are equal (Nosil and Harmon 2009;Chevin et al. 2014;Thompson et al. 2019), which prevents separation, and that loci are universally pleiotropic.
Beyond local adaptation, three simulation studies have examined reproductive isolation upon secondary contact following divergent selection of variable dimensionality in allopatry. Each study used orthogonal traits where each trait corresponded to fitness in a respective environmental dimension, hence dimensionality was varied along joint environment-trait axes. In one study, axes of divergent selection were added sequentially, with total selection increasing under the additive assumption (Nosil and Harmon 2009). Reproductive isolation, measured as the decrease in the average fitness of hybrids compared to perfectly adapted parental individuals, increased with higher dimensionality. The authors noted that the same effect might be achieved by increasing separation of optima on a single axis, but they did not attempt to isolate the effects of dimensionality from those of overall selection strength. The number of loci was held constant, with all loci influencing all traits. Therefore, selection per-locus also increased at higher dimensionality.
Two subsequent simulations based on Fisher's geometric model probed these issues more deeply (Chevin et al. 2014;Thompson et al. 2019). The overall dimensionality of selection was varied, rather than only divergent selection. The authors made an important distinction between sources of hybrid incompatibilities. Chevin et al. (2014) show that transgressive incompatibilities arise as recombination creates new combinations of alleles causing hybrid phenotypes to vary, not just along the axis separating adaptive optima but also on other axes, reducing fitness under stabilizing selection. Alternatively, "constitutive incompatibilities" can arise as some combinations of derived alleles are incompatible for reasons unrelated to the environment, further reducing hybrid fitness. Both transgressive and constitutive incompatibilities fit the Dobzhansky-Muller incompatibility (DMI) model in a general sense because they depend on negative interactions between derived alleles. Thompson et al. (2019) extend this to include the angle of divergence between two populations' adaptive trajectories, finding that transgressive phenotypes are more common with greater angles of divergence.
The transgressive incompatibility component depends on the set of mutations accumulated during divergence (Fig. 2). This is determined by the evolutionary trajectory of the two populations and not the distance between optima. As the dimensionality of the environment increases, more combinations of alleles in hybrids can generate phenotypes of low fitness under stabilizing selection. Thus, reproductive isolation increases at higher overall environmental dimensionality when all else is equal. Furthermore, under these assumptions, mutation accumulation between two diverging populations increases with the number of environmental dimensions. If some proportion of these mutations produce environment-independent fitness reductions when combined in hybrids, then constitutive incompatibilities are also expected to increase with the environmental dimensionality (Barton 1983;Chevin et al. 2014;Thompson et al. 2019).
Both the Chevin et al. and Thompson et al. models are sophisticated and provide an elegant picture of how environmental dimensionality can affect hybrid fitness in a variety of ways. However, there are important limitations. They do not vary divergence dimensionality, only the overall dimensionality of the environment. Critical assumptions (particularly many loci, each capable of influencing all traits and with free recombination) may be violated, potentially leading to trait-specific effects or an impact of genomic architecture. Furthermore, gene flow between diverging populations has not yet been considered in any model of dimensionality (Nosil and Harmon 2009;Chevin et al. 2014;Thompson et al. 2019). Most significantly, however, no model has yet distinguished the contribution of increasing genomic dimensionality from the effects of dimensionality of divergent selection.

Finding the Way Forward in Multiple Dimensions
Are there ways in which divergence dimensionality might impact local adaptation and speciation other than through increased overall selection, transgressive incompatibility, or increased genetic complexity of adaptation? Divergent selection can drive evolution beyond local adaptation toward speciation. Speciation may result purely from increasing local adaptation, perhaps enhanced by ecological character displacement, resulting in low fitness of migrants and hybrids, and so strong extrinsic barriers to gene exchange. However, for the progression and completion of speciation, secondary barriers to gene flow are likely to be required. How might multidimensionality affect the evolution of these additional barriers?
One possibility is that divergence along more phenotypic axes increases the chance that a trait under divergent selection also contributes to another component of reproductive isolation, such as assortative mating (i.e. that the set of traits includes a multiple-effect trait). Multiple-effect traits reduce gene flow both by reducing the production of hybrids and by reducing hybrid fitness and so they generate strong barriers to gene flow that are less prone to disruption by recombination (Gavrilets 2004;Servedio et al. 2011;Smadja and Butlin 2011). Other barriers may arise as pleiotropic effects or via indirect selection on loci linked to locally adaptive loci, including DMIs and prezygotic barriers such as assortative mating (Rundle and Nosil 2005;Maan and Seehausen 2011). Alternatively, these barriers might arise as an adaptive response, as in reinforcement (Smadja and Butlin 2011). Coupling of these different barriers can produce strong reproductive isolation (Butlin and Smadja 2018;Kulmuni et al. 2020) that can feed back to increase local adaptation by reducing the effect of gene flow. The likelihood that these additional barriers evolve is higher when overall selection is strong, hence they may be more likely to arise under additive multidimensional selection, but this is also affected by the dimensionality of the genomic re-sponse. With more loci under selection, it is more likely that loci underlying secondary barriers are impacted by divergent natural selection (as described for constitutive DMIs in Chevin et al. 2014).
Linkage disequilibrium between local adaptation loci and loci underlying secondary barriers is a key component of ecological speciation. Some of these effects have been described by unidimensional divergent selection baseline models of divergence hitchhiking that explore the impact of selection on multiple loci on divergence at a neutral locus. Where dimensionality is additive (fixed per-locus selection), more loci under selection lead to greater divergence in neutral loci (Feder and Nosil 2010;Flaxman et al. 2012), but under a diluted mode, the opposite is true as selection is applied more weakly across loci, many of which are unlinked with the neutral locus (Barton 1983;Flaxman et al. 2012). Understanding how this aspect of "genomic dimensionality" interacts with the dimensionality of divergent selection may radically alter our conclusions about multidimensional selection. However, thus far it has not been addressed, with models favoring the simplifying assumption that a fixed number of loci underlie adaptation regardless of environmental dimensionality (Nosil and Harmon 2009;Chevin et al. 2014).
Divergence in the presence of gene flow and the consequences of gene flow following secondary contact must also be considered because both are common in nature (Nosil 2008a,b;Seehausen et al. 2014;Schilling et al. 2018). In the presence of gene flow, divergent selection per-locus must be strong enough to build or maintain differences in allele frequencies to create local adaptation and reproductive isolation. If high dimensionality is associated with high numbers of loci and weak selection perlocus, it may impede local adaptation with gene flow. Gene flow also makes the spread of alleles with constitutive incompatibilities more difficult, because the incompatibilities are exposed to selection (Bank et al. 2012). However, gene flow also provides the opportunity for recombination, enabling genetic architecture to influence the evolution of local adaptation and speciation and it tends to align the G-matrix with environmental variability (Beldade et al. 2002). A greater number of loci increase the potential for linkage disequilibrium to build among loci to a point where indirect selection contributes to divergence and creates a strong barrier to gene flow (Flaxman et al. 2014;Nosil et al. 2017). It also increases the chances for linkage disequilibrium between loci under divergent selection and those underlying secondary barriers to gene exchange (Blanquart et al. 2012;Akerman and Bürger 2014;Seehausen et al. 2014;Tigano and Friesen 2016), so enhancing coupling and overall reproductive isolation (Butlin and Smadja 2018). Different demographic scenarios, such as divergence in allopatry, secondary contact, island models with migration or divergence along a cline, along with different genetic architectures, might interact to a greater or lesser extent with di-vergence dimensionality to alter the probability of local adaptation and speciation.
Arriving at satisfactory understanding of these issues will require more theoretical work, more comparative work, and metaanalyses performed using systems whose environmental dimensionality can be accurately measured (MacPherson et al. 2015;Muschick et al. 2020). However, to tease apart the processes and mechanisms of multidimensional selection, experimental speciation studies are also critical, taking up where Rice and others started. Using an experimental speciation approach, total strength of selection, number of dimensions of selection, migration, and population structure can be either controlled or manipulated under laboratory conditions (Fry 2009;White et al. 2020). It is less feasible to control for genomic effects, although evolve and resequence strategies (Kofler and Schlötterer 2014;Schlötterer et al. 2015) can be used to characterize them (Michalak et al. 2019) and it may be possible to choose traits that are likely to have more or less complex genetic architectures.

Concluding Remarks
Moving forward, there are several issues that require greater clarity. First, it remains necessary to arrive at a consistent concept and language of dimensionality. Taking any case in isolation, the additive versus diluted argument is irrelevant; stabilizing selection occurs over n dimensions and habitats vary on m of them, generating divergent selection with some total strength. However, when the effect of dimensionality is discussed, what is really being addressed? Are local adaptation and speciation more likely for greater n (niche dimensionality) or greater m (divergence dimensionality), all else being equal? Does genetic dimensionality matter? Alternatively, is it only the overall strength of divergent selection that matters? The available theory confirms that the overall strength of divergent selection is important, mainly through extrinsic isolation, and also shows that niche dimensionality can contribute to reproductive isolation through the transgressive and constitutive components of incompatibility (Nosil and Harmon 2009;Chevin et al. 2014;Thompson et al. 2019). However, it does not address whether divergence dimensionality or genetic dimensionality (i.e., the numbers of loci available for local adaptation or the dimensionality and orientation of the G-matrix) can also contribute. Covariation between these factors may be common, but it is not inevitable. Full understanding of mechanisms requires the isolation of each of these possible effects. Evidences from experiments, natural populations, and modeling all point to the significance of multidimensional selection but there is work to do to understand the mechanisms involved.