## Abstract

In contrast to our growing understanding of patterns of additive genetic variance in single- and multi-trait combinations, the relative contribution of nonadditive genetic variance, particularly dominance variance, to multivariate phenotypes is largely unknown. While mechanisms for the evolution of dominance genetic variance have been, and to some degree remain, subject to debate, the pervasiveness of dominance is widely recognized and may play a key role in several evolutionary processes. Theoretical and empirical evidence suggests that the contribution of dominance variance to phenotypic variance may increase with the correlation between a trait and fitness; however, direct tests of this hypothesis are few. Using a multigenerational breeding design in an unmanipulated population of *Drosophila serrata*, we estimated additive and dominance genetic covariance matrices for multivariate wing-shape phenotypes, together with a comprehensive measure of fitness, to determine whether there is an association between directional selection and dominance variance. Fitness, a trait unequivocally under directional selection, had no detectable additive genetic variance, but significant dominance genetic variance contributing 32% of the phenotypic variance. For single and multivariate morphological traits, however, no relationship was observed between trait–fitness correlations and dominance variance. A similar proportion of additive and dominance variance was found to contribute to phenotypic variance for single traits, and double the amount of additive compared to dominance variance was found for the multivariate trait combination under directional selection. These data suggest that for many fitness components a positive association between directional selection and dominance genetic variance may not be expected.

CHARACTERIZING the genetic basis of phenotypes and the form and resulting consequences of selection on these phenotypes is a major goal of evolutionary biology. Substantial effort has been devoted to estimating additive genetic variance in metric traits and fitness components (Falconer and Mackay 1996; Lynch and Walsh 1998), establishing that the majority of metric traits have additive genetic variance and a heritability in the range of 0.2–0.6 (Lynch and Walsh 1998). More recently, the necessity of examining multivariate patterns of additive genetic variance and selection has been emphasized (Walsh and Blows 2009) and shown to influence the multivariate response to selection in laboratory (*e.g.*, Mezey and Houle 2005; McGuigan and Blows 2009; Hine *et al.* 2014) and natural populations (Clements *et al.* 2011; Morrissey *et al.* 2012). In particular, additive genetic variance in all single traits often does not equate to genetic variance in all multivariate trait combinations (Hine and Blows 2006; Blows 2007; Walsh and Blows 2009), and a response to selection in trait combinations with low levels of additive genetic variance may be stochastic in nature (Hine *et al.* 2014).

In contrast to our growing understanding of patterns of additive genetic variance, very little is known about the relative contributions of nonadditive variance, particularly dominance variance, to multivariate phenotypes, despite the ubiquity of dominance and the attention given to the evolution of dominance for close to a century (Fisher 1928; Wright 1929; Wright 1934). The pervasiveness of dominance is demonstrated most clearly by inbreeding depression (Charlesworth and Charlesworth 1987), which is exhibited in almost all organisms to some degree (Husband and Schemske 1996; Lynch and Walsh 1998). The increased frequency of homozygous loci across the genome, caused by inbreeding, exposes recessive deleterious alleles that are typically held in a heterozygous state, resulting in the decrease in trait means, and initial fitness reductions that accompany inbreeding (Lande and Schemske 1985; Lynch 1991; Charlesworth *et al.* 1999). Similarly, selection for recessive deleterious alleles in heterozygotes results in the typically faster responses observed for downward than upward artificial selection, when the dominant allele confers a higher trait mean than the recessive (Falconer and Mackay 1981).

Remarkably high levels of inbreeding depression, averaging 50%, have been demonstrated for primary fitness components (*e.g.*, viability, fertility, egg production) in *Drosophila*, compared to levels of a few percent for morphological traits (Lynch and Walsh 1998), and asymmetric selection responses for fitness components have also been demonstrated in several cases (Frankham 1990). These observations are consistent with predicted patterns of dominance genetic variance for traits that are genetically correlated with fitness (Fisher 1930; Frankham 1990; Crnokrak and Roff 1995). Directional selection on fitness-correlated traits is expected to erode the additive genetic variance in these traits (Fisher 1930), resulting in lower heritabilites (Mousseau and Roff 1987; Price and Schluter 1991; Kruuk *et al.* 2000) and a higher proportion of dominance variance contributing to overall phenotypic variance (Merilä and Sheldon 1999). Consequently, fitness itself, a trait unambiguously under directional selection, is predicted to have the least additive and proportionately largest contribution of dominance variance to phenotypic variance. In addition to single traits, particular multivariate combinations of traits known to be under persistent directional selection through their association with sexual fitness have been demonstrated to have low additive genetic variance (McGuigan and Blows 2009; Hine *et al.* 2011; Sztepanacz and Rundle 2012), suggesting that similar patterns of additive and dominance contributions to phenotypic variance may also be observed for multivariate traits in certain circumstances.

In general, we know relatively little about proportions of dominance variance for any single traits, regardless of their association with fitness, and multivariate patterns of dominance for suites of functionally related traits are unknown (Merilä and Sheldon 1999). For single traits in general, compilations of dominance variance estimates by Crnokrak and Roff (1995), and more recently by Wolak and Keller (2014), indicate that dominance variance, as a proportion of total phenotypic variance, is in the range of 0.04–0.36 for undomesticated species. Overall, the current data suggest that dominance variance contributes, on average, 15% of total phenotypic variance of traits and that its contribution is significantly different from zero (Wolak and Keller 2014), supporting an important contribution of dominance variance to total phenotypic variance. These estimates of dominance variance are also consistent with estimates of nonadditive variance obtained by comparing broad-sense to narrow-sense heritability for a variety of traits (Hill *et al.* 2008).

