## Abstract

Determining how genetic variance changes under selection in natural populations has proved to be a very resilient problem in evolutionary genetics. In the same way that understanding the availability of genetic variance within populations requires the simultaneous consideration of genetic variance in sets of functionally related traits, determining how genetic variance changes under selection in natural populations will require ascertaining how genetic variance–covariance (**G**) matrices evolve. Here, we develop a geometric framework using higher-order tensors, which enables the empirical characterization of how **G** matrices have diverged among populations. We then show how divergence among populations in genetic covariance structure can then be associated with divergence in selection acting on those traits using key equations from evolutionary theory. Using estimates of **G** matrices of eight male sexually selected traits from nine geographical populations of *Drosophila serrata*, we show that much of the divergence in genetic variance occurred in a single trait combination, a conclusion that could not have been reached by examining variation among the individual elements of the nine **G** matrices. Divergence in **G** was primarily in the direction of the major axes of genetic variance within populations, suggesting that genetic drift may be a major cause of divergence in genetic variance among these populations.

## 1. Introduction

The effect of selection on the genetic variance in natural populations is central to several key questions at the interface between ecology and evolution. At a geographical scale, one explanation for the presence of restricted species ranges invokes the exhaustion of genetic variance in traits under selection in marginal habitats (Hoffmann & Blows 1994). Within populations, life-history trade-offs ultimately must be manifested as a lack of genetic variation in trait combinations involving different components of fitness (Lande 1982). In behavioural ecology, the maintenance of female preference under good genes models of sexual selection is critically dependent on the maintenance of genetic variance in male traits, and associated fitness benefits (Rowe & Houle 1996; Kirkpatrick & Barton 1997). The dynamics of genetic variance under selection is therefore associated with attempts to explain some of the most basic ecological observations concerning the distribution of species, and the life history and behaviour of individuals within populations.

How genetic variance changes under selection is a question that has vexed evolutionary genetics for the past 30 years (Lande 1979; Turelli 1985; Barton & Turelli 1987; Johnson & Barton 2005). Genetic variance in natural populations appears ubiquitous (Lynch & Walsh 1998), and yet it is maintained in the presence of strong directional and stabilizing selection (Kingsolver *et al*. 2001) that is expected to deplete it. Compounding this seemingly contradictory situation is the presence of high levels of genetic variance in life-history traits closely associated with fitness (Houle 1992), when genetic variance for fitness itself is predicted to be low under Fisher's fundamental theorem of natural selection. Although there are plausible mechanisms that may account for the maintenance of so much genetic variance (Turelli & Barton 2004; Burger 2005), it is hard to escape the conclusion that we have yet to come to an empirical understanding of the most basic aspects of how selection affects genetic variance in natural populations.

Evidence for how selection changes genetic variance is difficult to obtain, given that quantitative genetic designs usually need to be applied in more than one population to address this question (Agrawal *et al*. 2001). Responses to selection in artificial selection experiments for unidirectional increases in trait values have plateaued in some cases, suggesting a depletion of genetic variance (Hill & Caballero 1992; Blows & Hoffmann 2005). Mixed-model approaches to the analysis of selection experiments, in which pedigrees are known throughout the experiment, hold considerable promise for directly assessing changes in genetic variance during selection. However, these approaches have been applied in a few cases (Meyer & Hill 1991; Beniwal *et al*. 1992*a*,*b*; Heath *et al*. 1995), and analytical problems concerning the estimation of genetic variances in each generation of selection remain unresolved (Walsh & Lynch 1999).

Empirical investigations of spatial structure in genetic variances of traits under natural selection, such as the change in genetic variance along clines, are surprisingly rare (Barton 1999), particularly given the enormous number of studies devoted to assessing variation among populations in neutral markers. Clinal patterns in realized heritability (Blows & Hoffmann 1993), and temporal changes in genetic variance (Merila *et al*. 2001), suggest that selection may decrease genetic variance in some cases. However, selection may also increase genetic variance, as found in a natural selection experiment that replicated an increase in genetic variance in response to reinforcing natural selection found in field populations (Blows & Higgie 2003). This last example highlights the complexity of the problem; selection can either increase genetic variance (e.g. if favoured alleles segregate at low frequency in the starting population) or deplete it, depending on the underlying genetic details of the response to selection (Barton & Turelli 1987; Barton 1999).

One important factor contributing to our inability to find associations between selection and genetic variance in natural populations is likely to have been an empirical reliance on univariate measures of both selection and genetic variance (Blows & Hoffmann 2005; Johnson & Barton 2005; Blows 2007). It can be readily shown that univariate measures of these microevolutionary parameters can seriously mislead when traits are genetically correlated or experience correlational selection (Pease & Bull 1988; Walsh 2007). The extent to which we have been misled concerning the ubiquity of genetic variance is the subject of recent efforts to determine the dimensionality of genetic variance–covariance (**G**) matrices (Mezey & Houle 2005; Hine & Blows 2006; McGuigan & Blows 2007), and generally it seems likely that a number of trait combinations (dimensions) comprising a particular **G** matrix may lack detectable levels of genetic variance (Kirkpatrick 2008).

