Experimental evolution of a pheromone signal

Abstract Sexual signals are important in speciation, but understanding their evolution is complex as these signals are often composed of multiple, genetically interdependent components. To understand how signals evolve, we thus need to consider selection responses in multiple components and account for the genetic correlations among components. One intriguing possibility is that selection changes the genetic covariance structure of a multicomponent signal in a way that facilitates a response to selection. However, this hypothesis remains largely untested empirically. In this study, we investigate the evolutionary response of the multicomponent female sex pheromone blend of the moth Heliothis subflexa to 10 generations of artificial selection. We observed a selection response of about three‐quarters of a phenotypic standard deviation in the components under selection. Interestingly, other pheromone components that are biochemically and genetically linked to the components under selection did not change. We also found that after the onset of selection, the genetic covariance structure diverged, resulting in the disassociation of components under selection and components not under selection across the first two genetic principle components. Our findings provide rare empirical support for an intriguing mechanism by which a sexual signal can respond to selection without possible constraints from indirect selection responses.

constraints to their selection response and identify the mechanisms that can mitigate those constraints.
Understanding selection responses in sexual signals is challenging because signals are often composed of multiple components (Candolin, 2003;Higham & Hebets, 2013;Rowe, 1999). For example, mating songs can vary both in pitch and in rhythm (Blankers, Hennig & Gray, 2015;Tanner, Ward, Shaw & Bee, 2017;Wilkins, Shizuka, Joseph, Hubbard & Safran, 2015), color signals can be composed of multiple, functionally distinct patches (Cole & Endler, 2015;Grether, Kolluru & Nersissian, 2004), and sex pheromones are often blends of multiple chemical compounds (Ferveur, 2005;Linn, Campbell & Roelofs, 1987). To understand how multicomponent signals evolve, we thus need to consider the selection response in multiple dimensions simultaneously. Moreover, signal components can have a shared genetic or developmental basis, or can be subject to correlated selection pressures (Armbruster, Pélabon, Bolstad & Hansen, 2014;Cheverud, 1996). The resulting genetic correlations among signal components can influence how selection on the phenotype translates to changes in the underlying genotypes (Chenoweth & Blows, 2006). To understand how the genotype-phenotype map of sexual signals influences the selection response, we thus need to determine the genetic correlations between the different signal components.
Statistical frameworks in quantitative genetics, in particular, the (multivariate) breeder's equation, allow us to predict and quantitatively understand selection responses in correlated traits. In this framework, the response to selection is a function of the genetic (G) and phenotypic (P) variance-covariance matrix and the selection gradient: selection acts on the P matrix and the resulting response is constrained by the G matrix (Lande, 1979;Lande & Arnold, 1983;Lynch & Walsh, 1998). The difficulty in predicting selection responses of multivariate traits is that selection acting on multiple components may be counterbalancing, e.g., directional selection on one trait, but stabilizing selection on correlated traits, resulting in evolutionary constraints (Barton & Turelli, 1989). Counterbalancing selection is likely prevalent in the evolution of sexual signals, as choosing individuals may favor higher or lower values of some component of the signal, while changes in correlated components may result in reduced mate attraction.
In this study, we explored the response in phenotypic means and genetic (co)variances to artificial selection on the female sex pheromone of the moth Heliothis subflexa (Lepidoptera, Noctuidae). Like many other moths, H. subflexa females secrete a sex pheromone blend to which conspecific males are attracted. These sex pheromone blends are species specific and vary among species in both the presence/absence of components as well as in relative amounts (or ratios) of the components (Cardé & Haynes, 2004;Schneider, 1992).
To explore the selection response of the acetates, we selected for higher and lower amounts of acetates during 10 generations of truncation selection. Since the acetates vary geographically (Groot, Inglis, et al., 2009) and have a genetic basis that is partially independent of other components (Groot, Estock, et al., 2009), we hypothesized that the relative amount of acetates can evolve in response to univariate selection for higher/lower acetates but that genetic variance will be reduced after selection. However, since acetates also partially share their genetic basis with other components and since all sex pheromone compounds are produced through the same biosynthetic pathway, we also hypothesized that there will be indirect selection responses in other pheromone components, specifically in the unsaturated aldehydes, Z9-16:Ald and Z11-16:Ald, and alcohol, Z11-16:OH ( Figure 1). Since some of the correlations among components may reduce the selection response, we also expected genetic covariances to change during selection, specifically in a way that facilitates the phenotypic selection response.