Evidence supporting a larger contribution of dominance variance to phenotypic variance for fitness-correlated traits, compared to traits less correlated with fitness, however, is equivocal. There is some evidence to suggest that fitness-correlated traits have lower heritiabilites than morphological traits (Houle 1992; Kruuk *et al.* 2000), and this pattern does not appear to be due to proportionately smaller additive variances, but is a consequence of larger residual components of variance, which include dominance variance (Houle 1992). Consistent with the prediction of more dominance variance for fitness-correlated traits, life-history traits have previously been observed to have more dominance variance than morphological traits and also a higher proportion of dominance compared to additive variance (Crnokrak and Roff 1995). However, equal contributions of dominance variance to life-history and morphological traits and a similar contribution of dominance and additive variance to total phenotypic variance have been observed in a more recent comprehensive data set (Wolak and Keller 2014).

These conflicting results may be due, in part, to the way that levels of dominance were estimated and how the association of these traits with fitness was determined. For undomesticated species, dominance is typically estimated by inbreeding or by crossing inbred lines (*e.g.*, Hughes *et al.* 2002; Kearsey *et al.* 2003; Bilde *et al.* 2008). Throughout the genome, these processes create levels of homozygosity and heterozygosity, respectively, that are unlikely to be seen in natural populations. Given that dominance variance depends on gene frequencies (Falconer and Mackay 1981), whether estimates of dominance from line crosses are representative of segregating dominance variance is unclear. In addition, life-history traits have been broadly classified in these analyses as fitness related, based on presumed positive phenotypic correlations between these traits and fitness. However, phenotypic correlations between traits and fitness have been demonstrated in several instances to be unrepresentative of underlying genetic correlations, which ultimately govern patterns of genetic variance (Rausher 1992; Stinchcombe *et al.* 2002; Kruuk *et al.* 2002).

To unequivocally determine whether traits correlated with fitness exhibit a higher proportion of dominance variance and lower proportion of additive variance than traits more distantly related to fitness, a more direct method is to simultaneously examine the additive genetic covariance of traits with fitness and their respective dominance variance in outbred populations (Merilä and Sheldon 1999). It has been recognized, particularly for additive genetic variance, that neither single-trait variances nor bivariate covariances are likely to reflect the true nature of genetic variation (Lande and Arnold 1983; Blows and Brooks 2003; Blows and Hoffmann 2005) and that multiple functionally related traits are often the targets of selection (Lande and Arnold 1983; Phillips and Arnold 1989). While genetic correlations between fitness components may maintain additive variance in each trait due to contrasting selection on these traits, additive variance in total fitness can still be close to zero (Walsh and Blows 2009). Given that the relative contributions of additive and dominance variance to overall phenotypic variance are inextricably linked, a comprehensive understanding of patterns of dominance variance will require a multivariate approach.

Here, we used a large, multigenerational breeding design in an unmanipulated population of *Drosophila serrata* that enabled the estimation of additive and dominance variance for both multivariate wing-shape phenotypes and a comprehensive measure of fitness, within a single framework. Our experimental design enabled the characterization of additive and dominance variance for single morphological traits and fitness and the covariance between them. This allowed us to directly test the hypothesis that traits genetically correlated with fitness should exhibit high levels of dominance variance and the corollary that fitness itself should exhibit the highest level of dominance. By employing a multivariate approach we were not limited to examining these patterns on a trait-by-trait basis, but we also determined how dominance variance is distributed across multivariate trait combinations, and within this multidimensional space we examined the association between additive and dominance variance. Finally, by employing a multivariate form of the Robertson–Price identity (Stinchcombe *et al.* 2013), we were able to determine the levels of additive and dominance variance in the multivariate trait combination under the strongest directional selection, extending tests of the hypothesis, of high dominance variance in fitness related traits, to a multivariate framework.

## Materials and Methods

### Breeding design

Prior to phenotypic measurements, two generations of breeding were conducted using a previously described, outbred, and laboratory-adapted population of *D. serrata* (Hine *et al.* 2014). The first generation of breeding was a paternal half-sibling design in which 75 sires were each mated to three virgin dams. Following a 72-hr mating period, dams were moved to individual vials (day 0) and allowed to oviposit in these vials for 72 hr. On days 12–14, four virgin sons and four virgin daughters from productive half-sibling families were collected within 24 hr of eclosion using light CO_{2} anesthesia. These flies were held in separate sex vials at a density of two to four flies per vial, for 3–5 days prior to the second generation of breeding.

The second generation of breeding employed a double-first-cousin design. While half-sib designs are effective for estimating additive genetic variance, double-first-cousin breeding designs provide the coefficients of fraternity necessary for the estimation of dominance variance (Lynch and Walsh 1998). Here, breeding was done in a circular pattern: two pairs of full brothers from each family were mated to two pairs of full sisters, where the two pairs of sisters were not related to each other or to the mating partners of their brothers and sisters. Single male–female pairs were placed in individual vials and allowed to mate and oviposit for 72 hr. Virgin sons from productive families were collected within 24 hr of eclosion using light CO_{2} anesthesia, for use in either competitive fitness assays or for wing phenotypes. Males to be used for competitive fitness assays were collected on the second day of vial emergence and were held at a density of four flies per vial for 4 days prior to the fitness assay. Males to be used for wing phenotypes were collected on the first, second, and third days of emergence and were stored in familial groups in 1 ml Eppendorf tubes at −20**°**, prior to phenotyping. In total the two-generation pedigree comprised about 9500 individuals, and about 5000 males from 685 families were used for wing phenotypes and about 2800 males from 666 families were used for fitness assays.

### Fitness assay

Male competitive fitness was assayed using five virgin sons from each of 666 families. A single focal male (red eye) from the breeding design was placed in a vial with two randomly chosen females and one randomly chosen male from an outbred *D. serrata* population fixed for a recessive mutation conferring an orange-eye phenotype (day 0).