If we are to therefore understand how selection changes genetic variance in natural populations, we need not only to know what traits are under selection, but also to assess levels of genetic variance in a multivariate context. For example, multivariate approaches to associating either the direction of linear (Blows *et al*. 2004; Hine *et al*. 2004; Van Homrigh *et al*. 2007) or stabilizing sexual selection (Hunt *et al*. 2007) to **G** matrix orientation have suggested that persistent sexual selection at least may be very effective at depleting genetic variance. One limitation of such studies, however, is that they draw associations between the direction of selection and the orientation of genetic variance *within* a single population, and do not attempt to explain naturally occurring differences in genetic variance *among* populations.

Given the dependency of the change in genetic variance under selection on generally unknown genetic details (Barton & Turelli 1987), theoretical investigations of how **G** matrices evolve under genetic drift and selection have tended to take a simulation approach to the evolution of **G** (Reeve 2000; Jones *et al*. 2003, 2004, 2007). Much of the statistical focus on **G** matrix evolution has been directed at determining whether two such matrices display significant differences. Such approaches are designed to test specific aspects of matrix similarity by using vector forms (Roff 2000), a hierarchy of subspace and proportionality hypotheses (Phillips & Arnold 1999) or random projections (Cheverud & Marroig 2007). While further development of such hypothesis testing frameworks is badly needed (Steppan *et al*. 2002), we lack a robust framework that allows the variation in multiple **G** matrices to be simultaneously characterized. For example, in what trait combination has genetic variance diverged most among populations for a given set of traits? Can we identify independent components of divergence in genetic variance, which may have different underlying genetic causes, such as might be expected if different loci contribute to adaptive responses at different positions along clines (Barton 1999)? To address such questions, we advocate that a geometric understanding of how multiple **G** matrices vary is required to fully characterize among-population heterogeneity in patterns of genetic variance–covariance.

In this paper, we develop a tensor-based approach to the characterization of variation among populations in **G** matrices, which captures all the variance among populations in genetic covariance by maintaining the geometric integrity of the individual **G** matrices. Multiple **G** matrices can be considered as random second-order tensor variables, the variation among which can be characterized by a fourth-order genetic covariance tensor, **Σ _{G}**. We show how divergence among populations in genetic covariance structure in tensor form can be decomposed into independent components of divergence in genetic variance, which constitute the eigentensors of

**Σ**. Crucially, expression of the variation among

_{G}**G**matrices in tensor form allows divergence in genetic variance among populations to be naturally associated with divergence in selection acting on those traits using key equations from evolutionary theory (Lande & Arnold 1983; Zeng 1988). We illustrate our approach using estimates of

**G**matrices and the individual fitness surfaces for male sexually selected contact pheromones (cuticular hydrocarbons or CHCs) from nine geographical populations of

*Drosophila serrata*sampled along a latitudinal cline along the Australian east coast.

## 2. Material and methods

### (a) Tensors as a framework for characterizing divergence in G matrices

In asking questions about how either fitness surfaces (Arnold *et al*. 2001; Rundle *et al*. 2008) or **G** matrices (Roff 2000; Blows & Higgie 2003) vary among natural populations, we are addressing questions in which the data are now second-order random variables: second-order response surfaces in the case of individual fitness surfaces and second-order covariance matrices in the case of **G** matrices. Tensors from multilinear algebra extend the notion of vectors (first-order tensors) and matrices (second-order tensors) to higher-order structures that can be used to characterize variation in these lower-order variables. Although tensors have had limited application in evolutionary biology (Rice 2002, 2004), they are likely to become increasingly important in evolutionary genomics (Alter & Golub 2005; Omberg *et al*. 2007). Here, we adapt techniques for the singular value decomposition of fourth-order tensors (Basser & Pajevic 2007), which allow the characterization of variation in random second-order tensor variables such as **G** matrices.

To begin with, consider the second-order tensor, the **G** matrix that summarizes the variances and covariances among a set of traits within a single population. Each element of **G** is denoted by two indices (table 1), e.g. G_{xy} is the genetic covariance between traits x and y. Similarly, when there is more than one population, to summarize the variances and covariances among the set of **G** matrices, we require a fourth-order covariance tensor, **Σ _{G}**. The tensor must be fourth order because covariances of covariances are denoted by four indices, e.g.

**Σ**

_{G:xy,wz}is the covariance between G

_{xy}and G

_{wz}.

In practice, we can represent the fourth-order covariance tensor as a covariance matrix (**S**) of dimension *n*(*n*+1)/2, which is the number of unique elements in a **G** matrix for *n* traits. Figure 1 shows how to map the elements of **Σ _{G}** to

**S**, dividing the matrix into four quadrants. The top left quadrant contains variances of variances along the diagonal and covariances of variance pairs in the off-diagonal elements. The bottom right quadrant contains variances of covariances along the diagonal and covariances of covariance pairs on the off-diagonal, each multiplied by two. The last two quadrants are mirror images containing the covariances between covariance and variance elements of the original

**G**matrices, each multiplied by √2. The scaling factors 2 and √2 are necessary to map the appropriate number of occurrences of the various terms in the higher-order

**Σ**to the lower-order

_{G}**S**(Basser & Pajevic 2007).

