## Abstract

The genetic basis of stochastic variation within a defined environment, and the consequences of such micro-environmental variance for fitness are poorly understood. Using a multigenerational breeding design in *Drosophila serrata*, we demonstrated that the micro-environmental variance in a set of morphological wing traits in a randomly mating population had significant additive genetic variance in most single wing traits. Although heritability was generally low (<1%), coefficients of additive genetic variance were of a magnitude typical of other morphological traits, indicating that the micro-environmental variance is an evolvable trait. Multivariate analyses demonstrated that the micro-environmental variance in wings was genetically correlated among single traits, indicating that common mechanisms of environmental buffering exist for this functionally related set of traits. In addition, through the dominance genetic covariance between the major axes of micro-environmental variance and fitness, we demonstrated that micro-environmental variance shares a genetic basis with fitness, and that the pattern of selection is suggestive of variance-reducing selection acting on micro-environmental variance.

THERE has been a long history of interest in how genotypes vary in their environmentally generated phenotypic variance (Reeve and Robertson 1953; Waddington 1953; Lerner 1954) and in the fitness consequences of that variation (Schlichting 1986; Gavrilets and Hastings 1994; Wagner *et al.* 1997; Mulder *et al.* 2007; Tonsor *et al.* 2013). Recent analyses of long-term field data from wild populations (Garant *et al.* 2004; Charmantier and Gienapp 2013; Bouwhuis *et al.* 2014; Gienapp and Merilä 2014; Merilä and Hendry 2014) have highlighted the substantial covariance that the environment can generate between fitness and phenotypes (Rausher 1992). However, there remains considerable uncertainty as to whether environmental variance and fitness share a genetic basis, and therefore if environmental variance evolves as a consequence of selection. Improved understanding of these questions is key to understanding how heritable variation in environmental sensitivity might influence evolution (Lande 2009, 2015; Chevin *et al.* 2013).

There are three nested scales of environmental variation from which environmental variance (V_{E}) (as a component of phenotypic variance) can arise. The environmental variance generated by macro-environmental variation has been referred to as plasticity, while that from external and internal micro-environmental variation is referred to as canalization and developmental instability, respectively. While the difference between macro- and microenvironments is an arbitrary division of a continuous scale of environmental effects (Debat and David 2001), determining the genetic basis of macro- *vs.* micro-environmental variance requires different approaches, focused on trait mean and trait variance, respectively. This is because for macro-environmental variance the environmental parameter is known and can be manipulated, whereas for micro-environmental variance it is unknown (Falconer and Mackay 1996). Evidence for genetic variance in plasticity typically relies on analyzing the difference in mean phenotype of particular genotypes across defined macro-environments, with many studies providing clear evidence of genetic variance in phenotypic responses to macro-environmental variation (*e.g.*, Schlichting 1986; Lynch and Walsh 1998; Pigliucci 2005) (Figure 1A). In contrast, evidence for genetic variance in micro-environmental variance (both developmental instability and canalization) requires analyzing the variance in phenotype of particular genotypes within a given macro-environment (Figure 1B). This dependence on variances poses a twofold problem for empirical studies of micro-environmental variance. First, estimating variances requires greater sample sizes than estimation of means because of the greater sampling variance (Fuller and Houle 2002; Hill and Mulder 2010; Visscher and Posthuma 2010). Second, the trait of interest can no longer be adequately described by a single phenotypic measure made on an individual.

There are three general approaches that have been used to estimate the genetic variance in micro-environmental variance. The first two approaches have been predominantly used in evolutionary studies, and agricultural breeders have developed the third approach, which is particularly suited for outbred populations. Focusing first on the evolutionary approaches, inbred lines have been used to gain replicate phenotypic measures on genetically identical individuals reared within the same macro-environment, providing a combined estimate of the genetic variance in internal and external micro-environmental variance. As each genetic group has equal (zero) genetic variance (Falconer 1981), differences in the phenotypic variance among groups reflects genotype differences in total micro-environmental variance. Genetic variance in micro-environmental variance, measured in this way, has been demonstrated for a range of trait types in yeast, invertebrate, and plant models (Whitlock and Fowler 1999; Mackay and Lyman 2005; Hall *et al.* 2007; Ordas *et al.* 2008), with high broad-sense heritability estimates ranging from 36 to 81% (Ansel *et al.* 2008; Morgante *et al.* 2015). However, these broad-sense heritability estimates will be particularly sensitive to sample size, and cannot distinguish additive from nonadditive genetic contributions to micro-environmental variance. Data from recombinant inbred lines suggests that alleles underlying environmental sensitivity are largely recessive (Tonsor *et al.* 2013), consistent with the widely established decrease in environmental variance that occurs when inbred lines are crossed (Reeve and Robertson 1953; Lerner 1954; Falconer 1981), suggesting that nonadditive genetic variance may comprise a portion of the micro-environmental variance measured.

A second approach to estimating the genetic variance in micro-environmental variance, which has typically been used in evolutionary studies, has been to estimate levels of fluctuating asymmetry for replicated body parts: the within-individual component of micro-environmental variance (developmental instability). In contrast to the high broad-sense heritability of total micro-environmental variance, estimates of the narrow-sense heritability of fluctuating asymmetry (as a proxy for developmental instability) in outbred populations are typically low (Fuller and Houle 2002). However, fluctuating asymmetry is likely to be an incomplete measure of internal micro-environmental variation (Whitlock 1996, 1998; Leamy and Klingenberg 2005), and heritability might subsequently be much higher than inferred from these studies, comparable to heritability of life-history traits (Carter and Houle 2011). Therefore, it remains somewhat unclear whether narrow-sense heritability of micro-environmental variance is high or low in general. Investigations of developmental stability have also provided some evidence that nonadditive genetic effects partly determine levels of micro-environmental variance. There is some weak evidence that heterozygotes exhibit less fluctuating asymmetry than homozygotes (Leamy and Klingenberg 2005) and that developmental instability may be determined in part by epistasis (Pélabon *et al.* 2010). However, it is unknown whether nonadditive gene action might typically generate nonadditive genetic variance in micro-environmental variance.