Females were allowed to freely mate and oviposit for 48 hr, after which, all males and females were discarded from the vials. By day 14 the majority of offspring from these vials had eclosed, and these flies were transferred into clean plastic vials and stored at −20**°**. Subsequently, the numbers of orange-eye (competitor) and red-eye (focal) offspring produced in each vial were counted. The number of adult offspring sired by a male from the breeding design reflects his competitive mating success, the productivity of the female he mated to, and the survival to emergence of his offspring. In the following analyses, male competitive fitness is calculated as the natural logarithm of the odds ratio of focal male to competitor-male offspring sired (Reddiex *et al.* 2013): (1)This metric is biologically meaningful within the context of our experimental design. The ratio reflects relative fitness when both focal and competitor males sire offspring, and in cases where only one male sires offspring, it also reflects the total number he sired. Therefore, these data are more informative than a simple proportion, where the rankings within groups of sires that sire none or all of the offspring are arbitrary. Within the subset of focal sires that produced all (none) of the offspring, those siring a relatively higher (lower) number have a higher (lower) ranking. Taking the natural logarithm transforms this ratio, while preserving the rank of sires, into an approximately normal distribution that is suitable for standard multivariate analyses. This distribution is centered on 0, where the focal male sired an equal number of offspring as the competitor male, with positive and negative values indicating increasingly more focal-male and increasingly more competitor-male offspring sired, respectively.

### Wing phenotypes

To measure wing phenotypes, one wing from each male (either left or right) was removed and mounted on a standard microscope slide using double-sided tape. The number of wings collected from each family ranged from 0 to 15, with a median of eight wings per family, resulting in 5040 wings phenotyped from 685 families. Wings were photographed using a Leica MZ6 microscope with a Leica IC80 HD camera attachment, and nine landmarks corresponding to those previously described for *D. bunnanda* (Mcguigan and Blows 2007) were recorded using tpsDig2 software (Rohlf 2005). To characterize shape variation free from size variation, wings were scaled by centroid size and the landmark coordinates were aligned (Mezey and Houle 2005), resulting in nine X–Y coordinate pairs (18 traits) for each wing. To reduce the number of traits to make analyses computationally feasible, without losing information from regions of the wing, we calculated the Euclidian distance between the landmark coordinates (Mcguigan and Blows 2007). This reduced the 18 coordinate traits to 10 interlandmark distance traits that have been previously described (McGuigan and Blows 2007) and herein are referred to as: aLM1-8, aLM2-4, aLM2-5, aLM2-8, aLM3-7, aLM3-9, aLM4-5, aLM5-6, aLM5-8, and aLM7-9. Due to high multicollinearity between aLM1-8 and aLM2-8, and also between aLM3-7 and aLM7-9, traits aLM1-8 and aLM7-9 were removed from subsequent analyses.

### Statistical analyses

#### Genetic analyses:

A total of 24 multivariate outliers (0.5% of the total wing data) were identified and removed from the wing phenotype data using the multivariate Mahalanobis distance technique implemented in the software package JMP (v. 10, SAS Institute, Cary, NC). For the fitness measures, vials that did not produce any offspring (zeros), likely due to experimental error resulting in female damage or death, were removed prior to analyses. This resulted in 2883 fitness measures from 666 families, with an average of 4.3 fitness measures per family.

Wing phenotypes were measured from flies collected over 3 days of emergence, allowing for a potential effect of emergence day on shape variation. We did find statistically significant differences in the mean value of interlandmark distances from flies that had emerged on different days, for all traits except aLM3-9. However, in all cases the effect of “day” explained <0.8% of the phenotypic variance, indicating that the day effect was not an important biological phenomenon. Therefore, we pooled wings collected on different days, and interlandmark distances and fitness were subsequently standardized globally {∼*N*[0,1]}, prior to analyses.

#### Univariate analyses:

Univariate and multivariate genetic analyses were conducted using restricted maximum likelihood (REML) within an animal model framework, implemented in Wombat (Meyer 2007). The breeding design enabled us to fit a model that estimated two uncorrelated random effects at a time: additive genetic and dominance genetic (or additive by additive epistasis, or common environmental variance). The assumptions for the parameter means and variances were (2)where is the additive genetic variance, is the dominance genetic variance ( is the additive × additive epistatic variance; is the common environmental variance), and is the residual variance. We were able to include only two random effects at a time (additive variance + one other) in our analyses, because the nature of our breeding design, where full sibs were reared in a single vial, resulted in a partial confounding between common environmental variance, dominance variance, and additive × additive epistatic variance. Attempting to fit more than one of these three random effects in a single analysis would be redundant, resulting in a singular average information matrix (also known as Fisher information matrix, or a nonpositive definite Hessian).

The sparse numerator relationship matrix (**A**) was calculated and inverted within Wombat, whereas we calculated and inverted the dominance relationship matrix (**D**) and additive × additive epistasis matrix (**AA**) using the nadiv package (Wolak 2012) in R (R Core Team 2012). In cases examining common environmental variance, the vector of random effects was related to individuals in proportion to a blocked identity matrix. To calculate the dominance relationship matrix, the nadiv package (Wolak 2012) approximates the coefficients of fraternity () between individuals *i* and *j* by (3)where *k* and *l* represent the dam and sire of *i*, *m* and *n* the dam and sire of *j*, and *θ* is the additive genetic relatedness between individuals noted in the subscripts (Wolak 2012). This approximation assumes no inbreeding and ignores any dominance connections through grandparents, which were both inconsequential given the structure of our pedigree. The relationship matrix for additive × additive epistasis was calculated as the matrix product of the numerator relationship matrix with itself. The nonzero elements of the lower triangle of the inverse of these matrices were supplied to Wombat as a general inverse file with an additional codes file linking the running numbers of the individuals in the analysis with their pedigree identifiers.

Univariate analyses of wing traits and fitness employed the following mixed model, (4)where *μ* is the population mean, *a _{i}* is the additive genetic value of individual

*i*,

*d*(

_{i}*aa*) is the dominance (epistatic, common environmental) genetic value of individual