| Insect rearing
The laboratory population of H. subflexa was found from animals col- USA) at room temperature for approximately 10 days, after which larvae were reared in separate individual 37-ml cups filled with the same artificial diet and kept at 25°C and 60% relative humidity with 14 h:10 h light-dark cycle. Upon emergence, adults were provided sugar water. Pairs of males and females were housed in 375-ml paper cups covered with gauze and kept under the same conditions as the late instar larvae and virgin adults. The mating pairs were provided with sugar water. To stimulate oviposition, a freshly cut gooseberry fruit was placed on top of the gauze. Once the eggs began to hatch, the gauze and the eggs and larvae on it were transferred to Petri dishes, which were then placed at room temperature until larvae were transferred to individual cups. The females that produced fertile offspring were collected and phenotyped, if still alive. All matings were assigned unique numbers. This way, we obtained a full pedigree of all individuals in the selection and control lines.

| Phenotyping
Female sex pheromone can be extracted from mated and old females by injecting them with pheromone biosynthesis activating neuropeptide, that is, PBAN (Groot et al., 2005). Females were injected with a 7.5 pmol (2 µl of a 0.0146 µg/µl) PBAN solution to activate pheromone production post-mating. After a 90-min incubation time, female glands were extruded by squeezing the abdomen then fixed by firmly holding the abdomen with forceps just anterior of the gland. The gland was excised with microdissection scissors, and the moths were euthanized.
F I G U R E 1 Simplified biosynthetic pathway of the Heliothis subflexa sex pheromone, redrawn from Groot, Estock, et al. (2009) and based on Jurenka (2003). Desaturation and β-oxidation produce mono-unsaturated acyl-CoA precursors from 18 or 16 carbon acyl CoA derivatives, which are then modified to form acetates, alcohols, and aldehydes through specific enzymatic conversions. The ∆9 and ∆11 desaturases create a double bond between the 9th and 10th or 11th and 12th carbon of the 18-or 16-carbon acyl-CoA derivatives, respectively. β-oxidation shortens the chain length from 18 to 16 carbons. Note that (Z)9-16 acyl-CoA can be formed through two alternative routes. The compounds present in the pheromone glands of H. subflexa are in boxes that are color coded depending on their function in male response behavior Excess abdominal tissue and eggs that remained in the ovipositor were removed, after which the glands were submerged in 50 μl hexane containing 200 ng pentadecane as an internal standard. After 30 min, the glands were removed and the extracts stored at −20°C until analysis.
Relative amounts were calculated by dividing the absolute amounts by the total amount across all 11 components.

| Data transformation
Selection was performed based on relative amounts of acetates.
However, describing relationships among relative amounts is problematic because they sum to 100%, thereby mathematically constraining the (co)variation in the pheromone and biasing the analysis. We therefore transformed pheromone measurements to log-contrasts for all down-stream analyses in this study. This approach breaks the interdependency and normalizes the data. Since the divisor used in the contrast of the variable can no longer be part of any downstream analyses, we chose 14:Ald as the divisor because this component has a small but clearly detectable peak in the chromatogram, while it is irrelevant for male response behavior (Heath, Mitchell & Cibrian Tovar, 1990). Prior to downstream analyses, samples with a χ 2 -distributed Mahalanobis distance score (calculated using the "mahalanobis" function in the "stats" package) that exceeded a threshold value corresponding to a Bonferroni-corrected p < .05 were removed, which resulted in removing 64 of a total of 2861 samples. These samples showed abnormal pheromone ratios and were present across all selection and control lines and generations. We are therefore confident that these samples represent outliers, e.g., due to extraction or measurement errors, and are not representative of any relevant biological phenomena.