If nonadditive genetic effects comprise a large component of the broad-sense heritability of micro-environmental variance, determining whether micro-environmental variance can respond to selection is difficult, because the response to selection is determined by the additive genetic covariance between traits and fitness (Robertson–Price identity). Under certain simplifying assumptions regarding the shape of the population fitness function, theoretical models predict that that when an intermediate phenotype has the highest fitness and a population resides near its fitness optimum, selection should reduce environmental (and genetic) variance, favoring individuals that exhibit a trait value near the fitness optimum (Gavrilets and Hastings 1994; Wagner *et al.* 1997; Mulder *et al.* 2007; Tonsor *et al.* 2013) (Figure 1B). In contrast, under episodes of fluctuating, disruptive, or directional selection, individuals with more extreme phenotypes are more likely to exhibit the trait value with the highest fitness, and therefore, the magnitude of micro-environmental variance is predicted to increase (Scharloo *et al.* 1967; Hill and Zhang 2004; Mulder *et al.* 2007; Tonsor *et al.* 2013) (Figure 1B). For both correlative studies in natural populations and manipulative artificial selection experiments, some studies have reported patterns that are consistent with the predicted relationships between fitness and the micro-environmental component of phenotypic variance, while others have found no relationship or one counter to predictions (reviewed in Lens *et al.* 2002; Leamy and Klingenberg 2005; Pélabon *et al.* 2010). Overall, the empirical and analytical difficulties associated with measuring micro-environmental variance (Whitlock and Fowler 1997; Clarke 1998; Lens *et al.* 2002), estimating its additive genetic variance (Fuller and Houle 2002; Hill 2007; Visscher and Posthuma 2010), and finally in measuring the selection on this component of phenotypic variance, have all hindered a clear consensus of how micro-environmental variance evolves in response to selection (Mulder *et al.* 2015) and therefore of the evolutionary significance of this variation (Lande 2009, 2015; Chevin *et al.* 2013).

Here, we adopt the third approach that has been used in animal breeding studies, which enables an estimate of the additive genetic variance in the total micro-environmental variance. We conducted a large, two-generation breeding design in an outbred population of *Drosophila serrata*, measuring a suite of wing-shape traits and fitness. In outbred populations, an estimate of the additive genetic variance in micro-environmental variance can be gained by treating the micro-environmental component of phenotypic variance as a trait itself, using a pedigreed population and the appropriate linear model (Hill 2007; Hill and Mulder 2010). First, the phenotypic variance in a trait is partitioned to the genetic (V_{G}) and environmental (V_{E}) components, and then V_{E} itself is partitioned into genetic and environmental components. In contrast to the high broad-sense heritability estimates of micro-environmental variance from inbred lines, estimates of additive genetic variance in micro-environmental variance from outbred agricultural populations suggest the narrow-sense heritability is very low. It has been found to range from 0.00 to 0.048 for traits such as litter size in sheep (SanCristobal-Gaudy *et al.* 2001), pigs (Sorensen and Waagepetersen 2003), and rabbits (Ibanez-Escriche *et al.* 2008); and body weight in snails (Ros 2004), broiler chickens (Mulder *et al.* 2009; Wolc *et al.* 2009), and rainbow trout (Janhunen *et al.* 2012). Whether these estimates from populations typically subject to generations of artificial directional selection are representative of randomly mating, outbred populations is generally unknown; although recent evidence from a wild great tit population similarly found a low narrow-sense heritability of 0.005 for micro-environmental variance of fledgling weight (Mulder *et al.* 2016). If the narrow-sense heritability is indeed very low, as these studies suggest, the evolutionary implications of micro-environmental variance may be substantially different than predicted from micro-environmental variance estimates from inbred lines.

Employing a two-step method of analysis (Mulder *et al.* 2009; Wolc *et al.* 2009), our experimental design enabled us to estimate the additive genetic variance in the micro-environmental variance of wing shape, in addition to estimating the genetic variance for a comprehensive measure of fitness. Our multi-generation breeding design enabled us to reduce the nonadditive genetic and Mendelian sampling variance that contributed to the phenotypic variance of V_{E}, improving on previous studies (Mulder *et al.* 2009, 2016). Finally, by measuring both fitness and micro-environmental variance in the same breeding design we were able to directly test, for a population predicted to be under mutation–selection balance, whether micro-environmental variance could respond to selection (Mulder *et al.* 2015).

## Methods

The experimental design, based on a large, laboratory-adapted population of *D. serrata* housed under standard conditions at 25° (Hine *et al.* 2014), is described in detail in Sztepanacz and Blows (2015). Briefly, the first generation employed a paternal-half-sibling breeding design, where 75 sires were each mated to three virgin dams (total of 225 dams). The second generation employed a double-first-cousin breeding design, where two pairs of full brothers from each dam (*i.e.*, four sons per dam) were mated to two pairs of full sisters from a different sire, such that the two pairs of sisters were not related to each other or to the mating partners of their brothers and sisters. Upon emergence, virgin sons were collected from each of the 685 families that produced offspring in that generation. A median of five sons per family were allocated for use in a competitive fitness assay and a median of eight sons per family (brothers of the sons allocated for the fitness assay) were allocated for wing phenotyping.

