One Trait Case

Griffin et al.’s (2013) analysis focused on $$ r_{MF}$$ as an indicator of the degree of constraint. Consider the simple case of SAS on the expression of one gene. The predicted response to selection is:

The change in sexual dimorphism under antagonistic selection is:

The effect of $$ r_{MF}$$ can be highlighted by substituting $$ b=r_{MF}\sqrt{mf}$$, which yields,

Under concordant selection, the change in dimorphism is:

This can be put in the same form as equation (3), yielding:

Thus, concordant selection will change dimorphism whenever there are different genetic variances in the two sexes, as previously noted by many authors (see the Introduction). The response can either increase or decrease sexual dimorphism if selection is on a previously dimorphic trait.

Comparison of equations (3) and (5) shows that in the single trait case, the rate of change in dimorphism will be higher under concordant selection than antagonistic selection of equal strength, when $$ d<(2r_{MF})/(r_{MF}^2+1)$$. This condition is increasingly easy to satisfy the larger $$ r_{MF}$$is. At the average value $$ r_{MF}=0.75$$ found in Poissant et al.’s (2010) review of intersexual correlations, the condition is met when the ratio of the larger of m and f to the smaller is 1.5 or larger; when $$ r_{MF}=0.9$$ a ratio of 1.16 is enough for concordant selection to dominate changes in dimorphism.

Two Trait Case

When k = 2, differences in genetic variances between the sexes remain a key source of dimorphism. However, to focus on other features of the G matrix that can promote dimorphism, but are absent from the k = 1 case, we consider the special case when all trait variances are 1, yielding the G matrix:

If only the focal trait is under selection (e.g., for antagonistic selection $$ {\mathbf{\beta }}_A^T=[\begin{array}{cccc}\beta & 0& -\beta & 0\end{array}]$$), the above results hold for the selected trait. However, there is also an indirect response, which leads to a change in dimorphism in $$ z_2$$ under antagonistic selection on $$ z_1$$ of:

The first term in brackets is the average correlation between the traits in males and females, and the second is the average cross-sex correlation. Under concordant selection, sexual dimorphism changes at the rate:

The fact that there are indirect responses to selection makes the interpretation of the existing degree of dimorphism in particular traits more challenging. For example, it is quite possible for the change in dimorphism of the selected trait to be less than that of the unselected trait. The direct response to antagonistic selection, $$ {\Delta }_{A.1}$$, will be small when $$ r_{MF1}\approx 1$$, making it plausible that $$ {\Delta }_{A.2}>{\Delta }_{A.1}$$ if within-sex, cross trait correlations$$ r_{M1M2} \text{and} r_{F1F2}$$ are larger than cross-sex cross-trait correlations $$ r_{M1F2} \text{and} r_{M2F1}$$.

If both traits are under directional selection, there are two orthogonal, SAS vectors $$ {\mathbf{\beta }}_{A1}^T=\left[\begin{array}{cccc}\mathit{\beta }& \mathit{\beta }& \mathit{-}\mathit{\beta }& \mathit{-}\mathit{\beta }\end{array}\right]$$ and $$ {\mathbf{\beta }}_{A2}^T=\left[\begin{array}{cccc}\mathit{\beta }& \mathit{-}\mathit{\beta }& \mathit{\beta }& \mathit{-}\mathit{\beta }\end{array}\right]$$. Vector $$ {\mathbf{\beta }}_{A1}$$, for example, predicts a response vector:

The direct responses are given by the leftmost term in parentheses. The middle term gives the indirect responses due to the average within-sex correlations between traits, $$ {\overline{r}}_w$$. The final term gives the indirect responses due to the average off-diagonal (between-sex) correlations in B,$$ {\overline{r}}_b$$. The complementary antagonistic selection gradient, $$ {\mathbf{\beta }}_{A2}$$, reverses the signs of the indirect effects, so that $$ {\overline{r}}_w$$retards the evolution of dimorphism, whereas $$ {\overline{r}}_b$$ facilitates it. Thus, the effects of $$ {\overline{r}}_w$$ and $$ {\overline{r}}_b$$ are always antagonistic to each other. Dimorphism is always constrained by $$ r_{MF}$$.

Concordant selection, $$ {\mathbf{\beta }}_{C1}^T=[\begin{array}{cccc}\beta & \beta & \beta & \beta \end{array}]$$, predicts a change in dimorphism of:

In this special case, where all genetic variances are equal, there is no direct response in dimorphism. The first term inside the parentheses, $$ r_{\delta w}$$, is the indirect response in asymmetry due to differences between male and female within-sex correlations; the second term, $$ r_{\delta b}$$, is due to asymmetry of the off-diagonal between-sex, between-trait correlations within the B matrix. Changing the direction of the concordant vector reverses the sign of the effects of the differences.