_{i}; ce_{i}*i*, and

*e*is the residual error (no fixed effects were fit for any models). Univariate models were fit for each wing trait and fitness, individually, and statistical significance of each variance component was assessed using a log-likelihood ratio test comparing the fit of a model that included

_{i}*vs.*excluded the variance component of interest, with one degree of freedom. In all cases, estimation was carried out using a strict AI REML algorithm.

Initially, we explored whether there was any evidence for common environmental variance in these traits, finding little evidence that this was an important component of variance (see Supporting Information, File S1). Therefore, we excluded common environment as a random effect from all uni- and multivariate models. Next, we examined whether there was any evidence for additive × additive epistasis. This component of variance is typically ignored in standard quantitative genetic analyses that employ a half-sibling breeding design or parent offspring regression, for example, despite that additive variance in these cases is confounded with additive × additive epistasis. In these cases, however, additive × additive epistasis contributes only one-quarter the amount of additive variance and, therefore, is considered negligible. Unfortunately, double-first cousins and full-sibs, the only standard relationships that contribute to dominance variance, share an equal amount of dominance and additive × additive epistatic variance. An examination of our **AA** and **D** relationship matrices revealed that the coefficients of relatedness for dominance variance were completely nested within those for additive × additive epistasis; therefore, there was no way to separate these components in our analyses (Table 1). However, a thorough investigation comparing the likelihood and precision of models that include epistasis *vs.* dominance indicated that epistasis contributes marginally to phenotypic variance, compared to dominance variance for these traits (see File S1, Figure S1, and Table S1). Therefore, we excluded additive × additive epistasis from our models and subsequently examined additive genetic, dominance genetic, and residual variance only.

#### Multivariate analyses:

For multivariate analyses the more general matrix form of the mixed-model equations were used, (5)where *y* is the vector of observations on all individuals, and are the vectors of additive genetic and dominance genetic effects that are related to individuals proportionally to the numerator (**A**) and dominance (**D**) relationship matrices, respectively. is the vector of residual variance, related to individuals in proportion to the identity matrix.

Parameter constraints limited the number of traits we could include in a multivariate mixed model that simultaneously estimated additive, dominance, and residual variance. The number of estimated parameters increases exponentially with the number of traits included in multivariate analyses, a pattern that is exacerbated when additional random effects, such as dominance, are considered. This generates a more severe upper limit on the number of traits that can be included in a single analysis (Meyer 1992). In this case, the upper limit for the number of traits was driven by the ability to disentangle dominance from residual variance. Therefore, to obtain accurate estimates for all three components of variance, particularly dominance variance, the maximum number of traits we could include in an analysis was two. However, to obtain accurate estimates for additive variance these parameter constraints were not limiting because there was substantial power to disentangle additive variance from the rest of the phenotypic variance. For eight (nine) trait models that included additive, dominance, and residual variance, dominance and residual components could not be adequately separated from each other, rendering dubious estimates of these individual components; however, combined, they accounted for the correct proportion of variance. Additive genetic estimates were unaffected by the confounding of dominance and residual variance, and therefore these models can be thought of as estimating an additive covariance matrix and a joint dominance/residual matrix.

We used these models to determine the number of genetic principal components underlying the additive covariance matrix via a series of nested reduced rank models for both wing traits alone and wing traits with fitness. The number of principal components that underlie the additive genetic covariance matrix indicates whether the genetic association between traits confines genetic variance to certain combinations of these traits or, in other words, how many independent combinations of these traits have significant genetic variance. Covariance estimates from an analysis in which no structure was imposed on the covariance matrix at the additive genetic level were used as starting values for reduced rank models. Statistical support for the genetic dimensions underlying the additive genetic covariance matrix (**A**) was evaluated using nested log-likelihood ratio tests for the reduced-rank models (Hine and Blows 2006). Typically, reduced-rank estimation is carried out using factor analytic modeling (*e.g.*, Hine and Blows 2006; McGuigan and Blows 2010; Sztepanacz and Rundle 2012): these analyses are equivalent to factor analytic models in which all specific variances are assumed to be zero (Meyer and Kirkpatrick 2008).

For each of the models in this set of nested analyses that allowed us to test the statistical dimensionality of **A**, the dominance and residual covariance matrices were fit at full rank. As outlined above, interpreting individual estimates of these components would be inappropriate. In cases in which wings and fitness were estimated in the same model, the residual covariance between fitness and each of the wing traits was fixed at zero to reflect these traits being measured on different individuals. For reduced-rank models, likelihood estimation was carried out using a hybrid algorithm with initial iterates of an expectation maximization (EM) procedure followed by average information (AI) REML. In cases where convergence was difficult to achieve, this process was repeated for 100 cycles. In most instances this attained convergence, with the exception of reduced- rank models fitting zero, one, and two principal components at the additive genetic level. This is likely a consequence of restricting significant variance to such few dimensions (see *Results*). Assuming that all specific variances are zero when considering too few factors can result in biased estimates of residual components, with principal component models typically behaving more poorly for highly restricted ranks than factor analytic models (Meyer and Kirkpatrick 2008).

To simultaneously estimate additive and dominance variance in a single multivariate model, and within our parameter restrictions, we carried out analyses of all bivariate trait combinations using the multivariate mixed model described in Equation 3 and employing a strict AI REML estimation method. For bivariate analyses between fitness and each wing trait, the residual covariance between these traits was fixed at zero, to reflect these traits being measured on different individuals. These analyses appropriately partitioned additive, dominance, and residual variance from each other, yielding accurate covariance estimates for both additive and dominance random effects. Because bivariate analyses were carried out on all pairwise subsets of traits, they resulted in single estimates for the additive and dominance covariances between all traits, but several estimates for the additive and dominance variance in each trait. We used a maximum-likelihood approach to pool the estimates from these analyses, generating full additive and dominance (**D**) covariance matrices for the eight wing traits and wing traits plus fitness. This approach pools covariance matrices for all sources of variation simultaneously, while accounting for the typically strong negative sampling correlations between estimates from analyses with overlapping subsets of traits (Meyer 2012). The pooled additive genetic covariance matrix was nearly identical to the one estimated from the full multivariate model, so for consistency we present results from pooled additive and dominance matrices only, and use the “full” model for testing the significant genetic dimensions underlying the *A*-matrix, only. Although constructing a full dominance covariance matrix by pooling part analyses restricts our ability to establish the statistically significant genetic dimensions underlying this matrix, it does not prohibit us from examining the patterns and orientation of dominance variance in a multivariate context.