This experimental design enabled us to estimate the genetic variance in the mean of wing traits, in addition to the genetic variance in the micro-environmental variance of these same traits. Because we also measured fitness on the brothers of individuals that had their wings phenotyped, we were able to determine how both wing shape and its micro-environmental variance genetically covaried with fitness, and therefore whether wing shape and its micro-environmental variance could respond to selection. Under our experimental design, we consider variance among full siblings reared within the same vial. Variation in larval density among vials could potentially affect the phenotypic variance within families, and if density is largely determined by genetic variance in fitness, then this genetic variance in fitness can have a causal effect on phenotypic variance in wing phenotypes. A previous analysis of these wing data provides no evidence that common environment is an important component of the genetic variance in wings (Sztepanacz and Blows 2015), but the effect of common environment on micro-environmental variance is unknown. In total, we analyzed one wing from each of 5040 individuals in 685 families, and fitness from 2883 individuals in 666 families (Table 1).

Male wing phenotypes were measured for up to 15 males per full-sib family, with a median of eight wings per family. One wing from each male (either left or right) was removed and mounted on a standard microscope slide using double-sided tape. Wings were photographed and nine landmarks (Figure 2) corresponding to those previously described (McGuigan and Blows 2007, 2013) were recorded using tpsDig2 software (Rohlf 2005a,b). Landmarks were then aligned by generalized Procrustes least squares superimposition using tpsRegW (Rohlf and Slice 1990; Rohlf 1999, 2005a,b). Eight interlandmark distances were then calculated and used in subsequent analyses: aLM2–4, aLM2–5, aLM2–8, aLM3–7, aLM3–9, aLM4–5, aLM5–6, and aLM5–8 (Figure 2; Sztepanacz and Blows 2015). For simplicity, we refer to these traits as wing shape, however, because we did not fit size as a covariate in our analyses, the interlandmark distances are confounded with nonuniform allometric variation and therefore do not exclusively describe shape. By measuring only one wing from each male, we were unable to partition the within-individual variance in wings (fluctuating asymmetry) from the among-individual variance, and therefore these data result in a combined estimate of the internal and external micro-environmental variance in wings. As wings were chosen haphazardly from either the left or the right, the phenotypic variance in wing shape may also be confounded with directional asymmetry. While there is no evidence for directional asymmetry in *D. serrata*, it has been demonstrated in *D. melanogaster* (Klingenberg *et al.* 1998), and if present may therefore inflate the phenotypic variance estimated for our wing traits.

Male fitness was assessed for up to five sons from each family. A focal male from the breeding design (red eyed) was placed in a vial with two females and a competitor male (orange eyed), where flies were allowed to freely mate and oviposit for 48 hr before being discarded from the vials. After eclosion, the resulting offspring were scored for paternity using eye color, where orange is recessive to red. In this assay, 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. Vials that did not produce any offspring (zeros), likely due to female damage or death, were removed prior to analyses. In the following analyses, male competitive fitness was calculated as the natural logarithm of the odds ratio of the number of focal-male- to competitor-male-sired offspring (Reddiex *et al.* 2013; Collet and Blows 2014):Prior to analysis, a total of 24 multivariate outliers (0.5% of the total wing data, distributed randomly across families) were identified and removed from the wing-phenotype data using the multivariate Mahalanobis distance technique implemented in the software package JMP (version 10; SAS Institute, Cary, NC). The fitness and wing data were approximately normally distributed (summary statistics Table 1) and were not standardized; however, the data were centered on their trait-specific means and interlandmark distance values were rescaled (all multiplied by 1000) to ensure that the residuals, which were used as traits in part of the analysis, would not hinder model convergence due to their small scale.

### Genetic analyses

#### General overview:

The genetic analysis to estimate the additive genetic variance in the means of wing traits and in their micro-environmental variance involved a two-step approach that has previously been employed in several studies (*e.g.*, Mulder *et al.* 2009; Wolc *et al.* 2009; Janhunen *et al.* 2012). The first step employed a standard quantitative genetic model to estimate the additive genetic, dominance genetic, and residual variance for the means of the wing-shape traits () and for the mean of fitness (*w*). The second step used the residuals of the wing-shape traits obtained from the first model as the observations (phenotypic values) for the subsequent models, after transforming these residuals to represent the natural logarithm of squared deviations of individuals from their family mean Expressed in this way, these data enabled us to estimate the additive genetic variance in the micro-environmental variance of wing-shape traits. In this two-step approach, the transformed residuals that are used as phenotypic observations in the second-step models are a combination of the true micro-environmental effects of interest, and estimation error, which is not accounted for in the second-step model. The accuracy of the squared residuals is reflected in the diagonal elements of the “hat matrix,” which describe the leverages (influence) of observed values on the fitted values for the same observations (Hill and Mulder 2010). To determine whether any of our observations were particularly influential, we examined the leverages for every observation of each trait. The leverages were small with little variation for any of the traits studied (Supplemental Material, Figure S1 in File S1) and, therefore, accounting for the leverages by employing a double hierarchical generalized linear model (DHGLM) of analysis (Hill and Mulder 2010; Rönnegård *et al.* 2010; Mulder *et al.* 2013, 2016) would not be expected to have a large impact on our results. However, in situations where leverages are substantial, a DHGLM approach may be more suitable.

To determine whether wing shape and its micro-environmental variance were subject to selection, we used the measures of fitness that were obtained on the brothers of individuals that had their wings measured. First, we determined whether wing shape itself was under quadratic selection by estimating the additive genetic covariance between fitness (*w*) and the squared deviations of each wing-shape trait from their population mean (). Second, to establish whether the micro-environmental variance of wing shape was also subject to selection, we determined the additive genetic covariance between fitness (*w*) and the micro-environmental variance traits (). As described above, micro-environmental variance traits were, by definition, squared deviations from family means. Therefore, the genetic covariance between fitness and micro-environmental variance traits directly indicated whether these traits were also subject to quadratic selection.