In general, traits will experience a mixture of antagonistic and concordant selection simultaneously. The overall rate of change in sexual dimorphism is a function of all the components of the G matrix. As the number of traits increases, the conditions under which the effects of concordant selection on dimorphism can exceed that of antagonistic selection of equal strength become more varied, as the asymmetries of the off-diagonal elements of $$ {\mathbf{G}}_M,\hspace{0.17em}{\mathbf{G}}_F$$, and B, as well as unequal genetic variances, can contribute to the indirect responses that affect dimorphism.

Gene-Wise Analysis

Of the 12,071 genes investigated, we detected significant genetic variation in 10,489 genes, similar to the results obtained by Ayroles et al. (2009). We do not consider the nonsignificant genes further. Dimorphism of expression of the ith gene is quantified as $$ D_i$$, the difference between means of log₂-transformed male and female gene expression. Of the significant genes, 1,447 were male biased (MB: $$ D_i>1$$, a ratio of male to female expression greater than 2), 2,073 were female biased (FB: $$ D_i<-1$$, a ratio less than 0.5), and the remaining 6,979 genes had relatively unbiased expression (UB: $$ {\left|D\right|}_i\le 1$$, within a factor of two between the sexes).

Which Aspects of G Correlate with Dimorphism of Transcription?

Under the familiar hypothesis that SAS drives the evolution of dimorphism, equations (3) and (6) predict that the absolute value of sexual dimorphism, $$ \left|D\right|$$, will be positively related to $$ \overline{g}$$ and negatively related to $$ r_{MF}$$, as long as the G matrix is consistent over a relevant evolutionary time scale. In addition, $$ {\overline{r}}_w \text{and} {\overline{r}}_b$$ have conflicting effects on the change in dimorphism. Their net effect can be positive or negative. If sexually concordant selection (SCS) influences dimorphism, equation (4) shows that differences in male and female genetic variances should be positively correlated with $$ \left|D\right|$$, whereas equation (7) predicts that the absolute value of the average differences of within- $$ (\left|{\overline{r}}_{\delta w}\right|)$$ and between-sex correlations $$ (\left|{\overline{r}}_{\delta b}\right|)$$ should also be positively correlated with $$ \left|D\right|$$. In addition, several authors (Mank et al. 2008; Assis et al. 2012; Griffin et al. 2013; Dean and Mank 2016) have shown that the mean level of gene expression, $$ \overline{E}$$, and the degree of tissue-specificity, τ (Yanai et al. 2005), are related to $$ \left|D\right|$$.

The Pearson correlation matrix of $$ \hspace{0.17em}{\text{log}}_{10}(\left|D\right|+0.01)$$ and these predictors, calculated from the sample covariance matrix of sex-specific line means for all 10,489 genes, G*, is shown in supplementary table S1, Supplementary Material online. Note that statistical tests of association are invalid for these data, as the underlying expression data are not independently estimated for each gene. These relationships are nonlinear, and in some cases nonmonotonic, as shown in figure 1 for the relationship between dimorphism and a subset of these predictor variables. The predicted relationships of $$ r_{MF} \text{and}\hspace{0.17em}\hspace{0.17em}{\text{log}}_{10}(\overline{g})$$ with dimorphism under SAS are evident. In addition, $$ {\overline{r}}_w$$ has a strong positive relationship with dimorphism, suggesting the effects of the within-sex between-trait correlations swamp those of the between-sex between-trait correlations captured by $$ {\overline{r}}_b$$ (see eq. 6). Both $$ \hspace{0.17em}{\text{log}}_{10}(\left|m-f\right|)$$ and $$ \left|{\overline{r}}_{\delta w}\right|$$ have clear positive relationships with dimorphism as predicted under SCS (eqs. 4 and 7).

Fig. 1.

Exploratory distribution, density, and smoothed trend plots for $$ \hspace{0.17em}{\text{log}}_{10}(|\Delta |+0.01)$$ and potentially predictive aspects of the G matrix. We chose the ggpairs function from R (R Core Team 2020) package GGally (Wickham 2016; Schloerke et al. 2020) to visualize the distribution and correlations among variables. These functions use the grammar of graphics philosophy (Wilkinson 2005) and emphasize visualization over quantification. Plots on the diagonal are density plots in which the horizontal axis is on the scale shown at the bottom of the plot, whereas the vertical axis is scaled to the maximum density. All other plots use the axes shown to the left as the vertical axis label, and below the figure as the horizontal axis label. Density estimates in the panels on the diagonal are kernel density estimates calculated in the geom_density function in ggplot2. Plots above the diagonal are functions ±95% CIs calculated in ggplot2, which calls a generalized additive model in the R package mgcv. Plots below the diagonal were visualized via 2D kernel estimates calculated by function geom_point_2d in the R package MASS (Venables and Ripley 2002). Color scale is based on quantiles of local density, and consequently, the scale varies from panel to panel. All calculations performed using function default parameters.