#### Matrix comparison:

We employed three approaches to comparing patterns of additive and dominance genetic variance. First, to determine the amount of additive genetic variance in the eigenvectors of **D**, and hence the similarity between the eigenvectors of **A** and **D**, we projected each of these nine vectors through the additive genetic covariance matrix using the equation *d _{i}*

^{T}

**A**

*d*, where

_{i}*d*is the

_{i}*i*th eigenvector of the dominance covariance matrix (scaled to unit length) and

**A**is the additive genetic covariance matrix. Next, to examine whether the parts of the additive and dominance matrices that contained the majority of variance were similar to each other we employed Krzanowski’s common subspace approach (Krzanowski 1979), using the first four eigenvectors of

**A**and

**D**that accounted for 89 and 97% of the variance, respectively. This formal approach for comparing matrices identifies the similarity in orientation of two subspaces using the equation (6)where the matrices

**D**and

**A**contain the subset

*k*of the eigenvectors of the dominance and additive covariance matrices as columns and where

*k*≤

*n*/2 (Krzanowski 1979; Blows

*et al.*2004; Aguirre

*et al.*2014). Here, the sum of the eigenvalues of

**S**indicates the similarity in orientation of the subspaces, with 0 indicating that the subspaces are completely orthogonal and

*k*indicating that they completely overlap (Blows

*et al.*2004). High similarity between these subspaces would indicate that the multivariate traits with the most additive genetic variance are also the multivariate traits with the most dominance variance. Finally, to determine the relative contributions of wing traits

*vs.*fitness in generating the difference between additive and dominance matrices, we employed a genetic covariance tensor approach. For only two matrices the multivariate trait combination describing the greatest difference in genetic variance between them is given by the first eigenvector of the difference matrix (

**A**–

**D**), calculated here by subtracting the dominance covariance matrix from the additive covariance matrix (Hansen and Houle 2008; Sztepanacz and Rundle 2012). This trait combination is equivalent to the leading eigenvector of the first eigentensor of the fourth-order covariance tensor that characterizes variation among replicate matrices (Hine

*et al.*2009). While the Krzanowski subspace comparison indicates how different the multivariate trait combinations with the most additive and dominance variance are overall, the first eigenvector of (

**A**–

**D**) indicates which trait combination is most different between

**A**and

**D**. The trait loadings of this vector indicate the relative importance of each trait in determining the difference in orientation of additive and dominance variance.

## Results

Within a large multigenerational breeding design we measured the fitness of 2883 sons from 666 families (∼4 sons/family) and 5040 wings (∼8 sons/family) from their brothers. We subsequently carried out a genetic analysis using a suite of eight wing phenotypes (interlandmark distance traits: aLM2-4, aLM2-5, aLM2-8, aLM3-7, aLM3-9, aLM4-5, aLM5-6, aLM5-8) and fitness. Consistent with the levels of additive variance as a proportion of total phenotypic variance typically found for quantitative traits (Lynch and Walsh 1998), estimates of additive genetic variance (*V*_{A}) for the eight wing traits analyzed here ranged from 0.29 to 0.56 (Figure 1A and Table 2). Fitness, however, was strikingly different, with a *V*_{A} estimate of 0.008, 2 orders of magnitude less than wing traits (Figure 1A and Table 2). Statistical support for the presence of additive genetic variance in all eight wing traits was provided by univariate likelihood-ratio tests, which tested whether excluding *V*_{A} significantly worsened the fit of the model, with degrees of freedom equal to one (aLM2-4, χ^{2} = 75.7, *P* < 0.001; aLM2-5, χ^{2} = 50.4, *P* < 0.001; aLM2-8, χ^{2} = 101.4, *P* < 0.001; aLM3-7, χ^{2} = 79.0, *P* < 0.001; aLM3-9, χ^{2} = 33.0, *P* < 0.001; aLM4-5, χ^{2} = 68.0, *P* < 0.001; aLM5-6, χ^{2} = 63.8, *P* < 0.001; aLM5-8, χ^{2} = 50.4, *P* < 0.001). This statistical support was retained for all eight traits when the highly conservative sequential Bonferroni correction of the nominated significance level of 0.05 was applied (**α** = 0.05, *c* = 8). For fitness, the univariate test for the presence of *V*_{A} was not statistically significant (χ^{2} = 0; d.f. = 1, *P* = 1) (Table 2).

An eigenanalysis of the wing trait **A** revealed that the first seven (of eight) eigenvectors accounted for 99.7% of *V*_{A} in these traits. Despite the small proportion of variance (0.3%) accounted for by the last eigenvector of the additive matrix, genetic principal component modeling revealed statistical support for all eight genetic dimensions underlying **A** (reducing from eight to seven dimensions significantly worsened the fit of the model: χ^{2} = 14.00; d.f. = 1, *P* < 0.001). Although uncovering a full-rank genetic covariance matrix is uncommon, there is also evidence in *D. melanogaster* that the additive matrix for wing shape is full rank (Mezey and Houle 2005). Consistent with the lack of *V*_{A} for fitness found in the univariate test, statistical support was found for eight of nine genetic dimensions of the wing trait plus fitness **A** (reducing from eight to seven dimensions significantly worsened the fit of the model: χ^{2} = 13.83; d.f. = 1, *P* < 0.001).