#### Step 1: genetic variance in the mean of wing-shape traits:

All genetic analyses were conducted using restricted maximum likelihood within an animal model framework, implemented in Wombat (Meyer 2007), where the numerator relationship matrix was used to estimate the additive genetic values and the dominance relationship matrix was used to estimate the dominance deviations. The dominance genetic relationship matrix was calculated and inverted using the nadiv package (Wolak 2012) in R (R Core Team 2012) before being supplied to Wombat. Univariate quantitative genetic analyses of individual wing-shape traits and fitness employed the following mixed model:(1)where *y _{i}* is either fitness or a wing-shape trait for individual

*i*. The numerator relationship matrix was used to estimate

*a*, the additive genetic value of individual

_{i}*i*(population mean of 0, variance of ); the dominance genetic relationship matrix was used to estimate

*d*, the dominance genetic value of individual

_{i}*i*(mean 0, variance); and

*e*was the residual error (mean 0, variance). Initially, a common environment effect was included in the model to account for full siblings being reared in the same vial. However, common environment explained <1% of the variance in wing shape and did not differ significantly from zero for any of the traits analyzed, therefore this effect was excluded from these analyses (Sztepanacz and Blows 2015). Statistical support for the additive and dominance genetic variance in each trait was assessed by using a log-likelihood ratio test with one degree of freedom, comparing the fit of a model that included

_{i}*vs.*excluded the variance component of interest. These models were identical to the univariate models presented in Sztepanacz and Blows (2015), however, in the current analysis the data have not been standardized to unit variance. Residuals, expressed as deviations from the population mean, were obtained for each individual in the terminal generation of the pedigree. In addition to the univariate analyses of each interlandmark distance trait, we also estimated the additive genetic covariance matrix of all eight interlandmark distance traits,

**G**, as detailed in Sztepanacz and Blows (2015).

#### Step 2: genetic variance in the micro-environmental variance of wing-shape traits:

The residuals obtained from (1) quantify each individual’s deviation from the population mean phenotype, after accounting for the deviations resulting from additive and dominance genetic effects. Therefore, these data capture the external micro-environmental variance plus internal micro-environmental variance specific to each individual. To estimate the genetic variance in micro-environmental variation, we first transformed these residual data. Because we focus on analyzing the phenotypic variance generated by stochastic micro-environmental variation experienced among families, we first centered individual residual deviations on their family mean. Second, we square these residuals because we are interested in deviations about the mean but not in which direction they might occur. Third, we applied the natural logarithm to these residuals after they had been recentered and squared. Means and variances can be positively correlated, such that families with large trait means also tend to have large variances (Walsh and Lynch 2010). Taking the natural logarithm reduces the potential for such scaling effects to generate spurious results (Mulder *et al.* 2009). The phenotypes for micro-environmental variance used in this step of the analysis were therefore for each wing-shape trait, with these data describing the magnitude of an individual’s deviation from its family mean.

