Experimental evolution confirms signatures of sexual selection in genomic divergence

Comparative genomics has contributed to the growing evidence that sexual selection is an important component of evolutionary divergence and speciation. Divergence by sexual selection is implicated in faster rates of divergence of the X chromosome and of genes thought to underlie sexually selected traits, including genes that are sex-biased in expression. However, accurately inferring the relative importance of complex and interacting forms of natural selection, demography and neutral processes which occurred in the evolutionary past is challenging. Experimental evolution provides an opportunity to apply controlled treatments for multiple generations and examine the consequences for genomic divergence. Here we altered sexual selection intensity, elevating sexual selection in polyandrous lines and eliminating it in monogamous lines, and examined patterns of divergence in the genome of Drosophila pseudoobscura after more than 160 generations of experimental evolution. Divergence is not uniform across the genome but concentrated in “islands”, many of which contain candidate genes implicated in mating behaviours and other sexually selected phenotypes. These are more often seen on the X chromosome, which overall shows divergence above neutral expectations. There are also characteristic signatures of selection seen in these regions, with lower diversity and greater Fst on the X chromosome than the autosomes, and differences in diversity on the autosomes between selection regimes. Reduced Tajima’s D implies that selective sweeps have occurred within the divergent regions, despite considerable recombination. These changes are associated with both differential gene expression between the lines and sex-biased gene expression within the lines. Our results are very similar to those thought to implicate sexual selection in divergence in natural populations, and hence provide experimental confirmation of the likely role of sexual selection in driving such types of genetic divergence, but also illustrate how variable outcomes can be for different genomic regions. Impact Summary How does sexual selection contribute to the divergence of genomes? It is often thought that sexual selection is a potent force in evolutionary divergence, but finding ‘signatures’ of sexual selection in the genome is not straight-forward, and has been quite controversial recently. Here we used experimental evolution to allow replicate populations of fruit fly to evolve under relaxed or strengthened sexual selection for over 150 generations, then sequenced their genomes to see how they had diverged. The feature we find are very similar to those reported in populations of natural species thought to be under strong sexual selection. We found that genomic divergence was concentrated in small patches of the genome rather than widespread. These are more often seen on the X chromosome, which overall shows greatly elevated divergence. There are also characteristic signatures of selection seen in these regions, with lower genetic diversity and greater differences on the X chromosome than the autosomes. Selection was probably strong in these regions. The changes are associated with both differential gene expression between the lines and sex-biased gene expression within the lines. Many of the patches of divergence also contain candidate genes implicated in mating behaviours and other sexually selected phenotypes. Our results provide experimental confirmation of the likely role of sexual selection in driving such types of genetic divergence.

genes implicated in mating behaviours and other sexually selected phenotypes. These are more often seen on the X chromosome, which overall shows divergence above neutral expectations. There are also characteristic signatures of selection seen in these regions, with lower diversity and greater F st on the X chromosome than the autosomes, and differences in diversity on the autosomes between selection regimes. Reduced Tajima's D implies that selective sweeps have occurred within the divergent regions, despite considerable recombination. These changes are associated with both differential gene expression between the lines and sex-biased gene expression within the lines. Our results are very similar to those thought to implicate sexual selection in divergence in natural populations, and hence provide experimental confirmation of the likely role of sexual selection in driving such types of genetic divergence, but also illustrate how variable outcomes can be for different genomic regions.

Impact Summary
How does sexual selection contribute to the divergence of genomes? It is often thought that sexual selection is a potent force in evolutionary divergence, but finding 'signatures' of sexual selection in the genome is not straight-forward, and has been quite controversial recently. Here we used experimental evolution to allow replicate populations of fruit fly to evolve under relaxed or strengthened sexual selection for over 150 generations, then sequenced their genomes to see how they had diverged. The feature we find are very similar to those reported in populations of natural species thought to be under strong sexual selection.
We found that genomic divergence was concentrated in small patches of the genome rather than widespread. These are more often seen on the X chromosome, which overall shows greatly elevated divergence. There are also characteristic signatures of selection seen in these regions, with lower genetic diversity and greater differences on the X chromosome than the autosomes. Selection was probably strong in these regions. The changes are associated with both differential gene expression between the lines and sex-biased gene expression within the lines. Many of the patches of divergence also contain candidate genes implicated in mating behaviours and other sexually selected phenotypes. Our results provide experimental confirmation of the likely role of sexual selection in driving such types of genetic divergence.