These raw relationships are confounded, as suggested by the complex relationships among predictor variables evident in supplementary table S1, Supplementary Material online, and figure 1. To help disentangle these factors, we used multiple regression of sexual dimorphism, measured as $$ \hspace{0.17em}{\text{log}}_{10}(\left|D\right|+0.01)$$, on the other variables, bootstrapped at the level of inbred lines. The results are shown in table 1. The model explains 40% of the variation in dimorphism. Both mean expression $$ \overline{E}$$ and tissue specificity τ have consistent positive effects on dimorphism, as shown by the positive signs of the bootstrap quantiles. $$ \hspace{0.17em}{\text{Log}}_{10}\left(\overline{g}\right), r_{MF} \text{and} {\overline{r}}_w$$ explain 2.1%, 15.5%, and 6.7% of the variance in $$ \hspace{0.17em}{\text{log}}_{10}(\left|D\right|+0.01)$$, respectively. In addition, there is also a consistent positive relationship of $$ \hspace{0.17em}{\text{log}}_{10}(\left|m-f\right|)$$ with $$ \hspace{0.17em}{\text{log}}_{10}(\left|D\right|+0.01)$$, as expected under SCS, explaining 1.1% of the variance. Expression properties $$ \overline{E}$$ and τ also have significant positive relationships with $$ \hspace{0.17em}{\text{log}}_{10}(\left|D\right|+0.01)$$, explaining 3.2% and 3.5% of the variation.

Table 1.

Estimates of Slopes from Multiple Regression of $$ \hspace{0.17em}{\text{log}}_{10}(\left|D\right|+0.01)$$ on Expression Characteristics, Bootstrapped over Genes.

		All Genes			Biasedd			Unbiasedd
Parametera	Pred.b	Median	Quantilesc 2.5%, 97.5%	R 2	Median	Quantilesc 2.5%, 97.5%	R 2	Median	Quantilesc 2.5%, 97.5%	R 2
$$ \overline{E}$$		0.04	0.03, 0.05	3.2	0.03	0.02, 0.03	6.2	0.03	0.03, 0.04	1.8
τ		0.78	0.49, 0.96	3.5	0.99	0.82, 1.12	32.6	−0.01	−0.10, 0.10	0.0
$$ \hspace{0.17em}{\text{log}}_{10}(\overline{g})$$	A+	0.20	0.02, 0.34	2.1	0.05	−0.04, 0.11	0.8	0.14	0.10, 0.19	1.2
$$ {\overline{r}}_{MF}$$	A−	−0.63	−0.76, −0.48	15.5	−0.10	−0.15, −0.04	4.5	−0.28	−0.38, −0.19	2.3
$$ \hspace{0.17em}{\text{log}}_{10}(\left|m-f\right|)$$	C+	0.11	0.07, 0.16	1.1	0.00	−0.02, 0.02	0.3	0.04	0.02, 0.07	0.4
$$ {\overline{r}}_w$$	A?	4.81	2.99, 7.66	6.7	0.41	0.19, 0.60	0.3	−0.95	−2.13, 0.50	4.2
$$ {\overline{r}}_b$$	A?	0.03	−1.52, 1.41	0.1	0.09	−0.08, 0.30	0.2	0.15	−0.49, 0.96	0.0
$$ \left|{\overline{r}}_{\delta w}\right|$$	C+	3.72	−6.10, 14.23	1.8	0.01	−0.44, 0.52	1.3	1.17	−1.58, 5.59	0.6
$$ \hspace{0.17em}{\text{log}}_{10}(\left|{\overline{r}}_{\delta b}\right|+0.01)$$	C+	−0.09	−0.23, 0.20	0.0	−0.00	−0.08, 0.13	0.4	−0.07	−0.14, 0.01	0.1

Note.—A, antagonistic prediction; C, concordant prediction;?, both $$ {\overline{r}}_w$$ and $$ {\overline{r}}_b$$ are predicted to affect dimorphism under SAS, but the sign of these effects depends on the details of selection.