It has long been appreciated that **G** matrices could be decomposed into a set of eigenvalues and eigenvectors, where each eigenvector represents a linear combination of traits that contains variance independent of the variance in any other eigenvector (Dickerson 1955; Hill & Thompson 1978; Lande 1979; Pease & Bull 1988). Similarly, **Σ _{G}** can be decomposed into a set of eigenvalues and second-order eigentensors. The

*n*(

*n*+1)/2 eigenvalues of

**Σ**, are equal to the eigenvalues of its representative

_{G}**S**matrix (Basser & Pajevic 2007). The

*n*(

*n*+1)/2 elements of the

*k*th eigenvector of

**S**:(2.1)where 1≤

*a*<

*b*≤

*n*, can be rearranged to form an

*n*×

*n*matrix to give the

*k*th eigentensor of

**Σ**,(2.2)where, for example, and return the variance among the

_{G}**G**matrices in the genetic variance G

_{11}and the variance in the genetic covariance G

_{12}, respectively, which are explained by the

*k*th eigenvector of

**S**. Since the variances and covariances forming

**S**are calculated directly,

**S**will be a positive semi-definite matrix, i.e. all its eigenvalues will be greater than or equal to zero. By contrast,

**G**will often have negative eigenvalues due to sampling error as a result of the estimation process (Hill & Thompson 1978). The maximum number of non-zero eigenvalues possible for

**Σ**is the smallest of

_{G}*p*−1 and

*n*(

*n*+1)/2, where

*p*is the number of

**G**matrices, from different geographical populations, for example, included in the analysis.

The eigentensors arranged in the form of (2.2) represent mutually orthogonal aspects of how the original **G** matrices have diverged and can be interpreted in a fashion similar to the original **G** matrices. Depending on the sampling design from which the original **G** matrices were taken, the eigentensors of **Σ _{G}** have the potential to have a number of interesting biological interpretations. For example, the change in genetic variance represented by different eigentensors might represent different sources of selection acting on the original traits in natural populations. Alternatively, if the

**G**matrices have been sampled along a latitudinal cline, it is possible that different eigentensors might represent allelic clines at different loci that respond to selection at different positions along the cline (Barton 1999).

Each eigentensor itself can be subjected to a spectral decomposition to determine its eigenvectors and eigenvalues . Therefore, within the space of independent change in genetic variance among the constituent **G** matrices represented by an eigentensor, the eigenvectors describe genetically independent linear combinations of traits that have experienced a change in genetic variance. Although eigenvectors within an eigentensor are linearly independent, the changes in genetic variance in trait combinations they represent are not independent of each other. For example, if an eigentensor has a number of eigenvalues of substantial size, the independent change in the pattern of genetic covariance represented by this eigentensor is spread across a number of genetically independent trait combinations.

Alternatively, if the leading eigenvector of an eigentensor has a substantially larger eigenvalue than the remaining eigenvectors, this would suggest that the independent change in the pattern of genetic covariance represented by this eigentensor is explained primarily by a change in genetic variance in a single trait combination.

To this point, we have characterized how a set of **G** matrices vary using the genetic covariance tensor represented by **S**. It may then be desirable to determine the level of genetic variance in each population for the trait combinations that vary the most in genetic variance among populations. For example, it may be of interest to determine how genetic variance changes along a geographical cline (Barton 1999) in a set of functionally related traits. The analysis of such cases would be greatly facilitated by concentrating on the major independent aspects of divergence in genetic variance, rather than individual trait genetic variances.

In the simple situation of a single **G** matrix, characterizing the genetic variance in a given trait combination (** y**) is achieved by projecting

**through the**

*y***G**matrix using

*y*^{T}

**G**. The principle of linear projection can also be applied to determining the level of genetic variance within each population for those trait combinations that have diverged among populations to the greatest extent. However, there is a complication in using linear projections with

*y***Σ**, because it is possible that the divergence in genetic variance in a given combination of traits could be the result of multiple independent changes in the pattern of genetic covariance, which are accounted for by different eigentensors.

_{G}To understand how populations vary in relation to one particular independent change in genetic variance, it is necessary to partition the genetic variance within each population into parts associated with each eigentensor of **Σ _{G}**. To this end, the

*j*th population

**G**matrix can be expressed as a linear combination of the

*n*(

*n*+1)/2 eigentensors of

**Σ**,(2.3)where

_{G}*C*

^{j,k}are the coordinates of

**G**

_{j}in the basis of the eigentensors of , and can be calculated as the Frobenius inner product (the sum of the element-wise products, yielding a scalar) of and

**G**

_{j}. Since the describe ways in which a set of

**G**matrices differ, similar values of a given

*C*

^{j,k}for any pair of

**G**matrices would suggest that the two matrices do not differ from each other in the way described by the particular under consideration. Another way to view the

*C*

^{j,k}values is that their squared values represent proportions of the squared length of

**G**

_{j}in the space of the

*n*×

*n*

**G**matrices. Therefore, for an individual

**G**matrix,

**G**

_{j},

*C*

^{j,k}indicates how much genetic variance, combined across all traits, is involved in the change captured by the

*k*th eigentensor. To determine the change in genetic variance for a particular trait combination that results from an eigentensor, a final step then is to apply the projection,(2.4)which gives the genetic variance for any given trait combination (

**) in the**