| Selection
Since phenotyping consisted of extracting the sex pheromone gland invasively, selection was performed post-mating. Each generation, we continued with those families that had a maternal phenotype satisfying an increasingly stringent threshold. Initially, a so-called "high" and "low" line were started with offspring from females with a relative amount of acetates above 22% or below 16%, respectively.
These values were based on the distribution of the relative amounts of acetates in the starting population, representing the first and third quantiles. These thresholds were kept for the first three generations.
In subsequent generations, we increased these thresholds to >24% or <14% (generations 4-5), and >26% or <12% (generations 6-9). The high line consisted of a total of 2236 breeding females across nine generations, ranging from 189 to 295 matings per generation. After outlier removal, a total of 1234 high-line females were phenotyped, or between 89 and 171 per generation during the selection phase.
The low line consisted of 2250 breeding females (189-296 per generation) and, after outlier removal, phenotypes were measured for a total of 1180 females (57-169 per generation). Because maintenance of the lines required individual rearing (due to larval cannibalism) and phenotyping of so many individuals, we were unable to maintain replicate lines of the high and low line. Parallel results between the lines in the selection response thus lend confidence to the patterns observed. However, we avoid drawing conclusions based on differences between the high and low lines because these can both reflect differential selection effects or sampling variance. We stopped selecting in generation 10, but continued to phenotype ~25 females per line and per generation until generation 13. Throughout the selection experiment, we also maintained a control line from which we phenotyped between 6 and 111 females (median 20 females) per generation. For detailed sample sizes, see Table A1 in the appendix.

| Phenotypic selection response
To test the hypothesis that acetates can respond to selection and that other pheromone components change due to indirect selection, the selection responses as measured by the log-contrasts were visualized using the R-package "ggplot." In addition, we also calculated the average difference between selected and control lines in units of standard deviation of the starting population. To test whether log-contrast ratios for pheromone components had significantly increased or decreased during selection, we compared the starting generation with the final generation for the control, high, and low line using a Student's t-test.

| Genetic (co)variance selection response
To test whether genetic variances decreased with selection and whether the structure of the genetic variance-covariance matrix was affected by selection, we first ran multi-response animal models implemented in the R-package "MCMCglmm." We limited our analyses to the four components that have been shown to affect male response: Z9-16:Ald, Z11-16:Ald, Z11-16:OAc, and Z11-16:OH (in all cases, the log-contrast to 14:Ald was used). Model formulation roughly followed the MCMCglmm manual (Hadfield, 2010) and Jarrod Hadfield's course notes (Hadfield, 2012), as well as the animal model tutorial by Pierre de Villemereuil (Villemereuil, 2012). Highand low-line individuals had separate pedigrees, but shared ancestors for which we had pedigree data up to three generations prior to the onset of the experiment (in total 1065 breeding females). To minimize the influence from priors on the covariances, for which we had no prior expectations with high degrees of belief, we formulated flat, uninformative priors for the variance-covariance structure of the random effect ("animal"). To evaluate the models, we checked chain convergence and assessed effective sample sizes and levels of autocorrelation. To check whether the prior did not contribute unduly to posterior estimates, we also compared genetic (co)variance estimates among univariate, bivariate, and full (tetravariate) models and using different, more informative priors. We found models with flat priors to perform best.
We obtained estimates of genetic (co)variances for the starting populations (generations 0 and 1 combined), for generations 2 and 3 combined (early response), and for generation 9 (final response). ), which provides a standardized measure of the evolvability of a trait that is independent of other variance components. We also examined changes in the G matrix by inspecting the trait loadings on genetic principle components. The four genetic principle components were obtained for every posterior sample using the "eigen" function. We focused on the first two axes, gPC1 (also known as g max , the direction in phenotypic space containing the largest fraction of the genetic variance) and gPC2, because these axes jointly described >90% of the genetic variance (see Results). Comparisons of posterior distributions were done using the honest posterior density. Posteriors were considered statistically significantly different if 90% HPD intervals did not overlap.