a Parameter symbols explained in the text; $$ {\overline{r}}_w \text{and} {\overline{r}}_{b\hspace{0.17em}}$$ generalize the indirect selection parameters defined in equation (6) to include the mean of the average within-sex correlations of all traits with expression of the focal gene, whereas $$ \left|{\overline{r}}_{\delta w}\right| \text{and} \left|{\overline{r}}_{\delta b}\right|$$ generalize those in equation (7) to include the absolute values of the mean differences of all within-sex correlations with the focal gene, or the differences in the between-sex correlations involving the focal gene. $$ \left|{\overline{r}}_{\delta b}\right|$$ was log-transformed to minimize the influence of observations with exceptionally high values.

b Predicted sign of relationships with $$ \hspace{0.17em}{\text{log}}_{10}(\left|D\right|+0.01)$$ under concordant or antagonistic selection.

c Quantiles from 1,000 bootstrap resamples at the inbred line level. When the bootstrap 95% quantiles have consistent sign, we consider the effects to be statistically significant. Significant values are shown in bold face.

d Biased genes have $$ \left|D\right|>1$$; Unbiased genes $$ \left|D\right|\le 1$$.

Separate analyses of the biased $$ (\left|D\right|>1)$$ genes show that the G matrix properties are much less predictive of dimorphism in this subset. $$ \overline{E}$$ and τ remain significant, and the role of τ is much stronger in this class of genes, explaining 32.6% of the variation. Of the G matrix properties, $$ r_{MF}$$ and $$ {\overline{r}}_w$$ remain consistent predictors of dimorphism, although the variance explained is much reduced from that over the entire data set. Analysis of just the relatively unbiased $$ (\left|D\right|\le 1)$$genes shows that $$ \overline{E}$$ is again a significant predictor, but τ is not. G matrix properties expected to covary with $$ \hspace{0.17em}{\text{log}}_{10}(\left|D\right|+0.01)$$ under both SAS and SCS are significant in this subset; the relationships with $$ \hspace{0.17em}{\text{log}}_{10}\left(\overline{g}\right) \text{and} r_{MF}$$ suggest the action of SAS, whereas the positive relationship with $$ \hspace{0.17em}{\text{log}}_{10}(\left|m-f\right|)$$ is consistent with SCS. All of these variables explained a small proportion of the variance in dimorphism. On balance, these results indicate a strong signature of past SAS over the entire range of dimorphism values, and a more modest effect of SCS which is strongest for traits that have smaller amounts of dimorphism.

Innocenti and Morrow (2010) estimated selection on expression for a subset of the transcripts in our data set, but unfortunately did not report selection effect sizes. Transcripts with a significant sex by fitness interactions (their models used transcription as the dependent variable) experience some SAS, whereas those with a significant main effect of fitness experience some SCS. Note that selection may include both concordant and antagonistic components (Cheng and Houle 2020). When the identity of genes with a fitness main effect, or a sex-by-fitness interaction, were entered as predictors in the multiple regression model, these terms had small negative effects on $$ \left|D\right|$$ (median, 2.5 and 97.5 percentiles: main effect b= −0.034, −0.056, −0.011; interaction b= −0.057, −0.073, −0.034). Perhaps more importantly, the parameter estimates shown in table 1 were essentially unchanged by the fitness effect indicators. There is no clear predicted relationship between the selection measured in the lab by Innocenti and Morrow (2010) and dimorphism. One potential explanation for the negative relationship between selection and dimorphism is that Innocenti and Morrow’s study had greater power to detect selection in less dimorphic transcripts.

Quantitative Genetic Analysis of Gene Expression

Ayroles et al. (2009) measured gene expression in just 40 inbred lines, which precludes estimating a G matrix for more than a small number of traits. To choose these traits, we performed principal component analyses within the male-, female-, and relatively-unbiased gene classes, as described in the Methods. Through exploratory analyses described in supplementary results, Supplementary Material online, we chose the first two PCs of male-biased genes (symbolized MB), and of female-biased genes (FB), and the first four PCs for relatively unbiased genes (UB). These analyses revealed no evidence for genetic variance in female expression of MB genes, or for male expression of FB genes, suggesting that direct selection to alter expression in the low-expressing sex will be relatively ineffective. Accordingly, we dropped these traits from G. The transcriptional modules inferred by Ayroles et al. (2009) are not generally associated with the sex-bias class PCs used in the quantitative genetic analysis (see supplementary methods, results, and table S4, Supplementary Material online).