*y**j*th population that has been captured by the

*k*th eigentensor. If the trait combination represented by

**is an eigenvector of , this value is equal to .**

*y*### (b) Associating variation in selection with divergence in G matrices

The variation among populations in linear selection can be characterized by the covariance matrix among linear selection gradients (Felsenstein 1988; Zeng 1988) as(2.5)where the diagonal elements represent the variance among populations of the directional selection gradients (** β**) acting on the

*n*th trait, and the off-diagonal elements represent the covariance of directional selection gradients among populations for pairs of traits. The principle of linear projection can be applied to the association between variation in selection among populations and the orientation of

**Σ**. In this case, however, rather than projecting a first-order tensor (e.g. the trait combination

_{G}**) onto a second-order covariance tensor (e.g. ), a second-order covariance tensor (**

*y***B**) describing the variation among populations in linear selection is projected onto the fourth-order

**Σ**.

_{G}The variance among **G** matrices that can be explained by variance among ** β**s is equal to the length of

**B**in the space of

**Σ**. The Frobenius normalized

_{G}**B**, hereafter denoted as

**B**′ (or for an individual element of

**B**′), is calculated as:(2.6)which is then projected through

**Σ**(Basser & Pajevic 2007, eqn (26)),(2.7)where : represents the Frobenius inner product. In practice, this can be implemented by employing only scalar and vector values using (Basser & Pajevic 2007, eqn (29))(2.8)where

_{G}*λ*

_{i}is the

*i*th eigenvalue of

**B**′; is the

*j*th eigenvalue of the

*k*th eigentensor of

**Σ**;

_{G}

**η**_{i}is the

*i*th eigenvector of

**B**′; is the

*j*th eigenvector of the

*k*th eigentensor of

**Σ**; and

_{G}*σ*

_{k}is the

*k*th eigenvalue of

**Σ**. Note that Basser & Pajevic (2007) inconsistently defined

_{G}*σ*

_{k}to be the square root of the

*k*th eigenvalue of

**Σ**in the text preceding their eqn (11), but we use the definition presented in their nomenclature table. The projection value can be compared with the sum of the eigenvalues of

_{G}**Σ**to quantify how much of the variation among

_{G}**G**matrices shares a subspace with (and may therefore be attributed to) variation in linear selection.

### (c) Predicting divergence in G matrices using microevolutionary parameters

Just as the pattern of genetic covariance can have a substantial impact on the response of population means to selection (Lande 1979), it will also influence how **G** evolves. The within-generation change in **G**, which is predicted to result from a single generation of selection is given by (Lande & Arnold 1983; Phillips & Arnold 1989)(2.9)where ** γ** is the matrix of second-order sexual selection gradients and

**−**

*γ*

*ββ*^{T}is the curvature of the adaptive landscape (Estes & Arnold 2007). The presence of

**G**in (2.9) indicates that the pattern of genetic covariance will bias the direction of change in genetic covariance. For example, if selection acts on a trait combination that has zero genetic variance, genetic variance in this trait combination cannot (mathematically) become non-zero due to selection alone. To determine the potential impact of divergence in the predicted response of

**G**to selection, the variation among populations in Δ

**G**can also be characterized by the fourth-order tensor,

**Σ**. By applying (2.9) using an estimate of the ancestral

_{ΔG}**G**matrix for the set of

**G**matrices under consideration, the predicted change in

**G**generated by the specific selection regimes, characterized in the estimates of

**and**

*γ***, may then be compared with the observed changes in**

*β***G**. It should be emphasized that this approach is based on the assumption that the underlying distribution of allelic effects conforms to Lande's Gaussian model, which may not hold, particularly when genetic variances are observed to change substantially as a consequence of selection (Barton & Turelli 1987).

## 3. Results

### (a) Geographical variation in fitness surfaces and G matrices

Nine geographical populations of *D. serrata* were sampled over a 1450 km latitudinal range along the east coast of Australia, and laboratory populations were established (details in Chenoweth *et al*. 2008). The nine populations in this paper are numbered in order from north to south, where population 1 is the most northerly population (Cooktown, fig. 1 in Chenoweth *et al*. 2008). Sexual selection on male CHCs, arising from female mate preferences in each population, was characterized using standard two-stimulus mate choice tests, where a single female chose between two males (details in Rundle *et al*. 2008). Trait values from all populations were first standardized to the grand mean of the nine populations and a standard deviation of 1. Choice of scale can have a substantial influence on multivariate analyses (Hansen & Houle 2008). We chose to standardize to the grand mean of all nine populations because differences in population mean were important to retain in analyses that used the major axis of the variance–covariance matrix among population means (Lande 1979) to help in interpretation of the genetic covariance tensor. We chose to standardize trait variances as large differences among individual trait variances can cause convergence problems in the estimation of **G** using restricted maximum likelihood. Under this scale, no individual trait dominates the eigenstructure of individual **G** matrices.

Sexual selection gradients for each population were estimated using linear and quadratic regression (Lande & Arnold 1983). Genetic covariance in male CHCs was estimated from nine half-sib breeding designs conducted within each population four generations after they were established in the laboratory (details in Chenoweth *et al*. 2008; eigenstructure of these nine **G** matrices is given in the electronic supplementary material). We constrained each of the population **G** matrices to be positive semi-definite by applying a factor-analytic covariance structure at the sire level, and fitting all eight possible dimensions (Hine & Blows 2006).