| Selection response in sex pheromone
All three acetates responded to selection for higher/lower relative total amounts of acetates in the high and low line, respectively ( Figure 2). The line-specific means significantly increased in the high line for all three acetates between generations 0 and 9, and significantly decreased in the low line, while no significant differences were observed for the back-up line ( Table 1). The response was more pronounced in the more abundant biochemically related Z9 and Z11 isomers compared to the less abundant and biochemically separated Z7-16:OAc (Figure 2). The selection response was characterized by an immediate response in the direction of selection (from generation 0 to generation 1), followed by a gradual progression toward increasing/ decreasing relative amounts in the respective selection lines. In the low line, the selection response flattened after cessation of selection, while in the high line, the average log-ratio of Z11:16:OAc to 14:Ald was as high or higher in generations >10 compared to <10 (Figure 2d). The other pheromone components showed a significant decrease over time in both selected and control lines, but no differentiation among selection lines or between selected and control lines (Table 1; Figure 3a-c). The total amount of pheromone measured across the 11 biologically active components remained constant through time in the selected lines, indicating that changing ratios were independent of the pheromone titer (Figure 3d).

| Selection response in genetic (co)variances
In testing whether the selection response was associated with (i) a reduction in genetic variance in Z11-16:OAc (a change in the diagonal elements of the G matrix) and (ii) reorientation of genetic covariances among pheromone components (changes in the off-diagonal elements of the G matrix), we found no evidence for decreasing genetic variance in Z11-16:OAc or in the other components ( Figure   A1 in the appendix).
In contrast to genetic variances, we did find changes in the genetic covariance structure across generations. The first two genetic principle components showed a non-significant trend (differences in the mode, but overlapping 90% HPD intervals) toward change in magnitude of their eigenvalue, i.e., the amount of genetic variance they describe ( Figure A2). More of the total genetic variance across the four components tended to be in the direction of gPC1 in the high and low lines (posterior mode ranging between 80 and 90%) compared to the starting generation (70%). This means that more genetic variance of the sex pheromone was oriented toward a single dimension in phenotypic space.

| DISCUSS ION
In this study, we investigated how the selection response depends on and shapes the genetic architecture of a multicomponent sex pheromone signal. Through truncation selection, we gradually increased  blend of H. subflexa across 10 generations. As we hypothesized, we found that the acetates responded readily to selection, diverging three-quarters of a standard deviation away from the control line in 10 generations in both the high and low lines. However, in contrast to our hypothesis that selection would also result in indirect responses in the other traits, we found that the response was limited to the acetates only. In addition, we found that levels of genetic variance did not decrease, opposite to the expected erosion of genetic variance in response to selection. Lastly, in line with our expectations, the genetic covariance structure changed in a way that likely facilitated a selection response in the pheromone. We discuss each of these results further below.