The best-fitting model suggested that genetic variation in the 12 traits can be explained well by a 9 dimensional model, resulting in the G matrix shown in table 2. The resulting data set consisted of 12 traits, although only the four UB (UB1–UB4) traits were estimated in both sexes. Thus, only the portions of table 2 involving these traits (outlined in boxes) correspond to Lande’s (1980)B matrix. For these UB traits, $$ r_{MF}$$ averages 0.82. This is higher than the average value of 0.75 found in a review of previous studies (Poissant et al. 2010), and much higher than the average r_MF value of estimated in Griffin et al. (2013). Three of these correlations are less than 3 SE from a perfect correlation of 1, suggesting that there may be very little variation to generate opposing responses in the sexes for these traits (Sztepanacz and Blows 2017).

Table 2.

Genetic Correlation and Covariance Matrices from the 12 Trait Analyses.

Table 2. Genetic Correlation and Covariance Matrices from the 12 Trait Analyses.

Note.—Genetic variances are shown in bold face on the main diagonal. Genetic correlations are above the main diagonal, and genetic covariances are below. Sampling standard errors are shown in italic font on the line below the estimates. The boxes outline the G_M, G_F, and B submatrices for the four relatively unbiased traits where we estimated covariances between traits in both sexes. The diagonal elements of B are shown in bold italics.

Examination of table 2 also reveals several other correlations indicative of constraints on the evolution of biased gene expression. The most striking of these is female expression of FB1 and of UB1 where $$ r=-0.96$$. In addition, many of the correlations between MB traits and male expression of the relatively unbiased traits have absolute values in the neighborhood of 0.4–0.5, as does the correlation of FB2 with female expression of UB3. These elements show the potential for complex constraints between unbiased and biased genes.

Predicted Responses to Selection

We used G to predict the effects of current selection on changes in dimorphism from its current level. We predicted the responses to symmetrical antagonistic and concordant versions of five different selection gradients: directional selection on each of the four UB traits individually, and directional selection on all four simultaneously with the results shown in table 3. We show two different estimates of the amount of evolution under each selective regime. Evolvability, e, is the response in the direction of the selection gradient, whereas respondability, R, is the length of the total response to selection (Hansen and Houle 2008). As expected from the positive r_MF, the responses to SCS selection are overall larger than the responses to SAS.

Table 3.

Responses of Relatively Unbiased Expression Traits to Antagonistic (A) and Concordant (C) Selection of Equal Strength $$ (\Vert \mathbf{\beta }\Vert =1)$$.

	Antagonisticb		Concordantc		UB $$ \Vert \Delta \Vert $$d			Vector correlation. $$ \Delta {\overline{z}}_{\mathrm{M}}\hspace{0.17em}\text{vs}.\hspace{0.17em}\Delta {\overline{z}}_{\mathrm{F}}$$