### (b) Orientation of the G matrix and linear selection within multiple geographical populations

Linear projections have been used in this system to dissect the association between the direction of linear selection and the orientation of a **G** matrix within a population (Blows *et al*. 2004; Hine *et al*. 2004; Van Homrigh *et al*. 2007). These analyses have suggested that little genetic variance remains in the direction of linear sexual selection. We determined the generality of this association in *D. serrata* using the nine replicate populations (table 2). In eight of nine populations, the level of genetic variance in *β*_{j} was less than what would be expected if the **G** matrix was spherical. On average, *β*_{j} accounted for only 4.9 per cent of the genetic variance (calculated as the percentage of the trace of **G**, tr(**G**)), where it would be expected to be 12.5 per cent if the distribution of genetic variance was equal in all directions.

We conducted a randomization test to determine whether the orientation of *β*_{j} across the nine populations was associated with less genetic variance than expected by chance. For each population, 1000 random normalized ** β** vectors were constructed using random selection gradients drawn from a distribution that had the mean and variance of the 72 observed selection gradients from the nine populations. Genetic variance in the direction of each random vector was then determined by linear projection through

**G**

_{i}. The average genetic variance in

*β*_{j}as a proportion of tr(

**G**

_{i}), for each of 1000 sets of nine projections was then calculated, and this distribution was compared with the observed average genetic variance in

*β*_{j}. The observed level of genetic variance in

*β*_{j}(4.9%) fell between the first and second lowest estimates from the 1000 random estimates of the average genetic variance in

*β*_{j}, corresponding to a two-tailed Pr-value of 0.002. Therefore, the level of genetic variance in

*β*_{j}was significantly lower than expected at random given the observed orientation of the population

**G**matrices.

### (c) Divergence in G matrices characterized by the genetic covariance tensor

The eight eigenvalues of the covariance tensor **Σ _{G}**, for which non-zero values were possible, differed substantially in magnitude (table 3). The first eigentensor accounted for 70.5 per cent of the divergence among the

**G**matrices of the nine natural populations. This suggests that there is a single independent pattern of change in genetic covariance among the eight traits to which most of the divergence in genetic covariance among populations can be attributed. The second eigentensor accounted for a further 17.1 per cent of the divergence in genetic variance (table 3).

(see table 10 in the electronic supplementary material) was dominated by a single trait combination. The dominant eigenvector of had an eigenvalue of −0.977, the square of which (0.95) represents the proportion of the squared length of (equal to 1 because is normalized) that this vector explains. Therefore, 95 per cent of the variance represented by was explained by , indicating that the pattern of change in genetic covariance among populations represented by this eigentensor can be ascribed to a change in genetic variance in a single linear combination of the original traits. The loadings of (table 4) are of the same sign for all CHCs, a pattern that is indicative of the major axis of genetic variance (*g*_{max}) within eight of the nine populations (see the electronic supplementary material). By contrast, (see table 10 in the electronic supplementary material) had two eigenvectors with eigenvalues of substantial size and opposing signs. The first eigenvector accounted for 53.9 per cent of the variance represented by , and the second eigenvector accounted for 40.9 per cent of the variance. Therefore, the independent genetic change represented by occurs mainly in two genetically independent trait combinations; as genetic variance is increased in one trait combination within a population, it is decreased in the other trait combination.

The level of genetic variance (*C*^{j,k}) captured by each eigentensor in each of the nine population **G** matrices can be found in table 5 for and . The variation among **G** matrices accounted for by , and consequently, in the trait combination represented by , was dominated by particularly high levels of genetic variance within each of the two most northerly populations. By contrast, the variation among **G** matrices captured by was primarily caused by a peak in genetic variance in population six.

As a consequence of the large sampling variances known to be associated with individual **G** matrices (Hill & Thompson 1978), **Σ _{G}** will have substantial sampling error associated with its estimation. To assess the effect of sampling error of the constituent

**G**matrices on the eigenstructure of

**Σ**, we conducted a bootstrap resampling procedure that sampled sires within each population (with replacement) and then generated 100 bootstrapped replicates of the nine

_{G}**G**matrices. Only 100 bootstrap replicates were employed, given the computational requirements of estimating 900

**G**matrices employing restricted maximum likelihood. We then constructed

**Σ**for each of the 100 bootstrap replicates.

_{G}To determine the consistency with which the first two eigentensors of **Σ _{G}** were represented in the 100 bootstrap replicates, we projected the observed and into each bootstrapped

**Σ**to calculate the variance among

_{G}**G**matrices captured by each eigentensor. Figure 2 shows the relationship between the proportion of variance among bootstrapped

**G**matrices captured by (

*x*-axis) and (

*y*-axis). All projection values for and were above zero, indicating that the original eigentensors captured non-zero differentiation among the hundred bootstrapped sets of

**G**matrices. The projection value of through the bootstrapped

**Σ**was bigger than that of for all but three bootstrap replicates (figure 2). This suggests that the structure of

_{G}**Σ**is preserved to some extent among bootstrap replicates.