For wing traits, estimates of dominance variance (*V*_{D}) were lower, on average, than additive estimates, ranging from 0.029 to 0.317 (Figure 1A and Table 3). Univariate likelihood-ratio tests, testing whether excluding *V*_{D} significantly worsened the fit of the model (again, with degrees of freedom equal to one) provided statistical support for dominance variance in three of the eight wing traits (aLM2-5, χ^{2} = 4.6, *P* = 0.033; aLM3-7, χ^{2} = 4.8, *P* = 0.029; aLM3-9, χ^{2} = 6.5, *P* = 0.011). Applying the sequential Bonferroni correction (**α** = 0.05, *c* = 8) reduced statistical support to zero of eight traits; however, with only eight tests, this likely represents the conservativeness of the test rather than a true lack of significance, with the probability of obtaining a type I error for this family of tests equal to 0.34. Although estimates of *V*_{D} were, in general, only moderately lower than *V*_{A} (Figure 1A), due to the nature of our breeding design, statistical tests of dominance components suffer from reduced power. Here, all information on dominance comes from the 275 double-first-cousin pairs created in the second generation of breeding, 60% less sires than are used for additive genetic estimates. Despite the exceedingly low *V*_{A} for fitness, its dominance estimate of 0.23 was within the range of the majority of *V*_{D} estimates found here for wing traits and was statistically supported in a univariate likelihood ratio test (χ^{2} = 7.12; d.f. = 1, *P* < 0.01).

An eigenanalysis of the pooled **D** for wing traits revealed that the first five (of eight) eigenvectors accounted for 99.5% of the dominance variance, with the last eigenvector accounting for only 5 × 10^{−5%} of the variance in these traits. Although pooling bivariate analyses to generate the “full” dominance covariance matrix precludes us from testing the significant dimensions underlying it, the observation of such a low eigenvalue suggests that a null or nearly-null space may exist within the dominance matrix that does not exist within the additive matrix. Although fitness did have significant *V*_{D} in a univariate test, the amount of variance captured by the first five (of nine) eigenvectors of **D** that include wing traits and fitness was equal to the amount described by the first five (of eight) for **D** that includes only wing traits.

In total there was 47% less dominance than additive genetic variance for wing traits (as given by a trace of 3.66 *vs.* 1.94 for additive and dominance covariance matrices, respectively). When fitness was included, this difference in variance decreased to 41% less dominance than additive variance. Matrix projection of the nine eigenvectors of **D** through the full-rank **A** indicated that there was substantial additive genetic variance in each of these multivariate trait combinations (Figure 2). In fact, the proportion of additive genetic variance contained in the eigenvectors of **D**, and corresponding eigenvectors of **A**, were remarkably similar over most of the space. To determine the proportion of additive and dominance genetic variance in the multivariate trait combination under strongest directional selection we examined the genetic selection gradient, given by the eight-element vector of the additive genetic covariance between each wing trait and fitness (Table 2) (Stinchcombe *et al.* 2013). The projection of this normalized vector through **A** and **D** uncovers the proportion of additive and dominance genetic variance, respectively, in this multivariate trait combination. In contrast to patterns of *V*_{A} and *V*_{D} for the univariate trait under strongest directional selection, fitness, the proportion of *V*_{A} in this multivariate trait combination was 0.31, double the proportion of *V*_{D}, which was equal to 0.16.

The Krzanowski method of subspace comparison identified a 69% similarity in orientation of the subspaces defined by *a*_{max} to *a*_{4} and *d*_{max} to *d*_{4}** _{,}** indicating that the multivariate traits with the most additive genetic variance also have the most dominance genetic variance. This was demonstrated by a value of 2.76 for the sum of the eigenvalues of

**S**, which range here from zero to four, indicating orthogonal and coincident subspaces, respectively. Although the eigenvalues of

**A**decline less rapidly than those of

**D**, with the first four accounting for only 89% of additive variance but 97% of dominance variance in these traits, Krzanowski’s method is limited to comparing subspaces of dimension

*k*≤

*n*/2; therefore, a formal comparison of a larger subspace was not possible within this framework.

Despite the similarity between **A** and **D,** they are not identical, with notable differences in the distribution of variance over their respective eigenvectors shown by the more rapid decline of variance for the eigenvectors of **D** compared to **A** (Figure 1B). To determine which traits contribute most strongly to this difference we examined the leading eigenvector [(**A**–**D**)_{max}] of the difference matrix (**A**–**D**) (Table 4). The trait (**A**–**D**)_{max} accounted for the majority of the difference between **A** and **D** (0.62 of the total difference in genetic variance of 1.42), with the elements, or trait loadings, of this vector indicating the relative contribution of each trait to this difference. Surprisingly, fitness was one of the weakest contributors to this trait combination, despite being the only trait to have significant dominance variance, and essentially no additive variance. The strongest contributor to (**A**–**D**)_{max} was aLM5-6, with a trait loading of −0.605. This trait had a significant, moderate, estimate of *V*_{A} with a relatively low, nonsignificant estimate of *V*_{D}. Overall the trait loadings across (**A**–**D**)_{max} appear fairly even, with positive and negative values indicating that the major difference between **A** and **D** is primarily driven by contrasting contributions of anterior to posterior shape variation in the proximal region of the wing, to proximal to distal shape variation across the wing.

## Discussion

