Male‐biased gene expression resolves sexual conflict through the evolution of sex‐specific genetic architecture

Abstract Many genes are subject to contradictory selection pressures in males and females, and balancing selection resulting from sexual conflict has the potential to substantially increase standing genetic diversity in populations and thereby act as an important force in adaptation. However, the underlying causes of sexual conflict, and the potential for resolution, remains hotly debated. Using transcriptome‐resequencing data from male and female guppies, we use a novel approach, combining patterns of genetic diversity and intersexual divergence in allele frequency, to distinguish the different scenarios that give rise to sexual conflict, and how this conflict may be resolved through regulatory evolution. We show that reproductive fitness is the main source of sexual conflict, and this is resolved via the evolution of male‐biased expression. Furthermore, resolution of sexual conflict produces significant differences in genetic architecture between males and females, which in turn lead to specific alleles influencing sex‐specific viability. Together, our findings suggest an important role for sexual conflict in shaping broad patterns of genome diversity, and show that regulatory evolution is a rapid and efficient route to the resolution of conflict.

ngsTools 8 . The maximum root-mean-square deviation between iterations to assume convergence was 0.001.

Test for bias in estimating summary statistics due to unequal coverage between sexes
To test whether uneven depth between males and females, due to sex-biased expression, biases the estimation of per-gene Tajima's D or F ST between sexes, we simulated genomes for two populations using ms 11 . Specifically, we simulated two populations of equal and constant size diverging without gene flow 20, 200 and 2,000 generations ago. The same filtering criteria used for the F ST and Tajima's D analysis was imposed. We assumed an effective population size of 1000 and mutation and recombination rates of 5*10^(-8) per base per generation 6 . Under this model we simulated genes of 25kbp each, and a pseudochromosome of 10Mbp. We simulated 10,000 genes for each scenario. Sample sizes of the populations were set to 11 and 4 diploid individuals, representing male and female samples in our study respectively. From each simulated gene and chromosome, we simulated sequencing reads at the polymorphic sites using msToGlf from ANGSD 9 . We simulated three different sequencing scenarios with both sexes at high mean depth (30X), males at low depth (5X) and females at high depth (30X), males at high depth (30X) and females at low depth (5X). We then used the previously described pipeline to estimate F ST between sexes and Tajima's D for the whole population. For each scenario, Kendall's tau coefficient was used to assess the correlation between estimates of summary statistics. Kendall's tau measures association using a rank-based method accounting for ties.

Power assessment for detecting outliers of summary statistics due to unequal and low sample size between sexes
We simulated data with ms 11 to assess the power to detect outlier genes based on F ST and Tajima's D in case of high and even sample size compared to our data set with unequal and low sample size. For consistency with previous analyses, we assumed an ancestral effective population size of 1,000, mutation and recombination rates of 5*10^(-8) per base per generation, and an average gene length of 25kbp. We simulated sets of 1,000 genes along with a pseudo-chromosome of 10Mbp used to derive prior information.
In the case of F ST , we simulated 950 genes with a divergence time between the two populations of 40 generations. We then simulated 50 genes with no population structure and assessed the power to identify these outliers within the lowest 5 th percentile of the F ST distribution. Similarly but separately, we simulated 50 genes with a deeper divergence time (400 generations ago) and assessed the power to identify these outliers within the highest 95 th percentile of the F ST distribution.
In the case of Tajima's D, we simulated 950 genes with an effective population size of 10,000. We then simulated 50 genes with an effective population size of 100,000 and assessed the power to identify these outliers within the lowest 5 th percentile of the Tajima's D distribution. Similarly but separately, we simulated 50 genes with a lower effective population size (1,000) and assessed the power to identify these outliers within the highest 95 th percentile of the TD distribution.
We use the Matthews correlation coefficient 13  This metric is related to the classic F1-score but it additionally takes false negative rates into account.
We calculated Matthews correlation coefficient for detecting outliers in F ST and Tajima's D distribution for both large and even sample sizes (40 individuals in total) and small and uneven sample sizes (4 and 11 individuals per population).