Sel.a	E	R	e	R	A	C	ratio	A	C
UB1	35.6 (21.0–52.9)	60.6 (35.6–89.2)	90.8 (53.1–133.7)	104.7 (62.9–155.0)	27.4 (16.3–40.5)	31.2 (16.0–49.4)	0.87 (0.66–1.26)	0.27 (−0.22 to 0.60)	0.83 (0.65–0.97)
UB2	7.9 (4.1–12.7)	24.9 (14.0–38.2)	87.7 (55.1–132.0)	94.1 (60.0–140.6)	10.4 (5.7–16.1)	16.1 (6.6–29.7)	0.63 (0.34–1.46)	0.34 (−0.09 to 0.60)	0.89 (0.73–0.98)
UB3	4.2 (1.9–7.5)	15.2 (6.8–24.9)	84.4 (50.6–125.7)	88.6 (54.0–130.1	6.3 (2.8–10.6)	7.5 (3.0–16.5)	0.81 (0.29–2.24)	0.33 (−0.30 to 0.73)	0.97 (0.88–1.00)
UB4	7.0 (3.7–11.5)	16.2 (8.9–25.8)	73.0 (45.2–108.8)	77.4 (48.5–114.3)	6.2 (3.2–10.7)	7.2 (2.3–18.7)	0.86 (0.28–2.85)	0.39 (−0.52 to 0.80)	0.97 (0.82–1.00)
UB	27.6 (17.3–40.4)	48.9 (28.4–73.3)	70.4 (44.7–103.5)	75.9 (48.7–110.6)	23.8 (14.7–34.7)	13.2 (3.8–26.2)	1.80 (0.93–5.42)	0.10 (−0.46 to 0.46)	0.92 (0.80–0.99)

Note.—Values are medians (2.5–97.5% quantiles). UB, all four UB traits are simultaneously selected; e, evolvability, the response in the direction of selection; R, respondability, the total response to selection; A, C, total change in dimorphism under A or C selection; ratio, $$ \Vert {\Delta }_A\Vert /\Vert {\Delta }_C\Vert $$.

a Selection regime: symbols indicate trait subject to directional selection.

b Selected male traits have positive gradients, whereas female traits negative selection gradients.

c All selected traits have positive gradients in both sexes.

d Predicted change in length of dimorphism vector.

The key prediction concerns the total change in sexual dimorphism, the change in the length of the multivariate vector of differences between the sexes, $$ \Vert \Delta \Vert $$. Surprisingly, $$ \Vert \Delta \Vert $$ is often larger under SCS than under SAS. Although the confidence limits of the ratios of $$ \Vert \Delta \Vert $$ due to SAS and SCS are not significantly different from one for any of the selection gradients investigated, these results contrast with the generally assumed scenario that sexual dimorphism is the result of direct selection for dimorphism. As shown in the last two columns of table 3, responses of each sex to SAS are positively correlated, rather than negatively correlated, as selection favors.

To get a sense for how much dimorphism can be created by concordant selection, we can compare the direct response to concordant selection with the indirect response in dimorphism in table 3. These ratios range from 0.34 to 0.09 for the four UB traits. The total change in dimorphism when trait UB1 is concordantly selected, for example, is more than 1/3 as much ($$ \Vert \Delta \Vert =31.2$$) as the total concordant change (e = 90.8) in trait UB1. Thus, large changes in average gene expression can create smaller but still substantial total changes in dimorphism.

We also predicted the response of dimorphism of the relatively unbiased traits to selection on the sex-biased MB and FB traits with results shown in supplementary table S5, Supplementary Material online. The total changes in dimorphism under these scenarios are of comparable magnitudes to the changes arising from selection on the UB traits. Indirect responses of dimorphism to selection on sex-biased traits can have important effects on other traits.

The scenarios in table 3 and supplementary table S5, Supplementary Material online, assume that there is no selection on the traits not under directional selection. Alternatively, we can assume that traits not under directional selection are subject to such strong stabilizing selection that a response is only possible in the direction of the selection gradient (Hansen et al. 2003; Hansen and Houle 2008). Under this sort of selection, SAS and SCS again have similar effects on overall dimorphism as shown in supplementary table S6, Supplementary Material online, although the responses are much less than for the corresponding unconditional scenarios in table 3.

The results in table 4 show how these predictions are affected by the symmetry of G. To make these comparisons, we substituted the symmetrical version of $$ {\mathbf{G}}_{\mathrm{M}} \text{and} {\mathbf{G}}_{\mathrm{F}}$$ (their average, $$ \overline{\mathbf{G}}$$) and B (B_S) into the G matrix and compared the predicted response to those from the unaltered matrix. Substituting $$ \overline{\mathbf{G}}$$ reduces the change in dimorphism in response to SCS by an average of 30.4% over our five selective scenarios. Symmetrizing B can either increase or decrease the dimorphic response to SCS, but the average changes in dimorphic response are reduced by 28.9%, nearly as large as the effect of symmetrizing G_M and G_F. Interestingly, the asymmetries due to B, G_M, and G_F can have opposite effects. For example, simultaneous SCS on all UB traits leads to smaller changes in dimorphism relative to that under SAS than any other selection regime. The results in table 4 show that symmetrizing B increases dimorphism under SCS by 41%, while symmetrizing G_M and G_F decreases dimorphism by 41%. The change in dimorphism is small because the asymmetries work against each other. The predicted response in dimorphism in UB2, however, is large because both types of asymmetry affect dimorphism in the same direction. In addition, the effects of these asymmetries are highly nonlinear in combination, as eliminating both kinds of asymmetry eliminates the evolution of dimorphism under concordant selection, regardless of how each asymmetry separately affects dimorphism. As predicted, symmetrizing the matrices has no effect on the change in dimorphism under SAS. The B matrix as a whole is an important constraint to the evolution of dimorphism, as eliminating it always increases the change in dimorphism under SAS.

Table 4.

Ratio of Changes in Relatively Unbiased Transcript Dimorphism ($$ \Vert \Delta \Vert $$) Predicted from a Modified G Matrix Relative to Predictions from the Unmodified G Matrix.

	(UB $$ \Vert \Delta \Vert $$ Modified G)/(UB $$ \Vert \Delta \Vert $$ Unmodified G)b
	$$ \left[\begin{array}{cc}\overline{\mathbf{G}}& \mathbf{B}\\ {\mathbf{B}}^{\mathbf{T}}& \overline{\mathbf{G}}\end{array}\right]$$		$$ \left[\begin{array}{cc}{\mathbf{G}}_{\mathrm{M}}& {\mathbf{B}}_S\\ {\mathbf{B}}_S& {\mathbf{G}}_{\mathrm{F}}\end{array}\right]$$		$$ \left[\begin{array}{cc}{\mathbf{G}}_{\mathrm{M}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{G}}_{\mathrm{F}}\end{array}\right]$$		$$ \left[\begin{array}{cc}\overline{\mathbf{G}}& {\mathbf{B}}_S\\ {\mathbf{B}}_S& \overline{\mathbf{G}}\end{array}\right]$$		$$ \left[\begin{array}{cc}\overline{\mathbf{G}}& \mathbf{0}\\ \mathbf{0}& \overline{\mathbf{G}}\end{array}\right]$$
Sel.a	A	C	A	C	A	C	A	C	A	C
UB1	1.00	0.37	1.00	0.93	1.81	0.91	1.00	0.00	1.68	0.00
UB2	1.00	0.77	1.00	0.53	2.53	0.27	1.00	0.00	3.13	0.00
UB3	1.00	0.94	1.00	1.32	4.67	1.04	1.00	0.00	4.94	0.00
UB4	1.00	0.81	1.00	0.84	5.21	0.91	1.00	0.00	5.36	0.00
UB	1.00	0.59	1.00	1.41	1.33	1.38	1.00	0.00	1.48	0.00

a Selection regime (as in table 3).

b See text for explanation of modified G matrices.

Gene Expression Data

We reanalyzed the adult gene expression (Ayroles et al. 2009) in inbred lines of the Drosophila Genome Reference Project (Mackay et al. 2012). Each DGRP line was independently derived from a different inseminated female D. melanogaster sampled from a single outbred population, followed by 20 generations of brother–sister mating.

Ayroles et al. (2009) assayed whole-body gene expression using Affymetrix Drosophila Genome 2.0 microarrays. They assayed expression twice in each sex in each of 40 DGRP lines, for a total of 160 chips. The data were normalized and processed using the R Bioconductor package oligo (Carvalho and Irizarry 2010; Huber et al. 2015). In some cases, a single probe assayed expression of transcripts at more than one gene. We assigned results from such probes to just one of these genes, chosen arbitrarily. This left expression data on 12,701 genes. Gene expression was in log₂ units, so differences in expression are equivalent to log₂(ratios).

Gene-by-Gene Analyses

We tested for the presence of significant genetic variation at each gene using a mixed model analysis with sex, probe, and probe-by-sex effects fixed, and line, line-by-sex and line-by-probe effects treated as random. For 11,039 genes, a single probe was assayed, and the probe effects were omitted from the model. We compared the likelihood of the data for models that included or omitted all the random effect terms using a likelihood ratio test with 2 df for the genes with only one probe, and three df for genes with more than one probe. Genes were retained for further analysis if P < 0.01 for the total line effects. Mixed model analyses were fit using Proc Mixed in SAS/STAT software (SAS Institute Inc 2016).

For genes with significant line effects, we retained the least squares means from this model for each sex and line, and then calculated the covariance matrix of these means to form G*. Sexual dimorphism in expression of the ith gene, $$ D_i={\overline{z}}_{M·i}-{\overline{z}}_{F·i}$$, was measured as the difference between the least squares means of male and female log₂-transformed gene expression. We defined three categories of sex bias: male-biased (MB) genes with male expression more than twice as large as female expression, $$ D_i>1$$, female-biased (FB) genes with female expression more than twice as large as male expression, $$ D_i<-1$$, and relatively unbiased (UB) genes with $$ \left|D_i\right|\le 1$$.

From equation (3), we predict that if antagonistic selection plays a major role, $$ \left|D_i\right|$$ will be positively related to $$ {\overline{g}}_i$$, and negatively related to $$ r_{\text{MF.}i}$$. From equation (4), we predict a positive relationship between $$ \left|m_i-f_i\right|$$ and $$ \left|D_i\right|$$ under concordant selection. From equation (6), we can see that the averages of the corresponding off-diagonal elements of $$ {\mathbf{G}}_M^{*} \text{and} {\mathbf{G}}_F^{*}$$, and of off-diagonal elements above and below the diagonal of B* affect the evolution of dimorphism under antagonistic selection. To capture the effect of within-sex correlations for the ith gene, we calculated:

From equation (7), we can see that differences of the corresponding off-diagonal elements of $$ {\mathbf{G}}_{\text{M}}^{*} \text{and} {\mathbf{G}}_{\text{F}}^{*}$$, and of off-diagonal elements above and below the diagonal of B* promote the evolution of dimorphism under concordant selection. We quantified the effects of the differences in within-sex correlations as:

In addition, we also investigated the effects of the mean level of expression, $$ {\overline{E}}_i$$, and an index of the tissue specificity of gene expression, τ_i (Yanai et al. 2005), on the level of sexual dimorphism. We downloaded tissue-specific expression data from http://flyatlas.org/ Geo accession GSE7763 (last accessed January 5, 2021), and calculated τ for each transcript, then averaged them to obtain τ_i for gene i.

Values of $$ \left|D_i\right|$$ must be positive and are strongly right-skewed with a small minority of very large values. To reduce the influence of this large tail, while preventing values very near 0 from causing a left-skew to the transformed data, we transformed dimorphism as $$ \hspace{0.17em}{\text{log}}_{10}(\left|D_i\right|+0.01)$$. We log₁₀-transformed $$ {\overline{g}}_i$$, $$ \left|m_i-f_i\right|$$ and $$ \left|{\overline{r}}_{\delta b.i}\right|$$ before analysis for similar reasons.

Selection on Gene Expression in Drosophila melanogaster

The nature of selection on transcripts was derived from table S1 in Innocenti and Morrow (2010). For 15 genotypes with high and low male and female fitnesses, they tested whether fitness predicted the expression of each transcript, and for a sex by fitness interaction. They only report P values for tests judged to be significant, so it is unclear how many total transcripts were tested. When paired with the Ayroles et al. (2009) data set, 516 genes had least one transcript with a significant main effect, 1,300 had a significant interaction, and 348 had both effects significant.

Quantitative Genetic Analysis of Bias Classes

We performed separate principal component analyses (PCA) on the covariances of least-squares means of expression in each sex for the male-, female, and relatively-biased genes. For the unbiased genes, we performed a PCA on the covariance matrix of line-sex averages. For the two biased sets, we performed a PCA of expression in the dominant sex. We partitioned the variance in the principal component scores into genetic and nongenetic sources using restricted maximum-likelihood implemented in the program Wombat (Meyer 2006–2019). We assumed that all the inbred lines were unrelated. Whole-genome sequencing of these lines suggests that this is a good, albeit imperfect, approximation of the relationship among them (Huang et al. 2014). Estimation of the genetic variance–covariance matrix, G, was carried out for both full- and reduced-rank models (Kirkpatrick and Meyer 2004; Meyer and Kirkpatrick 2005, 2008), and we selected the best-fitting model on the basis of Akaike’s information criterion corrected for small sample size (AICc).

We assessed the fit of models with different numbers of traits drawn from the three expression bias classes, eventually settling on a 12-trait data set as described in the Results section. Sampling variances of matrix elements remained stable under models fit to even fewer traits. Once we obtained well-estimated matrices, we back-transformed estimates to the original scores and used the REML-MVN approach (Meyer and Houle 2013; Houle and Meyer 2015) to generate 1,000 replicate matrices drawn from the sampling distribution of the matrix G. These replicate estimates of G were used to generate the sampling distribution of the predicted responses to selection. In most cases, the distributions of G-derived statistics were asymmetrical, so we report medians and 2.5% and 97.5% quantiles.

Analyses of Modified G Matrices

To investigate which aspects of G matrix structure have effects on the evolution of dimorphism, we formed five modified matrices that highlight those aspects of G that affect the evolution of dimorphism. The k = 2 theory developed above, as well as the more general analyses of Cheng and Houle (2020) suggests that the evolvability of dimorphism to antagonistic selection is promoted by the average genetic variance in males and females, $$ \overline{\mathbf{G}}=({\mathbf{G}}_M+{\mathbf{G}}_F)/2$$, and restricted by the symmetrical component of $$ \mathbf{B},\hspace{1em}{\mathbf{B}}_S=(\mathbf{\text{B}}+{\mathbf{B}}^T)/2$$. When these matrices are substituted for $$ {\mathbf{G}}_M, {\mathbf{G}}_F,\hspace{0.17em}\mathbf{B} \text{and} {\mathbf{B}}^T$$ this eliminates the components of G that cause evolvability of dimorphism under concordant selection, namely the asymmetry between $$ {\mathbf{G}}_{\mathrm{M}} \text{and} {\mathbf{G}}_{\mathrm{F}},\hspace{0.17em}({\mathbf{G}}_{\mathrm{M}}-{\mathbf{G}}_{\mathrm{F}})/2$$ and the asymmetrical component of $$ \mathbf{B},\hspace{0.17em}{\mathbf{B}}_A=(\mathbf{B}-{\mathbf{B}}^T)/2$$. In addition, we explored the effect of eliminating B by replacing it with the zero matrix 0.