Introduction
The role of sexual selection in influencing evolutionary divergence and speciation is unclear (Panhuis et al., 2001;Ritchie, 2007;Maan & Seehausen, 2011;Servedio & Boughman, 2017). Associations between species diversity and proxies of sexual selection such as sexual dimorphism or mating system variation often imply that sexual selection can accelerate divergence, especially when acting alongside natural selection (Arnqvist et al., 2000;Gage et al., 2002;Ellis & Oakley, 2016). However, different indicators of sexual selection give contrasting results in such comparative studies, and a consensus is not clear (Kraaijeveld et al., 2011;Janicke et al., 2018). One potentially compelling source of evidence that sexual selection is involved in divergence is coming from the increasing number of comparative genomic studies available across a range of organisms. Many descriptions of genomes, including those of species thought to have undergone strong sexual selection such as the Hawaiian Drosophila or African cichlids, have found that genes associated with mating behaviour or sensory perception potentially involved in sexual communication are often outliers in measures of divergence (e.g. Mattersdorfer et al., 2012;Kang et al., 2016). It has also been known for some time that genes which diverge particularly rapidly and show stronger signatures of positive divergent selection are often sex-biased in expression (Pröschel et al., 2006;Ellegren & Parsch, 2007;Zhang et al., 2007). Sex-biased gene expression itself, especially malebiased expression, evolves rapidly and is associated with indicators of sexual selection such as increased sexual dimorphism in birds (Harrison et al., 2015;Wright et al., 2019). Additionally, divergence of sex chromosomes between species is usually much greater than autosomes, sometimes dramatically so (Counterman et al., 2004;Ellegren et al., 2012).
However, such patterns of divergence are not necessarily driven by elevated sexual selection on these genes or genomic regions. The evolution of sex-biased gene expression is often due to sexually antagonistic selection on gene expression, which is often an important factor in sexual selection. However, it is complicated by additional factors including dosage compensation, as well as turnover of sex-biased expression and resolution of conflict via sexlinkage or sex-limited expression (Mank et al., 2010a;Wright et al., 2019). The increased divergence of sex chromosomes is potentially influenced by many factors, only one of which is sexual selection. These include a greater role of genetic drift due to a smaller effective population size on X chromosomes compared to autosomes, dominance effects, and other consequences of sex-linkage such as dosage compensation (Vicoso & Charlesworth, 2006;Ellegren, 2009;Mank et al., 2010b). Hemizygosity results in a lower effective population size (N e ) on the X (N eX ) than on autosomes (N eA ). Under random mating the ratio is expected to be 3:4 and this should reduce neutral diversity and increase between-species divergence by the same proportion (Vicoso & Charlesworth, 2006). Hemizygosity should also result in an increased efficacy of selection for partially recessive beneficial mutations and against recessive deleterious mutations on the X. Finally, because of the inheritance patterns of X-linked loci (males transmit them only to daughters while females transmit them to both daughters and sons), sex-specific selection as well as sexual selection will influence their divergence (Mank et al., 2010a;Corl & Ellegren, 2012).
It is very difficult to infer the historical role of different evolutionary processes from patterns of contemporary divergence between populations and species, because they can result in similar genomic signals (Butlin et al., 2012). One way of directly addressing the role of sexual selection or mating system variation in genomic divergence is to examine the genomic consequences of experimental evolution under manipulated sexual selection regimes in the laboratory. A great advantage of this approach is that there are potentially fewer confounding variables involved than when making comparisons across species or natural populations.
However, a disadvantage is that the time scale over which divergence can be studied is typically much shorter than evolutionary time-scales in nature. Studies of experimental evolution and speciation are in their infancy, and general conclusions are, as yet, difficult to draw (White et al., 2020). Enforcing monogamy in otherwise polyandrous species will lead to both changes in the intensity of sexual selection and the balance of sexual conflict, as it effectively eliminates sexual selection and sexually antagonistic selection. A classic example of such manipulation is where D. melanogaster were kept under enforced monogamy for about 50 generations (Holland & Rice, 1999). Females from the monogamy treatment had reduced longevity compared to ancestral females, when exposed to ancestral males. This was expected because the reduction of conflict should favour less harmful males and females that are less resistant to male harm. Other experimental evolution studies under altered mating systems have been performed in dung flies Martin & Hosken, 2003), different species of fruit flies (D. melanogaster; (Gerrard et al., 2013;Hollis et al., 2014;Innocenti et al., 2014;Perry et al., 2016); D. pseudoobscura; (Crudgington et al., 2005), seed beetles (McNamara et al., 2020) and hermaphroditic flatworms (Janicke et al., 2016). Though aspects of the treatments differ amongst such experiments, some common patterns have emerged. Gene expression changes are seen, especially of genes that are initially sex-biased, though the details can vary between studies (Hollis et al., 2014;Veltsos et al., 2017). Moreover, gene expression changes can be more pronounced for genes expressed in reproductive tissues (Innocenti et al., 2014) , and genes involved in the post-mating physiological manipulation of female egg-laying and re-mating rates (Perry et al., 2016).
A feature emerging from genomic comparisons between diverging species is that details of genomic architecture complicate the assessment of patterns of divergence across chromosomes. Whole chromosomal regions can show correlated responses due to reduced recombination and hitchhiking effects, especially in species with segregating inversions. Early studies of species differences interpreted "islands" of divergence in the genome as resulting from divergent selection on genes within these regions with gene flow homogenising the genetic background (Turner et al., 2005;Nosil et al., 2009). More recently it has been appreciated that chromosomal inversions and other regions of low recombination or diversity can accentuate such clustered divergence (Noor & Bennett, 2009;Cruickshank & Hahn, 2014;Wolf & Ellegren, 2016;Ravinet et al., 2017). "Barrier loci", genomic regions under divergent selection that restrict gene flow (Butlin & Smadja, 2018), may occur within such clusters but the lack of recombination makes them difficult to localise precisely. In experimental evolution the amount of recombination will be determined by both genomic architecture and the number of generations completed during the study, which is often modest in studies of eukaryotes.
Also, in experimental evolution the lines can be kept effectively allopatric, so homogenising gene flow in regions not experiencing selection should be absent. The genomic divergence which occurs during experimental evolution is usually extensive, with widespread differences dispersed throughout the genome (Kawecki et al., 2012;Tobler et al., 2014;Michalak et al.,

2019).
Here we directly test the influence of sexual selection on genomic divergence. We examine replicated experimentally evolved lines of D. pseudoobscura in which sexual selection has been manipulated for over 160 generations. One set of 4 replicate lines were raised under enforced monogamy (M lines), which should eliminate both sexual selection and conflict.
Another 4 replicates were reared under elevated polyandry (E lines), with 6 males per female.
Polyandry mediates the strength of both intra-and intersexual selection and sexual conflict (Pizzari & Wedell, 2013) and elevated polyandry will increase both pre-and post-copulatory sexual selection via female choice and sperm competition (Snook, 2014). Previous studies of these lines have found divergence in some, but not all, of the types of traits predicted to diverge under sexual selection. Sperm morphology and heteromorphism, and testis mass did not diverge, but E males had larger accessory glands and a greater mating capacity (Crudgington et al., 2009), were more competitive in mating encounters (Debelle et al., 2016), and produced more attractive courtship song than M males (Debelle et al., 2017). Coe- volutionary changes have occurred in female song preferences (Debelle et al., 2014 ). Sexually dimorphic cuticular hydrocarbons have also diverged between the lines (Hunt et al., 2012).
Patterns of gene expression have changed between the lines. E females show an increase in expression of genes normally enriched in ovaries (Immonen et al., 2014). Sexbiased genes responded more strongly to the sexual selection treatment, but the direction of gene expression changes differed between sexes, tissues, and according to courtship experience (Veltsos et al., 2017). In most cases, the transcriptome was "feminised" under polyandry (i.e. female-biased genes were up-regulated or male-biased downregulated in E lines), in a striking contrast to a similar study with D. melanogaster (Hollis et al., 2014). Males changed in patterns of gene expression in the testes and accessory glands, and changes in gene expression in females following mating also diverged, especially in the female reproductive tract (Veltsos et al. in prep.).
Here we examine genomic divergence between these lines using a pool-sequence approach after more than 160 generations of experimental evolution. The relatively long timescale of this study should reduce linkage effects on allele frequency changes. We identify alleles that changed in frequency consistently across the replicates, to help reduce the potentially confounding effects of drift or replicate-specific selection. We find that divergent SNPs are not distributed randomly across the genome, but occur in distinct, obvious clusters. We examine what genes are involved and find several with mutant phenotypes related to mating and courtship behaviours. We found that the X chromosome has accumulated more divergence than the autosomes and explore if divergence is associated with recombination rate or changes in gene expression between the experimental lines.

Experimental Evolution
A full description of the experimental evolution methods is available elsewhere (Crudgington et al., 2005). Briefly, a population of D. pseudoobscura was established from 50 wild caught females, bred in the laboratory for four years then four independent monogamy (M) and elevated polyandry (E) lines were established. M females were housed with a single male and E females with 6 males, with females typically mating with two or three males. The effective population size was maintained around 120   sex ratio, thus reflecting the differential offspring production across families (Crudgington et al., 2005;Crudgington et al., 2009). Enforced monogamy is expected to eliminate sexual selection and sexual conflict while elevated polyandry increases both pre-and postmating sexual selection and sexual conflict (Crudgington et al., 2005;Bacigalupe et al., 2007;Crudgington et al., 2009).

Sequencing and Mapping
Sequencing was carried out after ca. 160 generations of selection (specifically, 164 for replicate 1, 163 for replicate 2, 162 for replicate 3, and generation 160 for replicate 4). Two were filtered to remove duplicate reads, reads with a mapping quality < 30, and any reads which were not properly paired, using samtools (v 1.3; Li et al., 2009following Schlotterer et al., 2014. Reads were locally re-aligned around indels using GATK (v3.7.0;McKenna et al., 2010;DePristo et al., 2011). The .bam files for each line were then merged using bamtools (Barnett et al., 2011) and the genome-wide coverage calculated from these merged files with bedtools (v. 2.26; Quinlan & Hall, 2010). SNPs were called using a heuristic SNP calling algorithm (PoolSNP; Kapun et al., 2020). Sites were considered only if the total coverage at the site was > 17 and < the 95 th percentile for each contig or chromosome. An allele was only called if the count for that allele across all pools was > 16. Nearly 2 million SNPs were called and used in downstream analyses (see Supplementary Material).

Genomic Analyses Identifying Consistent Allele Frequency Differences
Many evolve and resequence studies of Drosophila find that a multitude of SNPs have diverged, perhaps tens of thousands (Michalak et al., 2019). The number is inflated upwards at least in part due to segregating inversions and other areas of low recombination, and hitchiking (Barghi & Schlotterer, 2019). In order to focus on the loci most likely to have diverged due to the treatment, we only considered as significant SNPs which diverged consistently across all 4 replicate pairs of lines. We identified these using quasibinomial Generalised Linear Models, which are less prone than other statistical approaches to be influenced by strong divergence in only some replicates (Wiberg et al., 2017). The model structure applied was; where y is the allele frequency, treatment is the experimental evolution treatment regime, and e is a quasibinomially distributed error term. +1 was added to zero counts in any population. P-values were converted to q-values using the "qvalues" R package (v. 2.16.0; Storey & Tibshirani, 2003). A threshold of 0.05 was chosen to control the false discovery rate (FDR), thus we define "top SNPs" as those with q-value < 0.05 and the rest are referred to as "background" SNPs.

Genetic Diversity and Differentiation
We calculated genome-wide genetic diversity statistics (π and Tajima's D) for windows of 50kb (with a 10kb overlap) using available python scripts (Kapun et al., 2020). Similarly, we computed pairwise F st estimates between E and M line pairs for each SNP using the R package "poolfstat" (v. 0.0.1; Hivert et al., 2018), averaged in windows of 50kb (with a 10kb overlap between windows). Comparisons of parameters between selection regimes and genomic regions were tested using non-parametric Wilcoxon tests. Additionally, we estimated the F st expected from drift and differences in effective population sizes on X chromosomes (F X ) as in  using the equations of (Ramachandran et al., 2004) (equation 8), F X is given by: where, z is the ratio of the number of breeding males to females and F A is the observed F st on autosomes. We assumed z to be either 1 or 6 to represent extreme possibilities based on the

Linkage Disequilibrium (LDx)
Although haplotype information is not available from pool-seq data, short range linkage information is available from paired reads. We used LDx (Feder et al., 2012) to first compute the r 2 of SNPs located on the same read pairs. We only used SNPs with a minor allele frequency > 0.1, a minimum coverage of 10, a maximum read coverage of 400, and a phred score > 20. We binned pairs of SNPs into distance classes and then computed mean r 2 per distance class. We only used distance classes with a minimum of 5 SNPs. We estimated the decay of r 2 as a function of distance by fitting a linear model of r 2 as a function of the log of the distance between the SNPs. Thus, the slope measures the decay rate of linkage due to recombination (Feder et al., 2012), giving an indication of the distance over which LD is present. In regions of low recombination one would expect high overall values of r 2 but a weakly negative slope as LD is maintained over relatively longer regions of the genome.
Comparing the slope parameter across different genomic regions gives an indication of differences in the recombination rate (or extent of selective sweeps). This was performed for each chromosome, as well as for different regions on the 3 rd chromosome (see below).

Functional Genomics
In analyses of gene function and regulation we used the D. pseudoobscura annotation and a dataset of regulatory long non-coding RNAs (lncRNAs; Nyberg & Machado, 2016).
We identified genes or lncRNAs within a distance of 10kb up-or downstream of focal SNPs with bedtools (Quinlan & Hall, 2010) intersect (keeping any potential ties). Enhancer regions, transcription factor binding sites, and other regulatory regions can occur up to 1 Mb up-or downstream from a target gene in other species (e.g. Maston et al., 2006;Chan et al., 2010;Werner et al., 2010;Pennacchio et al., 2013) but typically lie within 2kb of a gene region in D. melanogaster (Arnosti, 2003), 10kb thus represents a compromise. We submitted the implicated genes to ModPhEA (Weng & Liao, 2017) for phenotypic enrichment analysis. We combined the phenotypic classes "courtship behavior defective" (FBcv:0000399) and "mating rhythm defective" (FBcv:0000401) into one phenotype group and also tested the phenotypic class "stress response defective" (FBcv:0000408) for enrichment. We chose these classes a priori because they were most likely to be involved in phenotypic differences between the treatments related to mating or courtship behaviour and responses. We also took advantage of recent sex-and tissue-specific expression data from the same experimental evolution lines (Veltsos et al., 2017;Veltsos et al., in prep.) to ask if SNPs co-localised with genes that are differentially expressed between the lines and if these also show different levels of diversity (Tajima's D) or differentiation (F st ) between E and M lines.
Full details of the differential expression analyses are given elsewhere (Veltsos et al., 2017 Veltsos et al., in prep.). For each tissue (head and abdomen) we performed an ANCOVA with chromosome (autosome, X-chromosome right arm, and X-chromosome left arm) as a co-factor, as well as mean Tajima's D across E lines and mean Tajima's D across M lines as co-variates. We also included the interactions between Tajima's D and chromosome. The full model is:

ΔSB EM ~ chromosome + TajD E + TajD M + TajD E :chromosome + TajD M :chromosome
We further extracted the 30bp up-and down-stream of each SNP from the reference genome using gffread from the Cufflinks package (v2.2.1; (Trapnell et al., 2010) (Kofler & Schlotterer, 2012). We considered SNPs to be associated with genes if they occurred within 10kb up or downstream of an annotated gene.
An empirical p-value distribution was produced from 1 million simulated SNP sets.
All statistical analyses were made with R (v. 3.6.3; R Development Core Team 2020) except where otherwise stated. Figures were drawn using the "ggplot2" package (v. 2.2.1; (Wickham, 2009) and associated packages (table S1).

Consistent Allele Frequency Differences
In total, 480 SNPs show significant consistent allele frequency differences due to the experimental evolution treatment (hereafter the "top SNPs"). These occur on all of the main chromosomes but many show striking co-occurrence into a few clusters of highly

Linkage Disequilibrium
A combination of hitchhiking and selective sweeps could lead to clustered genomic divergence, often with low diversity, especially in regions of low recombination such as telomeric regions. We examined patterns of linkage disequilibrium in the clusters and if this varied with treatment. Throughout the genome, the decay rate (a parameter) of LD is generally shallower (i.e. less negative) in the E treatment ( figure 4A). This is seen for chromosome 3 as well as both arms of the X chromosome ( figure 4A). A lower decay rate is indicative of more LD, due to less recombination and/or a potential for greater hitchhiking under positive selection. Contrary to predictions, we found a steeper rate of decay (less LD) within the differentiated region of chromosome 3 than outside it, especially in E lines ( figure   4B and C). Although statistically significant (F (2,13) = 4.6, p < 0.001), these differences are slight. The most striking pattern overall is greater overall LD on chromosome 3.

Gene functions and expression variation
The top SNPs were not significantly enriched in any GO terms after correcting for multiple testing, even at a 10% FDR (table S4). Similarly, there was no enrichment of genes with annotations for mating behaviour or stress response phenotypic classes. However, several genes within 10kb of a top SNP are potentially interesting candidate genes for traits evolving under sexual selection based on described functions (table S4). For example, the genes Odorant-binding protein 47a (Obp47), pickpocket 6 (ppk6), and Accessory gland protein 53C14c (Acp53C14c) all occur within 10kb of a top SNP and are genes potentially underlying sexually selected behaviours or traits. Two of these genes (ACP53C14c and Obp47a) are within the region of highly differentiated SNPs on the 3 rd chromosomes, which also includes several additional accessory gland proteins (Acp53Ea, Acp53C14b, Acp53C14a), and other genes (table S4), all of which are thought to influence mating and courtship behaviours or phenotypes based on known functions of similar genes in D.

melanogaster.
Previous studies have shown that there is divergence in gene expression patterns between E and M lines (Immonen et al., 2014;Veltsos et al., 2017;Veltsos et al., in prep.).
We therefore asked if these expression differences were associated with the top SNPs. Genes The regions immediately up-or down-stream of top SNPs are not enriched for TF binding motifs or lncRNAs, after correction for multiple testing, so there were no obvious differences between treatments in regions governing gene expression.

Discussion
There is much debate about the influence of sexual selection and sexually antagonistic selection on patterns of genomic variation (Mank, 2017;Sayadi et al., 2019) and how this may influence divergence between species (Wolf & Ellegren, 2016). Sex-biased gene expression, especially male-bias, evolves quickly and is related to phenotypic sexual dimorphism (Wright et al., 2019). Outliers in genome scans often implicate sexual selection as a diversifying force (Andres et al., 2008;Blankers et al., 2018). Sexual antagonism may be associated with genomic signatures of selective sweeps or balancing selection (Cheng & Kirkpatrick, 2016;Wright et al., 2019) and may be promoted by strong sexual selection.
However, inferences of the sources of selection on natural variation in genomic divergence are usually indirect and ambiguous, because multiple forces act in concert to produce variation seen at the genomic level in nature. Here we used experimental evolution to alter sexual selection intensity, elevating sexual selection in polyandrous lines and eliminating it in monogamous lines, and examined patterns of divergence in the genome after more than 160 generations of experimental evolution.
Many of the results we found recapitulate patterns seen in natural populations and between species. Divergence is not uniform across the genome but clustered in "islands" of divergence, some of which contain candidate genes for an involvement in mating success.
These clusters are more often seen on the X chromosome, which is a "hot spot" for divergence. There are signatures of selection within the islands of divergence, with marginally lower diversity (π) within clusters than the rest of the genome, but only in M lines.
Moreover, F st between E and M lines is greater within clusters. F st is also greater on the X than autosomes, and differences in diversity are seen in the autosomes between selection regimes. Tajima's D implies selective sweeps have occurred, but only within some of the divergent regions. These patterns of diversity and divergence are associated with changes in both differential gene expression between the lines and sex-biased genes. Overall, F st between the lines is high in all replicates, probably due to low overall effective population sizes, though effective population sizes are similar between E and M lines .
The concept of "islands" of divergence originated from omparisons of genomic divergence between species has given rise to (Nosil et al., 2009;Ravinet et al., 2017). These are usually thought to have arisen due to the combination of strong selection on barrier loci and genetic hitchhiking within genomic regions, with background gene flow reducing divergence outside of the islands. Here we find distinct clustered divergence akin to the islands seen in natural systems. Our system is effectively allopatric, so there was no background gene flow counteracting divergence outside of these clusters, which therefore must have arisen due to strong localised divergent selection across all replicates. Although D.
pseudoobscura has relatively well-characterised inversion polymorphisms (Sturtevant & Dobzhansky, 1936;Dobzhansky & Sturtevant, 1938;Wallace et al., 2011), the clusters we describe do not correspond to the most common inversions known for this species, which are often very large. Our short-read sequencing approach allowed some examination of LD and there was no suggestion of reduced LD within the clusters. In fact, the large peak at the right end of chromosome 3 (figure 4) surprisingly seems to be within a region of high recombination (which is often suppressed at telomeric regions). Interestingly, recombination is higher within this peak than the chromosome-wide rate, but also differs between the treatments, being greater in the M lines. Perhaps selection against recombination was reduced in monogamous individuals because of epistatic interactions in the region which were important in sexual selection or sexual conflict. There was no obvious difference in LD in the other clusters but their smaller size and hence "noisier" estimates makes robust inferences from pool-sequence data difficult. Indeed, the estimates of LD within the cluster on chromosome 3 also rely on relatively few SNPs at longer ranges compared to the rest of the chromosome, so inferences should be taken with caution. categories, however, they do include strong candidate genes for an involvement in mating system evolution. For example, the large region on chromosome 3 contains numerous accessory gland proteins. In D. melanogaster these are well known to influence male reproductive success, exert antagonistic effects on female fecundity and lifespan, and play a role in sperm competitive success (Chapman et al., 1995;Ram & Wolfner, 2007). Some of the evolutionary response in E lines is antagonistic, because M females have a higher fecundity when mated with E males. Moreover, when mated to E males, the reproductive schedule of M females is manipulated to the males benefit (Crudgington et al., 2010).
Accessory gland proteins show accelerated coding sequence and expression evolution across species (Swanson & Vacquier, 2002;Begun & Lindfors, 2005). Other genes within the clusters are involved in sexual chemical communication, which is also often implicated in outlier analyses in genome comparisons between species (Smadja & Butlin, 2009). For example, mutants of members of the pickpocket family in D. melanogaster show aberrant male mating success because of their involvement in the detection of female pheromones (Thistle et al., 2012;Toda et al., 2012). E males, subject to both intra-and intersexual selection, have diverged in aspects of courtship behaviour, such as time until initiation of courtship, have a higher intensity courtship song and have a higher competitive mating success than M males (Debelle et al., 2016;Debelle et al., 2017).
If strong selection has driven this clustered genomic divergence, an interesting question is whether the responses to selection are stronger in the E or M lines. Imposing monogamy on a naturally polyandrous species probably leads to relaxed selection on many genes involved in intra-or intersexual competition. Therefore, the response is likely to involve changes in both the intensity and direction of selection on some loci. Thus, perhaps the variation in signals of selection we see in Tajima's D and changes in LD are to be expected. Overall, we see stronger reductions in divergence in E lines, perhaps suggesting that directional selection was stronger when sexual selection was strengthened.
One pattern very commonly seen in studies of natural populations and species is more rapid divergence of the X chromosome (Vicoso & Charlesworth, 2006). We also see this here, the X having a higher prevalence of divergent clustered regions and higher F st between the lines. Remarkably, all SNPs with fixed differences between the lines occurred on the X. Faster X evolution can occur for many reasons, including greater genetic drift due to its smaller effective population size. We calculated expected X/A divergence ratios under a range of plausible sex ratios and the observed X/A divergence exceeded all of them, suggesting the accelerated X divergence is not due to drift effects alone, and that either selection or a combination of effects are involved. Genes under sexual selection are potentially more likely to be sex-linked, due to antagonistic, or sex-specific selection, or if sexually selected loci show dominance effects (Reinhold, 1998;Kirkpatrick & Hall, 2004;Grieshop & Arnqvist, 2018).
Links between genomic parameters and gene expression variation have been a somewhat contentious source of evidence of sexual selection, especially antagonistic forms of sexual selection (Kasimatis et al., 2019;Cheng & Kirkpatrick, 2020;Mank et al., 2020).
Genes that are male-biased in expression show accelerated divergence between species and sex-biased gene expression shows rapid evolution and turnover (Pröschel et al., 2006;Harrison et al., 2015). Whether sex-biased expression is expected to be related to sex-specific F st or signatures of balancing selection such as Tajima's D is open to debate, partly because of the potential resolution of antagonistic selection by the strengthening of sex-biased expression. Previously we found that gene expression differences have evolved between the lines, especially in sex-biased genes (Veltsos et al., 2017). Here we show that there is significant overlap between these genes and genes and the regions of genomic divergence of the lines found here. Thus, the expression divergence is associated with the broad patterns of genomic divergence. Also, genetic differentiation (F st ) is greater for the differentially expressed genes, once again recapitulating patterns from natural systems (sex-biased genes here are not more likely to be sex-linked, so this is independent of the large X effect seen).
We find no general difference in Tajima's D between these DE loci. However, there is one very intriguing pattern where the magnitude of change in sex-biased gene expression is related to Tajima's D. As ΔSB increases (more male-biased expression in E lines) Tajima's D in these lines becomes more negative. This pattern is potentially consistent with more resolved sexual conflict in the M lines, because males in M lines are released from sexual selection, and selection driving female-beneficial alleles to fixation or high frequency could result in sweeps and/or reduced balancing selection.
In conclusion, we have examined genomic divergence following >160 generations of experimental evolution under altered mating systems. We find that genomic divergence between the experimental lines is highly clustered in the genome, much greater on the X and is associated with changes in gene expression between the experimental lines. Associations

Author contributions
RAWW performed the data analysis. PV contributed data. The experiment was designed by MGR and RRS. All authors contributed to writing the MS

Data Accessibility
Raw reads have been deposited in the short read archive (SRA) of NCBI under the BioProject   B) The X chromosome to autosome ratio of π in the replicates of E and M lines and overall.