_{G}We found no clear association between the change in phenotypic mean of trait combinations along the latitudinal cline and the change in genetic variance captured by the first eigentensor of **Σ _{G}** (figure 3). For the trait combination represented by , the mean and genetic variance both reduce in a seemingly co-coordinated fashion in the three most northerly populations, but then change in an uncorrelated fashion thereafter (figure 3

*a*). Similarly, for the second eigentensor, the peak in genetic variance in the sixth population, uncovered for (figure 3

*b*), was not associated with a change in trait mean, and a substantial change in trait mean in population 4 was not associated with a change in genetic variance. However, for , the mean and genetic variance peak in the most northerly region in a co-coordinated fashion (figure 3

*c*).

An alternative approach to investigating the association between trait means and genetic variance is to determine the change in genetic variance in a particular trait combination of interest. In clinal studies, for example, the trait combination that has diverged the most in mean among populations may be of particular interest in characterizing clinal adaptation. This trait combination, *d*_{max}, is estimated as the major axis of the variance–covariance matrix of population means (Lande 1979; Zeng 1988; Blows & Higgie 2003; McGuigan *et al*. 2005). Here, *d*_{max} is sharply nonlinear, and the rapid increase in mean in northerly populations is matched by an increase in genetic variance (figure 4). The change in genetic variance in *d*_{max} among the populations was tested for significance by applying the vector to generate a univariate trait. We then implemented a mixed model (standard half-sib nested random effects) that allowed each population to have separate estimates of genetic variance in *d*_{max}, and then compared the fit of this model with the one in which all populations shared the same estimate of genetic variance in this trait using a log-likelihood ratio test. The model allowing for separate estimates of genetic variance for each population was a significantly better fit to the data (*Χ*_{9}^{2}=30.35, Pr<0.001).

The vector *d*_{max} can be projected onto using (2.4) to determine the change in genetic variance in this trait combination among populations accounted for by each eigentensor (figure 4). The large change in mean in *d*_{max} in the two most northerly populations is associated with an increase in genetic variance that is shared between and . Interestingly, the increase in variance in *d*_{max}, as a consequence of both eigentensors, is because *d*_{max} is highly correlated with (vector correlation of 0.90), but not completely orthogonal from (vector correlation of 0.48).

### (d) Associating variation in selection with divergence in G matrices

A geometric approach to assessing the impact of variation in selection on the variation in **G** matrices can be implemented by using equation (2.8) and projecting **B**′ through **Σ _{G}**. The projection of

**B**′ through

**Σ**resulted in a variance of 0.0368, suggesting that variation in linear selection accounted for only 0.74 per cent of the divergence in

_{G}**G**matrices. The same randomization procedure described in §3

*b*above was used to generate a distribution of 1000 random estimates of

**B**′, and indicated that the observed level of divergence in genetic variance in the direction of

**B**′ was significantly less than that expected at random (Pr=0.006, two-tailed test).

### (e) Predicting divergence in G matrices using microevolutionary parameters

We generated **Σ**_{ΔG} by first calculating Δ**G** for each population using the average within-population **G** matrix (**G**_{w}) from Chenoweth & Blows (2008), and then using these nine matrices to construct **S**. Ideally, an independent estimate of **G** from the ancestor of all nine populations would have provided the base from which to assess divergence in **G** (Steppan *et al*. 2002). Implicit in our use of **G**_{w} is the assumption that the population phylogeny is star-like in shape. This assumption has some empirical support, as neutral marker divergence among these populations suggested only weak population structure, and an overall pattern of isolation by distance (Chenoweth & Blows 2008). Therefore, **G**_{w} was considered a reasonable proxy for the ancestral **G** in the absence of further information about the phylogeny of these populations.

The first and second eigentensors of **Σ**_{ΔG} (see table 11 in the electronic supplementary material) accounted for 77.9 and 19.4 per cent of the variation in **Σ**_{ΔG}, respectively (table 3). The eigenanalysis of resulted in a single dominant eigenvalue (; table 4). All CHCs loaded onto with the same sign and similar magnitude, again consistent with the direction of greatest genetic variance (*g*_{max}) in CHCs within populations (see tables 1–9 in the electronic supplementary material). The eigenanalysis of yielded two dominant eigenvalues of opposing sign (table 4), reflecting that when genetic variance is predicted to increase in the direction of , it is predicted to decrease in the direction of (and vice versa) as a consequence of the independent change in response to genetic variance represented by .

Finally, we determined whether the eigenstructure of **Σ**_{ΔG} predicted how the **G** matrices had diverged among populations. The projection of these two eigentensors through **Σ _{G}** resulted in variances of 1.722 and 0.940, respectively, together explaining 53.8 per cent of the variance in

**Σ**. We assessed the significance of the association of divergence in the predicted change in

_{G}**G**with the observed divergence in

**G**by generating random

**matrices by drawing the diagonal and off-diagonal elements from separate distributions with their own mean and variance calculated from the observed selection gradients, in addition to generating random**

*γ***vectors as above, and then applying equation (2.9) to each replicate. In this way, we generated a null distribution of the association between**

*β***Σ**

_{ΔG}and

**Σ**, in which there was no natural orientation of selection acting on