To determine the additive genetic variance in the micro-environmental variance of individual wing-shape traits, univariate analyses employed the following mixed model:(2)whereis a micro-environmental variance trait for individual *i* on the ln scale [we simplify the notation of to just from this point forward to reduce the complexity of the subscripting below]. The numerator relationship matrix was used to estimate the additive genetic value of individual *i* for (population mean of 0, variance and is the residual variance of individual *i* for(mean 0, variance which accounts for random variation that is not explained by genetic effects on either mean shape or micro-environmental variance in shape. Initially, a dominance effect was also included in these models; however, the results indicated that dominance variance was not an important component of the genetic variance in micro-environmental variance (Table S1 in File S1). We therefore excluded the dominance genetic effect from our models.

To compare the levels of additive genetic variance for the mean of individual wing-shape traits that were obtained from (1) () to the micro-environmental variance of these traits that were obtained from (2) the variance components must be expressed as a proportion of the total phenotypic variance in wing shape. To do this, we used the heritability, as this metric has commonly been reported from other studies of genetic variance in trait means and micro-environmental variances. Heritability can be calculated as a proportion of using the parameter estimates derived from (2): To express the heritability estimates as a proportion of the phenotypic variance in the eight interlandmark distance traits that were measured, and thus as the same proportion as (1), we used the equations derived by Mulder *et al.* (2007). First, we calculated the additive genetic variance of micro-environmental sensitivity on the measured trait scale as: where is calculated as above and comes from (1). We then calculated the heritability of micro-environmental variance as a proportion of the phenotypic variance in the measured trait:(3)where is the phenotypic variance from (1).

#### Multivariate analysis of micro-environmental variance:

In addition to single-trait models, we carried out a multivariate genetic analysis of micro-environmental variance for all eight wing-shape traits. This analysis allowed us to estimate the additive genetic covariance matrix for micro-environmental variance in wing shape The multivariate linear model was:(4)where the phenotypic observations, for all individuals came from the univariate model (1) and were analyzed with a model where **Z** was a design matrix relating observations to the additive genetic effects **a**. **I** was an identity matrix, and **e** a vector of residual errors. The random effect (and residuals) were assumed to be normally distributed and elements of **a** were further assumed to be drawn from **a∼***N*(0, **A**), where was the additive genetic covariance matrix for micro-environmental variance in wing shape and **A** was the numerator relationship matrix. The statistical dimensionality of the covariance matrix, was evaluated using log-likelihood ratio tests for a series of nested reduced rank (factor analytic) models (Hine and Blows 2006).

#### Stabilizing selection on the mean of wing shape:

The form of selection acting on the mean of wing shape will have consequences for whether the genetic covariance between fitness and the micro-environmental variance of wing shape is predicted to be positive or negative. For a trait whose mean is under directional selection, the genetic covariance between fitness and the micro-environmental variance is predicted to be positive (variance-increasing selection), whereas under stabilizing selection a negative genetic covariance (variance-reducing selection) is predicted. Therefore it was important to first establish the form of selection acting on the mean of wing shape.

Previous analyses of these data indicated that wing-shape traits were subject to directional selection (Sztepanacz and Blows 2015), but whether wing shape was also subject to nonlinear selection was not explored. To quantify the multivariate nonlinear selection acting on the genetic variance in the mean of wing-shape traits, we estimated the additive genetic covariance matrix between fitness (*w*) and the products of the eight measured wing-shape traits when expressed as phenotypic deviations from their population mean (**φ**) (Rausher 1992; Delcourt *et al.* 2012). Because wing-shape traits were centered on their trait-specific means prior to estimating the additive genetic variance in their mean, they were already expressed as deviations from the population mean. Therefore, we simply needed to generate all cross-products prior to parameterization of **φ**. The parameterization of **φ** (Table S2 in File S1) was achieved by estimating each element of the matrix as the genetic covariance between fitness (*w*) and the corresponding trait deviation cross-products using bivariate models. A total of 36 bivariate models were run to estimate all elements of **φ**, with the standard format of the model shown here for one cross-product and fitness:where fitness and product of the *j*th and *k*th trait deviations were treated together as dependent variables, and were analyzed with a model where **Z _{1}** was a design matrix relating observations to the additive genetic effects

**a**,

**Z**was a design matrix relating observations to the dominance genetic effect

_{2}**d**,

**I**was an identity matrix, and

**e**a vector of residual errors. The random effects (and residuals) were assumed to be normally distributed and elements of

**a**were further assumed to be drawn from

**a∼**

*N*(0,

**G**

**A**), where

**G**was the additive genetic covariance matrix and

**A**was the numerator relationship matrix. Similarly, the elements of

**d**were assumed to be drawn from

**d∼**

*N*(0,

**D**

**Z**), where

**D**was the dominance genetic covariance matrix and

**Z**was the dominance genetic relationship matrix. The residual covariance between fitness and the products of trait deviations was fixed at zero to reflect these traits being measured on different individuals. Diagonalization of

**φ**yielded the eigenvectors and corresponding eigenvalues that characterize the form of nonlinear selection acting on wing shape, with negative eigenvalues being indicative of stabilizing selection (Delcourt

*et al.*2012).

#### Selection on micro-environmental variance of wing shape:

Once we had determined the form of selection acting on the mean of wing shape, the next step was to determine how the genetic variance in the micro-environmental variance of wing shape was associated with fitness. First, we focused on two multivariate trait combinations that were associated with statistically significant genetic variance, as determined by the factor analytic rank test (described above and in *Results*) of the eight-trait micro-environmental variance additive genetic covariance matrix, This enabled us to place statistical significance on the additive genetic and dominance genetic covariances between fitness and the major axes of micro-environmental variance, which captured 67% of the total genetic variance in micro-environmental variance. We scored all individuals for these two multivariate trait combinations using**Z**, where **Z** is a row vector of the observed micro-environmental variance traits for an individual and is either or (the first or second eigenvector of ). We then fit a bivariate model [of the same form as model (4)] to each of these composite traits with fitness to estimate both the additive genetic and the total (additive plus dominance) genetic covariance between them. The confidence intervals of these covariance estimates were obtained by employing the “sample” procedure in Wombat (Meyer 2007), which samples from the distribution given by inverse of the Fisher Information matrix, to approximate the standard error on (co)variance components and their functions (Meyer and Houle 2013; Houle and Meyer 2015). The elements of which describe the (co)variance among micro-environmental variance traits, are themselves squared deviations. Therefore, a negative covariance between fitness and an individual's score for the multivariate trait combinations is indicative of variance-reducing selection (Delcourt *et al.* 2012; Mulder *et al.* 2015).

In addition to focusing on the two composite traits above, which each explained a significant portion of the genetic variance in micro-environmental variance, we employed an alternate approach that enabled us to determine the total genetic covariance between all eight of the micro-environmental variance traits (capturing 100% of the genetic variance in micro-environmental variance) and fitness. Due to parameter constraints, we were unable to fit all nine traits simultaneously (eight wing-shape micro-environmental deviations plus fitness). Therefore, we fit bivariate models [as (4) above] between each pair of traits; these models simultaneously estimated additive genetic, dominance genetic, and residual (co)variances. We then used a maximum-likelihood approach, implemented in Wombat (Meyer 2007, 2012), to pool the estimates from bivariate analyses, generating full (nine trait) additive and dominance covariance matrices. The justification behind this pooling approach is given in Sztepanacz and Blows (2015), with details of the method provided in Meyer (2012). Briefly, this approach simultaneously pools the bivariate additive estimates to form a nine-trait additive covariance matrix and pools the bivariate dominance estimates to form a nine-trait dominance covariance matrix, while accounting for the typically strong negative sampling correlations between estimates from analyses with overlapping subsets of traits (Meyer 2012). After pooling, to obtain the matrix of total genetic (co)variance, we added together the additive and dominance covariance matrices. In the resulting matrix of total genetic (co)variance, the matrix elements of the column (or row) corresponding to fitness (excluding the diagonal element that represented the genetic variance in fitness) indicate the genetic covariance between each micro-environmental variance trait and fitness. These elements, when arranged as a vector, are analogous to the genetic selection differentials ** s_{g}**, or when scaled by the genetic selection gradient

**(Stinchcombe**

*β*_{g}*et al.*2013). However, as micro-environmental variance traits are squared deviations, this vector describes the predicted responses to quadratic, and not directional, selection.

### Data availability

File S2 contains the pedigree and corresponding phenotype information for each of the eight wing-shape traits and the log-odds fitness trait that were analyzed here. Descriptions of each column in the data are provided at the top of the file.

## Results

### Genetic variance in the mean of wing-shape traits

The heritabilities for the mean of individual wing-shape traits ranged from 0.296 to 0.564, and the heritability of fitness was three orders of magnitude smaller at 0.233 × 10^{−3} (Table 2). We found no statistical support for the presence of additive genetic variance in fitness, but there was for each wing-shape trait (Table 2). We also found statistical support for dominance variance in three of eight wing traits and for fitness (Table 2).

### Genetic variance in the micro-environmental variance of wing-shape traits

Univariate analyses of the micro-environmental variance in each wing-shape trait demonstrated statistically significant additive genetic variance in for five of the eight micro-environmental variance traits (Table 3). The heritabilities of the micro-environmental variance of wing-shape traits were <1% in all cases, and differed markedly from the heritabilities for the mean of wing-shape traits, being on average two orders of magnitude smaller (Table 2 and Table 3). The multivariate analysis of micro-environmental variance traits revealed strong correlations among many of them, with positive covariances predominating (21 of 29 covariances in were positive; Table 4). Factor analytic modeling revealed statistical support for only the first two eigenvectors of accounting for 67% of the total genetic variance (reducing from two to one dimensions significantly worsened the fit of the model *χ*^{2} = 15.43, d.f. = 7, *P* = 0.03). The contrast between the total additive genetic variance explained by the first two eigenvalues of (67%) and the summed variance of the two individual traits with the largest variance (38%) highlights the extent to which the micro-environmental variance wing-shape traits shared a common genetic basis.

The genetic covariance between mean wing shape and its micro-environmental variance will determine whether the response to selection acting on micro-environmental variance is constrained by selection acting on mean wing shape. However, estimating the genetic covariance between the mean of a trait and its micro-environmental variance has proven challenging for traits such as litter size in pigs and rabbits (Hill and Mulder 2010; Yang *et al.* 2011). In a preliminary analysis of a single wing-shape trait (aLM2–4), the parameter estimate for the covariance between the mean and micro-environmental variance of this trait was 0.248 and statistically significant (*χ*^{2} = 16.96, d.f. = 1, *P* < 0.001). However, when micro-environmental variance phenotypes were randomized across the pedigree, we also found evidence for a statistically significant covariance between mean wing-shape and its micro-environmental variance (parameter estimate = −0.13, *χ*^{2} = 7.016, d.f. = 1, *P* = 0.008), suggesting the possibility of a spurious covariance in the unrandomized data and precluding us from obtaining a quantitative estimate of this evolutionarily important covariance. As a proxy, we determined whether trait combinations associated with high additive genetic variance for micro-environmental variance in wing-shapes were also associated with high additive genetic variance for mean shape by projecting each of the first two eigenvectors of through **G** (Figure 3). Projecting vectors, including selection gradients, through **G** is a common approach for determining the genetic variation available for a response to selection in that direction of trait space (Blows *et al.* 2004; Hansen and Houle 2008). Here, the projection demonstrated that the first two eigenvectors of were each associated with relatively high levels of additive genetic variance: 18 and 11% of the total additive genetic variance in wing-shape trait means, respectively (Figure 3). Although this result does not provide any definitive proof that the same genetic variation contributes to trait mean and to trait variance, it does indicate that wing-shapes associated with high levels of additive genetic variance (and high evolutionary potential) for micro-environmental variance are also associated with high levels of additive genetic variance (evolutionary potential) in the mean. Such a relationship suggests the potential for selection on wing trait means to bias any response to selection on wing-shape micro-environmental variance.

### Stabilizing selection on the mean of wing-shape

Consistent with previous evidence in *Drosophila* (McGuigan *et al.* 2011; McGuigan and Blows 2013), we found a considerable amount of the selection on wing shape to be stabilizing. An eigenanalysis of the genetic covariance matrix between fitness and the products of the eight wing traits (expressed as deviations from their population means) (**φ**) returned four positive and four negative eigenvalues (Table S2 in File S1). The negative eigenvalues, indicative of stabilizing selection, accounted for 61% of the total quadratic selection present; therefore, there was evidence for stabilizing selection because the negative eigenvalues explained >50% of the total quadratic selection. Because **φ** was composed from bivariate models, it was not possible to determine statistical support for the genetic dimensions of this matrix within the framework of the multivariate mixed model employed here.

### Selection on the micro-environmental variance of wing-shape

In addition to the stabilizing selection on wing-shape, the additive genetic effects (breeding values) for the micro-environmental variance of wing-shape were correlated with those for fitness. There was a negative additive genetic covariance between fitness and each of the first two axes of additive genetic variance in the micro-environmental variance of wing-shape of −0.064 and −0.069, for and respectively. A negative correlation is indicative of variance-reducing selection (Delcourt *et al.* 2012; Mulder *et al.* 2015). However, due to the low additive genetic variance in fitness, the SE of these covariances was relatively large, such that the confidence intervals of these estimates overlapped zero [ −0.19, 0.07; −0.18, 0.04]. Therefore, while the additive genetic covariance between fitness and each of the major axes of micro-environmental variance suggests weak variance-reducing selection, the evidence for this association was marginal. The projection of the normalized vector, whose elements were the additive genetic covariance between each micro-environmental variance trait and fitness (Table 5) through also supported only a very weak association between the breeding values for fitness and those for micro-environmental variance (Figure 4).

While only the presence of an additive genetic covariance between fitness and a trait provides evidence that the trait can respond to selection, a genetic covariance between fitness and the nonadditive genetic variance that contribute to trait variation suggests that there are alleles affecting both the trait and fitness, and consequently that they share a genetic basis. Here, the total (additive plus dominance) genetic variance in fitness was positively correlated with micro-environmental variance. The genetic correlation between fitness (whose genetic variance was primarily comprised of dominance variance) and the first two eigenvectors of (whose genetic variance was comprised of additive variance) was 0.049 [−0.131, 0.229] and 0.179 [−0.031, 0.389], for and respectively. It is possible that the shared genetic basis between fitness and micro-environmental variance in wing-shape does not fall along either of these two orthogonal axes () of micro-environmental genetic variance. To explore this, we projected the normalized vector whose elements were comprised of the covariance between the total genetic variance (additive plus dominance) in fitness and in each micro-environmental variance trait (Table 5) through The projection indicated that 29% of the genetic variance in was associated with fitness (Figure 4). Accordingly, this vector fell between the first and second eigenvalues of within a subspace that was found to have statistically significant genetic variance (Figure 4). Therefore, we can infer an association between the genetic variance in micro-environmental variance and in fitness; however, this association seems to be driven primarily by dominance variance in fitness.

## Discussion

Despite considerable interest over the past seven decades (*e.g.*, Reeve and Robertson 1953; Waddington 1953; Lerner 1954; Schlichting 1986; Gavrilets and Hastings 1994; Wagner *et al.* 1997; Mulder *et al.* 2007, 2015, 2016; Tonsor *et al.* 2013), the evolution of micro-environmental variance remains unresolved. Two major empirical challenges have slowed progress on assessing the relative evolutionary significance of this variation, namely the difficulty of estimating the additive genetic variance in micro-environmental variance in randomly mating, outbred populations, and in characterizing the selection acting on this variation. Here we demonstrated that the micro-environmental variance of a suite of wing-shape phenotypes in a randomly mating, outbred population of *D. serrata* had a genetic basis that was shared with fitness, and therefore that micro-environmental variance may respond to selection.

Univariate analyses uncovered significant additive genetic variance in the micro-environmental variance for five of the eight wing-shape traits. However, the heritabilities of the micro-environmental variance traits (as a proportion of the total phenotypic variance in the measured wing-shape traits) were low (<1%). Heritabilities of this magnitude are a marked departure from standard quantitative trait heritabilities, which are typically between 20 and 60% (Lynch and Walsh 1998), and between 30 and 57% for the wing-shape traits in this study. Although the average reported heritability of micro-environmental variance in agricultural populations is higher than found here (3.8%; Hill and Mulder 2010), a similarly large difference (an order of magnitude) in heritability between a trait mean and its micro-environmental variance is typical of agricultural studies (*e.g.*, Mulder *et al.* 2009; Wolc *et al.* 2009, 2011; Janhunen *et al.* 2012).

The low narrow-sense heritability of micro-environmental variance in *D. serrata* wings and for life history or production traits in agricultural studies is in contrast to the broad-sense heritability of micro-environmental variance that has been estimated from inbred lines. For traits such as abdominal bristle number (Mackay and Lyman 2005), sternopleural bristle number, chill coma recovery, startle response, and starvation stress resistance in *D. melanogaster* (Morgante *et al.* 2015), and gene expression in yeast (Ansel *et al.* 2008), heritabilities of micro-environmental variance were an order of magnitude larger than reported here, in the 36–81% range. Inbred line experiments may result in overestimates of the heritability; inbred lines have double the additive genetic variance compared to outbred populations, with estimates also confounded by dominance variance (Falconer 1981). For quantitative traits in general, dominance is predicted to be an important mechanism for masking the effects of deleterious alleles, and therefore for maintaining invariant phenotypes (Meiklejohn and Hartl 2002; Bagheri-Chaichian *et al.* 2003; Bagheri and Wagner 2004; Agrawal and Whitlock 2011). Although we found little evidence that dominance variance contributed to the levels of genetic variance in micro-environmental variance, our power to detect dominance genetic variance in micro-environmental variance was low. Additionally, a lack of dominance variance does not indicate a lack of dominant gene action (Hill *et al.* 2008; Huang and Mackay 2016), it simply demonstrates that the variance in dominance deviations for all loci affecting the trait of interest are, on average, zero.

In general, our results suggest that the heritability of micro-environmental variance might typically be lower in natural populations than is implied by estimates from inbred lines, an observation that, if general across other types of traits and populations, suggests that the response to selection on environmentally generated trait variability might be limited. Carter and Houle (2011) applied artificial directional selection to fluctuating asymmetry of *D. melanogaster* wings, estimating the realized heritability of this component of micro-environmental variance. Comparable to our heritability estimates of total micro-environmental variance, the heritability of fluctuating asymmetry was <1%. Similarly, in a ferent experiment the lack of selection response of fluctuating asymmetry after 10–14 generations of artificial selection also suggests that heritability is low (Tsujino and Takahashi 2014). Inconsistent selection responses for traits associated with very low additive genetic variance suggests that rare alleles contribute the variance, and consequently the random sampling of rare alleles in small populations could explain the lack of evolutionary response (Hine *et al.* 2014). Overall these studies of outbred populations are consistent with our results, suggesting that genetic variance for micro-environmental variance is likely much lower than suggested by analyses of inbred lines.

The first eigenvector of the eight-trait additive genetic covariance matrix for the micro-environmental variance traits () accounted for 48% of the genetic variance present in micro-environmental variance, more than double the amount accounted for by any single trait, indicating that the covariance structure among micro-environmental variance traits concentrates genetic variance onto specific axes. All but one trait contributed concordantly to (Table 5), reflecting the predominantly positive additive genetic covariance among micro-environmental deviations (Table 4). The mechanisms buffering phenotypes against stochastic micro-environmental variation are poorly known, and it remains unclear whether the levels of micro-environmental variance are determined by trait-specific or common-buffering mechanisms (Breuker *et al.* 2006; Siegal and Leu 2014). Our results are consistent with a common-buffering mechanism shared among wing-shape traits, conferring robustness against micro-environmental perturbations. However, considering that all of the traits here are aspects of wing shape, which likely share a developmental pathway, their sensitivity to environmental perturbations may be expected to be positively correlated (Leamy and Klingenberg 2005). Nonetheless, our results highlight the potential of multivariate quantitative genetic analyses to determine the extent to which different types of traits share a common micro-environmental buffering mechanism.

Although there has been considerable theoretical attention given to the effects of selection on micro-environmental variance (Zhang 2005; Zhang and Hill 2005, 2007), there is currently limited empirical data addressing whether micro-environmental variance is under selection (Mulder *et al.* 2015, 2016). Carter and Houle (2011) observed greater response to artificial selection to increase fluctuating asymmetry than to decrease it. Such a biased response to selection is typical of fitness traits, suggesting that genetic variance is due to rare alleles kept at low frequency due to their deleterious effects on fitness (Frankham 1990; Hill and Caballero 1992; Falconer and Mackay 1996; McGuigan and Blows 2007; Sztepanacz and Rundle 2012). In general, understanding the evolutionary consequence of selection requires a genetic approach. The Robertson–Price identity predicts the response to selection from the genetic covariance between fitness and the trait of interest. This approach overcomes the challenges to inference posed both by environmentally generated covariances between trait and fitness (Rausher 1992; Kruuk *et al.* 2002; Stinchcombe *et al.* 2002), and by the “missing fraction,” individuals removed from the population by viability selection prior to their adult phenotypic traits being measured (Hadfield 2008; Mojica and Kelly 2010; Polak and Tomkins 2013). Here, our results provide evidence that fitness and the environmental variance in wings in this population of *D. serrata* share a genetic basis that is suggestive of variance-reducing selection operating on micro-environmental variance.

In our experimental population, which had been kept under constant conditions in the laboratory for ∼50 generations prior to the experiment, we found significant (and moderately high) dominance variance, but no evidence for additive genetic variance in fitness (Table 2) (Sztepanacz and Blows 2015). This suggests that persistent directional selection on fitness had eroded much of its additive genetic variance and, therefore, that this population may be evolving under variance-reducing (stabilizing) selection in the vicinity of its phenotypic optimum. Our results were consistent with stabilizing selection acting on the mean of wing-shape phenotypes to some extent, with the multivariate genetic analysis of fitness and squared phenotypic deviations of each wing-shape trait from the population mean (**φ**) (Table S2 in File S1) indicating that 61% of the quadratic selection present on wings was stabilizing. However, estimating **φ** from a series of bivariate models precluded us from testing how much of this selection was statistically significant.

Given the presence of stabilizing selection on wing shape, we would predict the micro-environmental variance of wing shape to also be subject to variance-reducing selection (Zhang and Hill 2008; Hill and Mulder 2010). The two-dimensional subspace of that was found to have statistically significant genetic variance in micro-environmental variance was found to consistently display negative genetic correlations with the additive genetic variance in fitness, indicative of variance-reducing selection. However, in both cases the upper confidence intervals of these genetic correlations overlapped zero, limiting our ability to draw a strong conclusion from these associations. In contrast, the genetic correlations between dominance variance in fitness and the additive variance in the same subspace of the micro-environmental variance were consistently positive, indicating that alleles generating dominance deviations with respect to fitness may increase the levels of micro-environmental variance. Although we did not find evidence that dominance variance was an important component of micro-environmental variance, alleles resulting in dominance deviations with respect to fitness may act additively for micro-environmental variance (Lynch and Walsh 1998). Therefore, fitness variation that is presumably generated by recessive alleles for a population at mutation–selection balance may also increase levels of micro-environmental variance. This result is consistent with data from inbred line crosses and recombinant inbred lines in *Arabidopsis thaliana* which indicates the alleles underlying environmental sensitivity are largely recessive (Hall *et al.* 2007).

In summary, our results provide evidence for an additive genetic basis of micro-environmental variance for a suite of morphological traits that describe wing-shape variation in a randomly mating, outbred population of *D. serrata*. While there was limited evidence for variance-reducing selection on the micro-environmental variance in this population, there was a stronger association between dominance genetic variance in fitness and the micro-environmental variance. Consistent with previous studies of developmental stability and environmental canalization, our data suggest that dominance genetic effects may mask much of the genetic variance in sensitivity to environmental perturbations. Additional empirical studies of this nature will be needed to assess how the genetic basis of fitness is associated with the micro-environmental variance of quantitative traits, and to determine how trait means and their micro-environmental variance evolve under variance-reducing selection.

## Acknowledgments

We thank Miranda Vidgen, Donna Petfield, Amanda Denton, Grace Stanton, Julie Collet, Emma Hine, and Jessica McCarroll for help in the laboratory.

## Footnotes

Supplemental material is available online at www.genetics.org/lookup/suppl/doi:10.1534/genetics.116.199075/-/DC1.

*Communicating editor: C. D. Jones*

- Received December 11, 2016.
- Accepted June 12, 2017.

- Copyright © 2017 by the Genetics Society of America