Dominance is widely recognized to be pervasive and has implications for the maintenance of genetic variance (Barton and Keightley 2002) and the genetic dynamics of small populations (Whitlock and Fowler 1999; Turelli and Barton 2006; Taft and Roff 2011) and may influence responses to selection (Keightley 1996). However, data on levels of naturally segregating dominance variance is limited, and the association between dominance variance and fitness is equivocal (Crnokrak and Roff 1995; Roff and Emerson 2006; Wolak and Keller 2014). A direct test of the prediction that dominance variance should be positively correlated with directional selection, and the multivariate extension of this prediction was possible by employing a breeding design that enabled the simultaneous estimation of additive and dominance variance for a suite of functionally related traits and fitness. We demonstrated that fitness itself had low additive and high dominance variance, as predicted. However, both single and multivariate morphological traits were found to have significant additive genetic variance that contributed a similar proportion as dominance variance to overall phenotypic variance, regardless of the association between these traits and fitness.

Consistent with population genetic theory which predicts low additive genetic variance in fitness for populations at equilibrium (Falconer and Mackay 1996), we did not find statistical support for additive genetic variance in fitness, with the *V*_{A} estimate of 0.008 (as a proportion of the total phenotypic variance) 2 orders of magnitude lower than the heritability typically observed for quantitative traits (Lynch and Walsh 1998). Although many fitness components have been shown to have high levels of additive genetic variance, when scaled to reflect their large residual components (Houle 1992), fitness components are often negatively correlated (Falconer and Mackay 1996; Brooks 2000; Chippindale *et al.* 2001) and may be subject to antagonistic selection arising through either direct (Godin and McDonough 2003) or pleiotropic (Hunt *et al.* 2004; McKean and Nunney 2008) effects of underlying alleles on other fitness components. Antagonistic selection may maintain additive genetic variance in fitness components, but additive genetic variance in net fitness can still be close to zero (Charlesworth and Hughes 2000; Walsh and Blows 2009).

How to best measure fitness is a controversial topic (Brommer *et al.* 2004), with its measurement often an empirically challenging task. The *D. serrata* population used here was maintained on a laboratory transfer schedule, in which adult lifespan is only a few days once sexual maturity is reached, for ∼50 generations prior to this experiment. Our fitness assay was, therefore, designed to capture the significant components of fitness under these conditions and was similar to other assays that have been used to characterize fitness in *D. serrata* (*e.g.*, Delcourt *et al.* 2012; Reddiex *et al.* 2013; Collet and Blows 2014). In this species, female remating rates are among the highest in *Drosophila*, with previous evidence indicating that during 48 hr of mating opportunities a single female will produce offspring from 2.9 sires, on average (Frentiu and Chenoweth 2008). Here we have housed two females with a single orange-eyed competitor male and a single red-eyed focal male for 48 hr, with the resulting offspring likely produced by both females that had each mated several times. Therefore, our measure of fitness was comprehensive, including a male’s mating success, the productivity of the female he mated to, and the survival to emergence of his offspring, and it was within a competitive arena that is likely to occur in nature (Frentiu and Chenoweth 2008).

Low additive genetic variance for net fitness has also been demonstrated in other populations of *D. serrata* (Delcourt *et al.* 2012), in red deer (Kruuk *et al.* 2000), and in North American red squirrels (McFarlane *et al.* 2014); however, we cannot exclude the possibility that *V*_{A} exists here for fitness but we lacked the power to detect it. We were able to detect statistically significant variance of a magnitude smaller than 0.008 for the eigenvectors of the wing-shape additive covariance matrix; however, we had reduced power to detect additive genetic variance in fitness for two reasons. First, while wing phenotypes and fitness were collected from approximately the same number of families, there was, on average, half the number of fitness measures compared to wing measures per family. Second, competitive fitness was measured in relation to random, genetically variable competitor males, thus increasing the within-family variation of a trait that is already particularity sensitive to environmental perturbations (Price and Schluter 1991; Whitlock and Fowler 1999).

In contrast to the low estimate of *V*_{A} for fitness, the dominance estimate of 0.253 was moderate and statistically supported, despite the reduced power (275 double-first cousins) we had for detecting dominance variance. Therefore, dominance variance was the predominant component of genetic variance in fitness. Published estimates of dominance variance for fitness are generally lacking, although several agricultural studies have examined *V*_{D} for traits that are under persistent artificial directional selection. These artificially selected traits may be analogous to fitness in natural populations, as they are always under directional selection (with selection on other traits generally reduced), and their trait values primarily determine whether an individual contributes offspring to subsequent generations. The *V*_{D} observed here for fitness is consistent with dominance variance accounting for 14–22% of the phenotypic variance for body weight at harvest of artificially selected trout populations (when excluding common environmental effects) (Pante *et al.* 2002), and with dominance variance accounting for 19–52% of the phenotypic variance for productive life and lifetime in selected cattle populations (Fuerst and Sölkner 1994).

Two alternative scenarios may explain our observation of low additive and moderate dominance variance for fitness. First, if many loci are overdominant with respect to fitness (and overdominant alleles vary in their dominance coefficients), there will be a range of intermediate gene frequencies for which additive variance is close to zero and dominance contributes most to genetic variance (Lynch and Walsh 1998). Pleiotropy may generate overdominance for fitness (Van Dooren 2006), resulting in the maintenance of heterozygotes due to selection on each fitness component (Mitchell-Olds *et al.* 2007); however, in general, selection experiments provide little evidence that overdominance is common (Eisen 1980; Falconer and Mackay 1996; Lynch and Walsh 1998). Although, there is limited evidence in *Caenorhabditis elegans* that induced mutations affecting fitness may exhibit overdominance (Peters *et al.* 2003). Second, if selection has fixed beneficial and purged moderately deleterious alleles, and therefore eroded much of the additive genetic variance in fitness in the population (Fisher 1930), then remaining segregating mutations affecting fitness would be mildly deleterious (Eyre-Walker and Keightley 2007). The dominance variance estimate of 0.253 would then suggest that, on average, these segregating alleles are recessive to some extent. This would be consistent with data from the yeast deletion project, which suggests that the typical segregating mutation affecting fitness is moderately but not completely recessive (Agrawal and Whitlock 2011). Because dominance variance characterizes the departures from additivity for the summed effects of all loci affecting a focal trait (Lynch and Walsh 1998), without knowing the distribution of dominance coefficients for all alleles affecting fitness, distinguishing between mechanisms that may generate dominance variance is difficult.