_{G}**G**. The observed projection was smaller than all random projections in this distribution, equating to a probability level of Pr<0.001. Therefore, even though variation in

**Σ**

_{ΔG}explained over 53 per cent of the variation in

**Σ**, this was still significantly less than expected at random.

_{G}## 4. Discussion

### (a) The use of the genetic covariance tensor in uncovering geographical patterns in genetic variance

In this paper, we have provided a geometric framework using higher-order tensors, which enables the empirical characterization of how **G** matrices have diverged among populations. This tensor-based approach has two advantages in describing **G** matrix divergence. First, the entirety of divergence in **G** matrices is elegantly captured, and can be decomposed into independent aspects of divergence (the eigentensors), which in turn can be interpreted in a fashion similar to the original **G** matrices themselves. The eigentensors, therefore, hold considerable potential to enhance our understanding of the genetic basis of adaptation in natural or experimental populations. Second, divergence in **G** matrices is characterized while maintaining the geometric integrity of the individual population matrices, which is a crucial determinant of the response to selection within each population. As a consequence, the tensor framework lends itself naturally to further investigations of how microevolutionary processes may have led to the observed divergence in **G** matrices using the key equations for evolutionary change.

In the example of **G** matrix divergence with which we have illustrated this method, the analysis of the genetic covariance tensor **Σ _{G}** provided a number of insights into the evolution of genetic variance in CHCs of male

*D. serrata*. Much of the divergence in genetic variance occurred in a single trait combination (the first eigenvector of the first eigentensor), a surprising conclusion that could not have been reached by examining variation among the individual elements of the population

**G**matrices. Furthermore, this major axis of divergence among

**G**matrices was very similar to the major axis of genetic variance (

*g*_{max}) within eight of the nine geographical populations investigated here, and in previous work (Blows

*et al*. 2004). Although a theory of how

**G**matrices diverge as a consequence of genetic drift is yet to be fully developed (Phillips

*et al*. 2001), the level of genetic variance represented in

**G**has been predicted to decrease in proportion to the inbreeding coefficient in a population (Lande 1980

*b*). Hence, the largest changes in genetic variance generated by genetic drift would be expected in those directions with the greatest genetic variance. While our finding that the first eigenvector of the first eigentensor is similar to

*g*_{max}within populations is consistent with genetic drift, we cannot exclude selection from playing a role in this association, as the response of the genetic variance to selection may be greatest for those trait combinations that initially contain the greatest genetic variance.

It is worth highlighting here a potential point of confusion between the commonly employed *Q*_{ST}−*F*_{ST} method of inferring adaptive genetic change and the genetic covariance tensor. In *Q*_{ST}−*F*_{ST} comparisons, the among-population variance in phenotypes raised under common-garden conditions is referred to as the additive genetic variance among populations (Merila & Crnokrak 2001), and in its multivariate form, the among-population **G** matrix (Kremer *et al*. 1997; Chenoweth & Blows 2008). This method partitions genetic differences in mean within and among populations, where the within-population variation is subsequently partitioned into genetic and non-genetic variance components in a breeding design. By contrast, the genetic covariance tensor allows the characterization of among-population variation in genetic variance, and not simply genetic differences in mean. The projection of *d*_{max} through **Σ _{G}** essentially returns how genetic variance in the major axis of the among-population

**G**matrix from Chenoweth & Blows (2008) differs among populations. What is clear from the application of the two approaches to this set of nine populations is that the major axes of genetic differentiation among populations do not equate to the major axes of divergence in genetic variance. This can be attributed to the well-known complexity of the association between the mean and genetic variance (Barton & Turelli 1987), a point we now explore in more detail below.

Our approach had mixed success in uncovering any clear associations between the change in phenotypic mean and genetic variance along the latitudinal cline from which the populations were sampled. For example, in a theoretical investigation of the behaviour of genetic variance along clines, Barton (1999) found that for a change in genetic variance, as a consequence of allele frequency change, the ratio of the widths of the change in mean and genetic variance should indicate the number of loci involved. We found no obvious association between mean and genetic variance for the major axis of divergence captured by the first eigentensor of **Σ _{G}**. This may be because no such association between the mean and genetic variance would be apparent for that part of the divergence in genetic variance that is consistent with genetic drift. In contrast, there was some association between a substantial change in mean in northern populations and an increase in genetic variance in the trait combination that had diverged most in mean among the populations. Since the divergence in mean in this trait combination is likely to be associated with some form of climatic selection on CHCs (Chenoweth & Blows 2008), this may explain why an association between mean and genetic variance is apparent in this trait combination. It is important to note, however, that the associations we have drawn between genetic variance and the putative forces of genetic drift and climatic selection are not independent. That is, the change in genetic variance in

*d*_{max}was captured by both the first and second eigentensors of

**Σ**, and therefore, as suggested above, genetic drift and selection may be operating in non-orthogonal directions.