Test for bias in identifying sex-biased genes due to unequal number of samples between sexes
Male guppies show a remarkable variety of colouration patterns in the wild and our male samples exhibit high phenotypic diversity. Therefore, we predict higher variability in expression across males than females, which can in turn influence the identification of sexbiased genes. We chose to use more male individuals (eleven samples) than females (four samples) to mitigate this and ensure the full transcriptional variation seen in males was reflected in this study. However, to exclude the possibility that the unequal number of male and female samples has biased the identification of sex-biased genes, and therefore measures of Tajima's D or F ST across sets of genes, we conducted pairwise rank order correlations of autosomal expression across samples. We calculated pairwise Spearman's rank correlation across female samples, and calculated average rho and p-value using the R 12 package Hmisc. 95% confidence intervals were calculated using bootstrapping with 1000 replicates. We repeated this for male samples. Finally, we sampled four males at random, 1000 times, and calculated average rho.

SUPPLEMENTARY RESULTS
We performed population genetic analyses using RNA-seq data and the guppy genome that we previously assembled and annotated 1 . Briefly, we mapped RNA-seq data from male and female tails (S1 Table) 14 in order to estimate allele frequencies and genotype likelihoods using ANGSD 9 and ngsTools 8 . ANGSD uses an empirical Bayes approach to account for the uncertainty in genotype calling from next-generation sequencing data when estimating neutrality test statistics. Estimates of site frequency spectra from next-generation data can be heavily biased by differences in coverage in next-generation sequencing data and ANGSD circumvents this problem by directly analyzing genotype likelihoods 9 . This approach has obvious benefits for our analyses given the large variance in coverage associated with expression data.

Population structure
We performed a Principal Component Analysis on the genotype likelihoods of male and female guppies to investigate population structure. In particular, population structure between male and female individuals could confound our estimation of Tajima's D and associated statistics. The first principal component explains 9.47% of the variance across individuals and the second component 8.8% (S1 Figure). We repeated the analysis after excluding sites on the sex chromosomes, which are subject to distinct evolutionary pressures relative to the autosomes 15 , and found similar patterns (S1 Figure). We conclude there is no genetic population structure between males and females in our analysis, consistent with a single large, outbred laboratory population We estimated Tajima's D for all genes within the guppy genome. For autosomal genes with Poecilia formsa orthologs (9582 genes, Sign test p<0.001, median=0.609), median Tajima's D is significantly greater than 0 (S2 Figure, Panel A), indicating a population contraction, consistent with initial collections used to establish the captive laboratory population.

Inbreeding coefficients
We also calculated inbreeding coefficients across our population using unlinked sites. Mean inbreeding coefficient across all individuals is 0.074 and we therefore assume Hardy Weinberg equilibrium for subsequent analyses. To further test for deviation from Hardy Weinberg Equilibrium, we calculated the likelihood of all samples being inbred, with inbreeding coefficients varying from 0.1 to 0.2. Our null hypothesis is that each sample has its estimated inbreeding coefficient while the alternate hypothesis is that all samples are highly inbred. A likelihood ratio test largely fails to reject the null hypothesis (p-value < 0.001, chi-square test, d.f=15) suggesting that all samples are either non-inbred or exhibit a negligible inbreeding coefficient.

Test for bias in identifying sex-biased genes due to unequal number of samples between sexes
We found significant (p<0.001) Spearman's rank correlation across all pairwise female comparisons with an average rho of 0.935 (95% CI = 0.936-0.981). We found significant (p<0.001) Spearman's rank correlation across all pairwise male comparisons with an average rho of 0.919 (95% CI = 0.917-0.939). Although marginally non-significant, the lower correlation across male samples relative to females is consistent with higher phenotypic diversity in the tail due to a high variation in colour patterns. Furthermore, when we resampled males with 1000 repetitions, we found a similarly high correlation across four randomly sampled individuals, with an average rho of 0.919. Together, this suggests that our statistical power to identify sex-biased genes is not limited by the unbalanced sample number. S1 Table: RNA-seq information for each sample 1 . Female-biased genes are defined as genes with log 2 fold change < -1 and significant p-value, male-biased genes are defined as genes with log 2 fold change > 1 and significant p-value.

Raw paired reads
High Tajima's D was defined as > 0.893 and low Tajima's D was defined as < 0.272 to account for the inferred population contraction within our population (Supporting Results).
High F ST was defined as > 0.047 and low F ST was defined as < -0.008 (Supporting Results). Test for Tajima's D bias due to unequal coverage between sexes. Upper panel. Correlation (left panels) between Tajima's D values calculated from even high coverage data and values estimated from uneven coverage, under a scenario of recent divergence and male (upper panel) or female (lower panels) samples at low depth. Histograms of the two distributions are presented on the right panels.