Except for fitness itself, we did not find any evidence that levels of dominance variance are associated with directional selection for single traits. Dominance variance for wing shape traits was similar to the *V*_{D} observed for fitness and within the range 0–35% of the phenotypic variance previously observed for morphological traits (Wolak and Keller. 2014). However, contrary to the prediction that traits more strongly correlated with fitness should have more dominance variance contributing to phenotypic variance, there was no relationship between the additive genetic covariance of traits with fitness and their dominance variance. One caveat is that the 95% confidence limits (given by the sampling variance) for the additive genetic covariance between each trait and fitness overlapped zero in all cases (Figure S2), indicating that directional selection on these traits may be limited, although the low *V*_{A} estimate for fitness likely weakened our power to detect a significant covariance. Nevertheless, estimates of *V*_{D} for single traits were also comparable to estimates of *V*_{A} (Figure 1A), which were consistent with the additive genetic variance previously observed for wing landmark coordinate positions in *D. melanogaster* (Mezey and Houle 2005). While there is some previous evidence that life-history traits, presumably under stronger directional selection, have more dominance compared to additive variance than morphological traits presumably under weaker directional selection (Crnokrak and Roff 1995; Roff and Emerson 2006), more recent data suggest that this may not be a general pattern (Wolak and Keller 2014). Previous studies examining levels of additive and dominance variance in single traits have relied mainly on data from line crosses to estimate these parameters (*e.g.*, Roff and Emerson 2006). Gene frequencies of 0.5, characteristic of inbred lines, might be expected to result in a higher proportion of dominance *vs.* additive variance for traits, compared to populations that have more dispersed gene frequencies (Lynch and Walsh 1998; Hill *et al.* 2008), providing a potential explanation for the higher proportions of dominance *vs.* additive variance previously observed.

In addition to single traits, we also examined the similarity between additive and dominance genetic covariance matrices for multivariate wing phenotypes, followed by the contributions of additive *vs.* dominance variance for the multivariate trait under directional selection. Multivariate analyses revealed that the distribution of additive genetic variance across the eigenvectors of **A** declined in a typical exponential fashion (Hine and Blows 2006; Kirkpatrick 2009; Walsh and Blows 2009; Sztepanacz and Rundle 2012), and despite the last two eigenvectors accounting for only 0.8 and 0.3% of the additive genetic variance in these traits, principal component modeling revealed significant additive genetic variance in all dimensions. This result is consistent with the full-rank additive covariance matrix observed for wing shape in *D. melanogaster*, in which bootstrapping was used to determine significance (Mezey and Houle 2005). While observing a full-rank genetic covariance matrix is unusual (Kirkpatrick 2009), with only the first few eigenvectors often having statistical support in principal component (factor analytic) models (*e.g.*, Hine and Blows 2006; McGuigan and Blows 2007), our large sample size of 685 sires provided substantial power for detecting small levels of additive genetic variance. We were unable to determine the number of significant genetic dimensions underlying the dominance covariance matrix because it was generated using a series of bivariate models.

The lack of association between directional selection and dominance variance observed for single traits was also observed for the multivariate trait combination under directional selection, ** s_{g}**. Contrary to predictions, this trait combination actually had a higher

*V*

_{A}estimate of 0.31, compared to

*V*

_{D}of 0.16, given by the projections of this vector, standardized to unit length, through the additive and dominance covariance matrices, respectively. This vector of selection differentials is unbiased by environmentally induced covariances between traits and fitness (Robertson 1966; Lynch and Walsh 1998; Delcourt

*et al.*2012); however, predicted evolutionary responses characterized by

**do not distinguish between the effect of selection acting directly on the traits of interest and selection acting through correlated traits (Stinchcombe**

*s*_{g}*et al.*2013). Although

**is under directional selection, the combination of positive and negative trait loadings contained in this vector indicate that additive genetic variance may be maintained by antagonistic selection on the individual traits comprising this multivariate trait combination. In addition, for traits aLM2-4, aLM4-5, and aLM5-8 the genetic selection gradient (**

*s*_{g}**), which scales**

*β*_{g}_{,}**by**

*s*_{g}**A,**was close to zero; however, the genetic selection differentials were 0.46, 0.51, and −0.26, respectively. This indicates that evolutionary responses would occur in these traits as a consequence of indirect selection on correlated traits and not selection acting directly on them (Stinchcombe

*et al.*2013), again indicating that pleiotropy may maintain additive genetic variance in the multitrait combination,

**.**

*s*_{g}In summary, we have found that for fitness, a trait unequivocally under directional selection by definition, dominance variance contributed substantially more than additive genetic variance to overall phenotypic variance. For morphological traits, however, there was no relationship between the strength of directional selection on single or multivariate traits and *V*_{D}, with similar proportions of *V*_{A} and *V*_{D} observed for single traits. Using a direct, quantitative, test of the hypothesis that dominance variance should increase with trait–fitness correlations, in an outbred population, we failed to find any association between dominance variance and directional selection for quantitative traits. A positive association between directional phenotypic selection and dominance variance may, therefore, be an exceptional, rather than general, observation.

## Acknowledgments

We thank Miranda Vidgen, Donna Petfield, Amanda Denton, Grace Stanton, Katrina McGuigan, Julie Collet, Emma Hine, and Jessica McCarroll for help in the lab and David Aguirre for useful discussions concerning the analyses of our data.

## Footnotes

Supporting information is available online at http://www.genetics.org/lookup/suppl/doi:10.1534/genetics.115.175489/-/DC1.

*Communicating editor: C. D. Jones*

- Received February 11, 2015.
- Accepted March 11, 2015.

- Copyright © 2015 by the Genetics Society of America