_{G}Theory indicates that the extent to which genetic variance will change under selection depends heavily on the genetic details of the selection response. In the presence of a large number of loci of small effect under selection, genetic variance will change little (Lande 1980*a*). However, with fewer genes, some of which have relatively large effect, the change in genetic variance can be substantial (Barton & Turelli 1987; Reeve 2000). The substantial change in genetic variance in *d*_{max} is consistent with a relatively simple genetic basis, although given that the increase in genetic variance occurred across only two of our populations, it is difficult to ascertain with these data if Barton's (1999) prediction of relative widths of change in mean and variance supports this view. Furthermore, the presence of any such gene of major effect is not sufficient for a substantial difference in genetic variance to remain among populations. Fixation of the favoured allele would result in only a transient effect on the **G** matrix (Agrawal *et al*. 2001), and therefore any major gene that contributes to *d*_{max} is unlikely to have gone to fixation in these northern populations.

### (b) Microevolutionary predictions of divergence in genetic variance

The microevolutionary analysis of the evolution of **G** using the pattern of divergence in linear sexual selection was spectacularly unsuccessful in explaining virtually any of the divergence in **G** matrices. Instead, random patterns of selection were significantly more likely to explain more of the divergence in **G** matrices than current patterns of sexual selection. Such an extreme association between the orientation of sexual selection and the divergence in **G** matrices may be attributed to how sexual selection affects genetic variance in each population. There is, on average, a non-random association between the linear selection gradients and the **G** matrices within these nine populations, reinforcing a number of experiments that have indicated that persistent linear sexual selection is effective at depleting genetic variance (Blows *et al*. 2004; Hine *et al*. 2004; Van Homrigh *et al*. 2007). Since the directions of linear selection in groups of these populations are associated to some extent based on whether they were allopatric or sympatric with *Drosophila birchii* in the field (figure 1*a*; Rundle *et al*. 2008), the depletion of genetic variance tends to be in similar directions within these two groups of populations. Therefore, the directions in which linear selection operates are unlikely to contribute substantially to the variation among populations in **G** as these directions generally have low genetic variance. In summary, the overriding consequence of sexual selection on male CHCs is to deplete genetic variance in the trait combinations under selection, and since the direction of sexual selection tends to be similar among groups of populations, sexual selection does not generate divergence in the genetic covariance structure.

A complete microevolutionary prediction of the divergence in **G**, however, requires taking into account the orientation of the ancestral **G**, as the evolution of **G** may be constrained considerably by the pattern of genetic covariance itself. Although we have set out how this may be accomplished using **Σ**_{ΔG}, there are a number of caveats that apply to the particular dataset we have used to illustrate the method. Most importantly, although more than half the variance in **Σ _{G}** was explained by the first two eigentensors of

**Σ**

_{ΔG}, this was actually less than what would be expected given random selection matrices acting on

**G**

_{w}. Therefore, the strength of the association between

**Σ**

_{ΔG}and

**Σ**can be considered spurious.

_{G}The seeming ability to predict a reasonable proportion of divergence among **G** matrices is likely to be an artefact of two related causes. First, the calculation of Δ**G** by equation (2.9) constrains Δ**G** to be in the space of **G**_{w}. Most of the genetic variance is restricted to a subspace of **G**_{w}, with the first four eigenvectors of **G**_{w} accounting for over 94 per cent of the estimated genetic variance within populations (Chenoweth & Blows 2008). Second, given that the orientation of **G**_{w} is similar to the population **G** matrices, particularly with regard to the major axis of genetic variance (and hence little genetic variance remains in the direction of selection as shown above), it is hardly surprising that the divergence in predicted response captured by **Σ**_{ΔG} is heavily biased in the direction of this major axis of genetic variance. So, our substitution of **G**_{w} for ancestral **G** and the use of estimates of sexual selection, which are strongly associated with the depletion of genetic variance, are key limitations of the specific application we present here. Ideally, an independent estimate of ancestral **G**, coupled with estimates of selection that are relevant to divergence along a latitudinal cline (e.g. climatic natural selection), would represent a better system to investigate the ability of microevolutionary parameters to predict divergence in **G** along a latitudinal cline.

### (c) Conclusion

The division between microevolutionary processes that generate evolutionary change, and the macroevolutionary patterns in divergence among taxonomic groups in nature, has yet to be fully bridged (Arnold *et al*. 2001). At the level of traits means, retrospective selection analyses have attempted to determine whether the observed difference in trait means could be predicted by the multivariate breeder's equation when it was assumed that divergence in trait means equated to the historical selection gradient (Lande 1979; Turelli 1988). Extensions of this approach, based on the theory of population divergence in multiple traits (Felsenstein 1988; Zeng 1988; Hansen & Martins 1996), have allowed the observed divergence in population means to be predicted from microevolutionary estimates of both selection and genetic variance (Chenoweth *et al*. submitted). Finally, the integration of phylogenetic comparative methods with quantitative genetic theory (Felsenstein 2008; Hohenlohe & Arnold 2008) promises further insights into the roles of selection and genetic constraints on population mean divergence. Tensor-based approaches to the characterization of divergence in genetic covariance structure have the potential to play a central role in the development of a comparative quantitative genetic framework for understanding divergence in genetic variance.

## Acknowledgments

We are grateful to Peter Basser and Austin Lund for advice on tensors, Russ Lande and an anonymous reviewer for comments on the manuscript and the Australian Research Council for funding.

## Footnotes

One contribution of 14 to a Theme Issue ‘Eco-evolutionary dynamics’.

- © 2009 The Royal Society