| Univariate responses in multicomponent signals
In nature, sexual signals are often found to be subject to multivariate stabilizing and directional selection (Bentsen, Hunt, Jennions & Brooks, 2006;Blankers et al., 2015;Brooks et al., 2005 responses. In contrast, we thus observed independent evolution of coupled traits, which is not an uncommon result of artificial selection experiments (Hill & Caballero, 1992;Saltz, Hessel & Kelly, 2017), including for sexual signals (Ritchie & Kyriacou, 1996).

| Evolution of the genetic covariance structure
Changes in genetic covariances are rarely documented in natural or laboratory selected populations. Here, we found divergence (differences between before and after selection) in the loadings of the pheromone components on the genetic principle components, in particular in the magnitude of the correlation between the component under selection and gPC1 and gPC2. The observed divergence in the genetic covariance structure of the female sex pheromone may result from various mechanisms.
First, the covariance structure may change due to changes in genetic variances because stabilizing or directional selection is expected to erode genetic variance (Barton & Turelli, 1989;Estes & Arnold, 2007). Although our truncation selection increased and decreased the phenotypic mean by three-quarters of a standard deviation in the high and low line, respectively, we observed no significant decrease in the coefficient of additive genetic variance. One explanation is that several genetic factors may underlie the trait under selection and that selection therefore leaves no detectable signature on the genetic variance associated with that trait (Bulmer, 1971;Johnson & Barton, 2005). Genetic mapping studies have revealed several QTL underlying difference in sex pheromone composition in H. subflexa and even a single conversion step in the biochemical pathway, such as from Z11-16:OH to Z11-16:OAc, could be catalyzed by multiple enzymes (Groot et al., 2013). In addition, variation in pheromone composition in H. subflexa is explained by fitness variation (Blankers et al., 2021), indicating that many different genetic factors likely contribute small additional fractions to the total variance. Another explanation comes from population genetics: reorientation of genetic covariances in response to bottlenecks may free up additive genetic variance, thereby paradoxically increasing levels of genetic variance (Carson, 1990;Templeton, 2008). Since truncation selection is effectively a non-random bottleneck of the population, this may also explain the maintenance of genetic variance that we observed. Lastly, to avoid inbreeding depression, we maintained large populations and avoided mating first-and second-degree relatives. Such a mating scheme likely counteracts some loss in genetic variability due to selection (Du, Bernstein, Hoppe & Bienefeld, 2021).
Second, the genetic covariance structure can evolve due to changes to interactions among unlinked genetic loci affecting the coexpression of multiple traits (Wolf et al., 2005). Since there is a variety of enzymes involved in the conversion between pheromone components (Groot, Dekker & Heckel, 2016;Jurenka, 2003;Tillman, Seybold, Jurenka & Blomquist, 1999) and variation in different pheromone components is due to effects at different chromosomes (Groot, Estock, et al., 2009;Groot et al., 2014), it is likely that a significant portion of genetic covariance among sex pheromone components results from interactions among unlinked loci. It may therefore be unsurprising that we see an immediate change in the covariance structure following the onset of selection. This finding is of broad significance to predicting selection response from quantitative genetics because these predictions often depend on assumptions of a constant matrix over multiple generations. Our findings lend support to both theoretical (Arnold et al., 2008;Eroukhmanoff, 2009;Roff, 2000) and empirical (Björklund et al., 2013;Hine, Chenoweth, Rundle & Blows, 2009;Uesugi et al., 2017) research on the potential instability of the G-matrix.
Interestingly, we found a relationship between phenotypic evolution and change in the genetic covariance structure. The component under selection, Z11-16:OAc, was decoupled from the other components as observed by divergence in the loadings on gPC1 and gPC2. This observation thus fits the prediction that genetic variances and covariances may be reoriented to align the phenotype with the dominant direction of selection (Melo & Marroig, 2014).
These changes in the G matrix can likely facilitate a more specific response where only Z11-16:OAc evolves in the absence of indirect selection responses. The observed evolution of the G matrix is biologically relevant because selection for heterospecific mating avoidance favors higher rates of acetates in the female sex pheromone of H. subflexa (Groot et al., 2006) while intra-specific (sexual) selection is expected to act stabilizing on the rates of the other major and critical secondary sex pheromone components.

| CON CLUS ION
From our results, we conclude that: (1) pheromone components in H. subflexa can evolve in response to selection independent of other components of the sex pheromone; (2) selection need not erode the genetic variance in order to drive this phenotypic change; and (3) selection alters the genetic correlations among pheromone components. Our study thus shows that univariate selection responses in multicomponent sexual signals are possible, even though genetic correlations were high both prior to and during selection. Our study also shows that sexual signal components can respond to selection without reductions in genetic variance, but with rapid changes in the genetic covariance structure. These results correspond to geographic variation in the sexual signal found in this species and help to explain how sexual signals in general can respond to selection without running into constraints from indirect selection responses or depleted genetic variation. This may be one reason why sexual signals have evolved rapidly and repeatedly throughout the evolutionary history of animal taxa, thereby contributing to the speciation process.

ACK N OWLED G EM ENTS
We thank Dennis van Veldhuizen for his invaluable contribution to the rearing of our moths. We also thank Laura Caton, Frederica Lotito, and Lilian Seip for their help with the selection lines. We further thank Jet ten Berge for digitizing the pedigree, Michiel van Wijk for commenting on a previous version of the manuscript, and all members of the Groot lab for their helpful and insightful discussions.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data and analysis scripts are archived in the Dryad repository https://doi.org/10.5061/dryad.79cnp 5hxt.