Introduction

The evolutionary role of sexual reproduction is the generation of variation between individuals. In diploid organisms this is achieved through recombination during meiosis and the generation of haploid gametes. Traditionally, the type of the gametes defines whether an individual is called a biological female or a male. This binary distinction is a hallmark of most sexually reproducing species and it is therefore often considered as a corner stone of evolutionary conservation. However, already Darwin has pointed out that the phenotypes of the sexes should evolve fast, when they are subject to sexual selection, e.g. play a role in mate choice (Darwin, 1871). In some cases, such as the famous peacock’s tail, sexual selection may even override optimal adaptation to the environment. But also the maintenance of an equal sex-ratio should be subject to continuous evolution, due to the fact that an excess of one sex in a population would generate a selection effect towards favoring the rarer sex (Hamilton, 1967).

The sexual phenotypes of adult individuals are generated by hundreds to thousands of genes with sex-biased expression (Mank & Rideout, 2021; Parsch & Ellegren, 2013; Tosto, Beasley, Wong, Mank, & Flanagan, 2023). The combination of the gene variants of these genes is expected to determine the range of phenotypes of individual males and females (Tosto et al., 2023). For example, human height is partially determined by genes with sex-biased expression, resulting in overall smaller females than males (Naqvi et al., 2019). However, height distributions of males and females are highly overlapping, such that knowledge of the height of a given individual would not be sufficient to classify it as either female or male. In fact, such overlapping distributions are characteristic for sex-related phenotypic characters (Maney, 2016).

Different organs are known to be shaped by different sets of sex-biased genes (Naqvi et al., 2019; Rodríguez-Montes et al., 2023), making the possible source of sex-specific variation between individuals even larger. Sex-biased gene expression is also known to evolve fast, supporting the notions sexual selection theory (Tosto et al., 2023). However, the evolutionary comparisons of sex-biased gene expression have so far focused on comparisons between relatively distant species and relatively few individuals per species only (Naqvi et al., 2019; Rodríguez-Montes et al., 2023), making it difficult to gauge the true evolutionary rates. Here we use the house mouse radiation to address the evolution of sex-biased genes and their expression. The system is an excellent model to study evolutionary turn-over patterns at a shallow evolutionary population level as well as between sub-species and species (Harr et al., 2016; Zhang, Xie, Ullrich, Zhang, & Tautz, 2021). For our study we use outbred individuals from two Mus musculus subspecies, M. m. domesticus (DOM) and M. m. musculus (MUS) as well as two sister species, M. spretus (SPR) and M. spicilegus (SPI). To avoid confusion with species status, we designate all four with the neutral term “taxon” and use only the three-letter abbreviations. The taxa cover an evolutionary distance of about 2 million years (Figure 1). Given the short reproduction cycles and the fast molecular evolution in rodents, this divergence corresponds approximately to the equivalent of gibbon to human evolution. Hence, the system allows to study the continuum from microevolutionary to macroevolutionary patterns.

Figure 1:

We use individuals derived from natural populations but kept under common environmental conditions and in an outbreeding regime to retain the natural variation (Harr et al., 2016). We have prepared gonadal and non-gonadal tissues from both sexes to quantify the variation and evolutionary turnover of sex-biased gene expression. It has recently been shown that the major onset of sex-biased gene expression is at the developmental transition to sexual maturation (Rodríguez-Montes et al., 2023) and we are therefore focusing on sexually mature stages. The use of a set of individuals from each taxon allows to compare within-group variation with divergence between groups, a contrast that has been lacking so far in most studies on sex-biased gene expression, except for humans (Khodursky et al., 2022).

We find thousands of sex-biased genes, but with marked differences in numbers between the non-gonadal organs in the different evolutionary lineages. Intriguingly, only a small percentage of them is conserved in their sex-biased expression patterns among all lineages and the variances of expression between individuals are higher than for non-sex-biased genes.

To represent the individual variation in expression of sex-biased genes, we have developed a sex-biased gene expression index (SBI) that combines the sex-biased gene expression for each individual to reflect its overall femaleness or maleness in each organ. In the non-gonadal organs, we find variable patterns of distinction and overlap between the sexes, while they are always well separated in the gonads.

We apply the SBI also to human data from the GTEx resource (Aguet et al., 2020) and find overall much fewer sex-biased genes when applying the same criteria as for the mouse data, and stronger overlaps between the sexes than in mice. Only very few genes show a conserved sex-biased expression across mice and humans, most only in a subset of the organs or cell types.

Our results show that sex-biased gene expression evolves even faster than thought so far, implying that sexual conflict is one of the major drivers of evolutionary divergence already at the early species divergence level. The different organs show a large individual variation in sex-biased gene expression, making it impossible to classify individuals in simple binary terms. Hence, the seemingly strong conservation of binary sex-states does not find an equivalent underpinning when one looks at the genetic makeup of the sexes.

Results

Gene expression data were collected from taxa and organs depicted in Figure 1. Nine age-matched adult females and males each were chosen among the outbred individuals from each of the four taxa and the same set of organs were retrieved for each of them for RNASeq analysis. We used five non-gonadal organs and three gonadal organ parts for each of the sexes (Figure 1), yielding a total of 576 samples (suppl. Table S1).

Figure 1 provides a graphic overview on the numbers of sex-biased genes identified in the different organs of the mouse taxa; detailed numbers are provided in suppl. Table S2. The numbers range between 13 - 5,645 for the non-gonadal organs and 3,444-13,383 for the gonadal pairwise comparisons (see Methods for the scheme of gonadal pairwise comparisons).

Female-biased genes outnumber male-biased genes in all comparisons. Among the non-gonadal organs, brain has the lowest numbers of genes with sex-biased expression, kidney the highest. Most notably, the numbers of sex-biased genes are subject to fast turnover between the taxa for the non-gondal organs. For example, DOM has an unusually high number of sex-biased genes in the heart, compared to the other taxa. On the other hand, the kidney bias is much stronger in DOM and MUS than in SPRE and SPI and liver has the highest number of female-biased genes in SPI. The gonadal organs have the largest numbers of sex-biased genes, but these generate evidently also very different tissues in the two sexes. The overall numbers remain relatively constant between the taxa.

Evolution of sex-biased gene expression

The evolutionary turnover of genes that become subject to sex-biased expression in any given taxon is extremely high, even between the closely related taxa that are studied here. Figure 2 shows which fractions of genes with sex-biased expression are shared between the taxa for each organ in any pairwise comparison. The turnover is particularly high for the non-gonadal tissues. In brain and heart, fewer than 10% of the genes occur in more than one taxon, for the other organs, fewer than half occur in more than one taxon. The fractions of shared genes are higher for the gonadal tissues (Figure 2), but even these tissues show a substantial turnover of sex-biased genes. Interestingly, many of the sex-biased genes reverse their state in at least one taxon from female-bias to male-bias or vice versa. This includes 596 genes for non-gonadal organs and 3,895 for the gonadal organs (suppl. Table S2). Hence, although the overall numbers of sex-biased genes to not change much between the taxa for the gonadal organs, the actual composition of the gene sets evolves substantially.

Figure 2:

Variance patterns

Given the fast turnover of genes with sex-biased expression, one can ask whether the variances of such genes are also high within each of the populations. As a measure of variance, we used the interquartile range IQR (difference between the 75th and 25th percentiles of the data, representing the non-parametric version of the coefficient of variation), divided by the median of the data range. This measure buffers against outlier values that can frequently be observed in the data. Sex-biased genes show indeed in most cases higher variances than non-biased genes, but again to different extents in the two sexes and in the different taxa (Figure 3A).

Figure 3:

Previous studies have suggested that protein sequences of male-biased genes are particularly fast evolving (Grath & Parsch, 2016; Harrison et al., 2015), while the opposite pattern was found for the avian brain (Mank, Hultin-Rosenberg, Axelsson, & Ellegren, 2007). To test for elevated positive selection in our population comparisons, we used the McDonald-Kreitman (MK) test that compares the level of polymorphism versus fixation of sites to infer the frequency of positive selection events (McDonald & Kreitman, 1991). However, its results can be influenced by confounding factors, such as the level of segregation of slightly deleterious variants, as well as population demography. We have therefore used here the asymptotic MK test that takes care of these problems (Haller & Messer, 2017; Messer & Petrov, 2013), but requires the testing across a large number of loci to become reasonably robust. We have therefore combined the whole set of sex-biased genes of the non-gonadal organs for a given taxon in the test. We restrict the comparison to the DOM and MUS samples, since these are closest to the reference sequence to ensure appropriate alignments of the coding regions. The polymorphism data were taken from previous whole-genome sequencing studies of these populations (Harr et al., 2016).

The results are shown in Figure 3B. Both, female-biased and male-biased genes show relatively higher alpha values compared to non-biased genes. To ask whether genes that are sex-biased in one taxon, but not sex-biased in another taxon evolve differently, we generated reciprocal sets of orthologous genes that show no sex-biased expression in the respective other taxon (DOM versus MUS). Intriguingly, the reciprocal values are indeed much lower (Figure 3B), implying that the enhanced positive selection pressure is triggered by their status of being sex-biased in either taxon.

A sex-biased gene expression index

Given that the sex-biased expression of genes is the basis for the generation of the sex-specific phenotypes, we were interested in assessing the individual variation in the cumulative expression of sex-biased genes for each organ. To do this with a single factor for every individual, we have developed a sex-bias index (SBI) based on the normalized expression values of the genes (transcripts per million = TPM). This is calculated as:

SBI = (MEDIAN of all female-biased tpm) - (MEDIAN of all male-biased tpm)

For the male-biased genes we exclude the ones encoded on the Y-chromosome, since they have no equivalent in females. The MEDIANs are used to reduce the influence of outlier values in some individuals.

When the sexes are well separated, one would expect two different distributions for each sex. The centers of the distribution can move along the X-axis, depending on the relative expression levels on males and females. If they are similar, the females should show a range of positive values representing their individual variation and the males should show an analogous range of negative values. But the variation ranges could also overlap, which would indicate that the phenotypes generated by the set of sex-biased genes are not completely sex-specific, similar as height or weight measures are not completely sex-specific.

The distributions of individual sex-bias indices are plotted as density functions across all mouse organs and taxa in Figure 4. We find indeed overlapping patterns in most comparisons of the non-gonadal organs. However, the organs differ with respect to their overlap patterns and these differences evolve fast between the taxa. This is particularly evident for kidney and liver. Kidney shows non-overlapping distributions in DOM and MUS, while in SPR and SPI, the male-biased genes show a very broad distribution that overlaps with the distributions of female-biased genes. Strong sexual dimorphism in mouse kidneys was previously described and is linked to corticoid signaling pathways, especially in males (Laulhe et al., 2021). In liver, DOM and MUS are also well separated, while SPR and SPI are overlapping. For mammary, the overlaps are also substantial in all four taxa, although one might have predicted that that they could be distinct. Such a distinctness is evident for all comparisons between gonadal organs, but these generate also very different morphological structures, i.e. an overlap is not expected.

Figure 4:

Given that the same sets of organs were retrieved from each individual, it is possible to ask whether SBIs might correlate between organs of the same individual. For example, an individual with a strong femaleness score in one organ might also have a strong femaleness score in another organ. However, this is not the case. There is no significant correlation in any pairwise comparison between organs for a given sex (Spearman rank correlation test, all p-values with multiple testing correction > 0.05) (suppl. Table S3). This implies that the femaleness and maleness status of an individual is not homogeneous throughout the body - each organ can be somewhat different in this respect.

Sex-biased gene expression in humans

The SBI can be applied to any comparative transcriptome data of population samples from both sexes. The GTEx consortium has generated such data for humans (Aguet et al., 2020) and these have previously been analyzed with respect to sex-specific expression patterns (Gershoni & Pietrokovski, 2017; Khodursky et al., 2022; Oliva et al., 2020). Unfortunately, these data are much more heterogeneous than the mouse data, with variable age distributions, data quality and death reasons. Various statistical procedures are therefore usually employed to control for confounding variables (Aguet et al., 2020; Khodursky et al., 2022; Oliva et al., 2020; Wolf et al., 2023). The largest problem is the overdispersion of the data, including the frequent occurrence of outlier values. The procedure that we have used for the mouse data addresses these problems stringently (see Methods) and we apply it therefore also to the human data. For comparison, we use also the lists of the sex-biased genes generated by the GTEx consortium that were generated with a more permissive procedure across all samples (Oliva et al., 2020). We focus here on individuals younger than 49 who did not die due to a long disease. These individuals are relatively rare in the GTEx data, given the strong bias for older individuals in the dataset. But we could retrieve nine females and males each for most organs and organ subsets (N =27; suppl. Table S4), making the patterns directly comparable to the mouse patterns. Our overall results are qualitatively very similar to the previously published results on these data, but given our more stringent filtering due to including an FDR step and an explicit cutoff, the overall numbers are lower, especially in comparison to (Oliva et al., 2020).

Compared to the mouse, we find generally fewer sex-biased genes in most human tissues (suppl. Table S5). Among 27 tissues included, only 10 have at least five sex-biased genes in either sex (Figure 5A). The tissues with the largest numbers of sex-biased genes are “Adipose Subcutaneous” and “Breast Mammary Tissue” (Figure 5B - note that the Y-axis scale of this figure is 15-fold expanded compared to Figure 5A). Further, the comparison between ovary and testis (OvaTes) between males and females show large numbers of sex-biased genes, but the tissues have evidently also very different cell compositions for these organs.

Figure 5:

Calculation of the SBI makes is limited when there are only very few sex-biased genes in one or both tissues. We are therefore plotting the SBI only for tissues with more than five sex-biased genes for either sex (Figure 5C). All of the SBI distributions are more or less overlapping between males and females, including also the breast tissue, similar as in mouse. Only the comparison between ovaries and testis (OvaTes) are very distinct, as to be expected.

When using the lists of sex-biased genes from Oliva et al. (Oliva et al., 2020) to calculate SBIs for the same individuals and same organs, we find qualitatively similar results but an even higher overlap between the sexes (Figure 5D). This is due to the inclusion of many additional genes that fall below our threshold of 1.25-fold change (compare numbers in suppl. Table S6).

Sex-biased expression in human single-cell data

Gene expression in organs is measured from complex aggregations of diverse cell types, making it difficult to distinguish between sex differences in expression that are due to regulatory rewiring within similar cell types and those that are simply a consequence of developmental differences in cell-type abundance (Darolti & Mank, 2023). It has actually been shown that only a subset of cells in a given tissue may express the sex-biased genes (Rodríguez-Montes et al., 2023). In complex organs with many cell types, such as the brain, it is therefore possible that the sex-bias signal is blurred by different sets of genes expressed in different cell types. Unfortunately, suitable datasets from single cell studies that include sufficient numbers of individuals where one could approach this question are still rare. We have analyzed here data from a study with patients that have developed Alzheimer’s disease. This study included also non-disease control individuals and we have used the datasets from these individuals, yielding comparisons between six individuals of each sex for dorsolateral prefrontal cortex (DLPFC) samples and seven individuals of each sex for middle temporal gyrus (MTG) samples (Gabitto et al., 2023). Sex-biased expression is still at a low level in these data (Table 1, suppl. Tables S7 and S8), similar as found for brain tissues in general and too low to generate meaningful SBI distribution plots. While this may be partly due to the age of the individuals, and/or insufficient coverage of low level expressed genes in single-cell datasets, it still implies that there is not any particular cell-type with very strong male-female differentiation.

Table 1:

Conserved genes with sex-biased expression

While most studies on sex-biased genes have reported that only a subset of them is conserved across larger evolutionary distances, they still report sometimes substantial numbers of genes with conserved sex-biased expression (Harrison et al., 2015; Khodursky, Svetec, Durkin, & Zhao, 2020; Naqvi et al., 2019; Rodríguez-Montes et al., 2023). However, given the high turnover of genes that we observe already within mouse lineages that are only at most 2 million years apart, it seems well possible that the seeming conservation among lineages with larger evolutionary distances is due to convergence, i.e. independent recruitment of genes into sex-biased expression in different lineages.

We have therefore sought to generate a set of genes with conserved sex-biased expression using stringent criteria. First, we retrieved all genes that show a consistent sex-biased expression in all four mouse taxa and in at least two of the analyzed tissues. Based on this set, we asked which orthologues of the genes show also sex-biased expression in human tissues. Interestingly, only very few genes are left when applying these filters (Table 2). There are eight orthologous genes that are female-biased in two or more tissues and three that are male-biased. Note, however, that the latter are all genes that are encoded on the Y-chromosome. The only two genes with consistent sex-biased expression across all organs and species are the female-biased gene Xist and the male-biased gene Ddx3y.

Table 2:

Four of the conserved genes (Kdm5c, Kdm5d and Kdm6a and Uty) are involved in histone methylation and demethylation, thereby triggering gene regulatory cascades that are known to be involved in setting up the sex-specific expression networks (Deegan, Karbalaei, Madzo, Kulathinal, & Engel, 2019; Samanta et al., 2022). Hence, the balancing of the interplay between these epigenetic regulators may well contribute to the fast evolution of sex-biased gene expression in their target genes. However, while all four are broadly expressed in adult tissues, they are themselves subject to differential sex-biased expression in a subset of the tissues (Table 2), implying that they depend also on sex-bias controlling expression networks.

Discussion

The sexual selection theory assumes that males and females have different interests in their reproductive strategies, which results in a continuous evolutionary conflict due to divergent trait optima expressed in either sex (Darwin, 1871; Hamilton, 1967; Mank, 2017a; Parsch & Ellegren, 2013; Price, Parkus, & Wright, 2023). The resulting sexual dimorphisms can occur at many levels, not only as prominent major phenotypic differences, but also as very many cryptic differences, as has e.g. been shown in a systematic survey of the phenotypes of mouse knock-out strains (van der Bijl & Mank, 2021). Many of these cryptic differences are expected to be caused by the sex-biased expression of genes (Grath & Parsch, 2016; Mank, 2017b). Hence, there is much room for the fast evolution of sex-biased gene expression due to the existence of continuous sexual antagonism. Our data show that sex-biased gene expression evolves extremely fast, as to be expected from a continuous conflict situation. Different sets of genes can get easily associated with the sex-biased expression networks and are also easily lost. Intriguingly, many of them reverse their role by switching between their sex-bias to the other sex. Similar results were found for the evolution of sex-biased genes in the Drosophila brain (Khodursky et al., 2020). Intriguingly, this fast turnover of roles is also reflected in the protein evolution patterns. Genes that are sex-biased in a given lineage show higher adaptive amino acid replacement rates than their non-sex-biased orthologues of a sister lineage.

Sex-biased genes do not only evolve fast between lineages, but tend also to be more variable within lineages. This implies that they contribute more to individual variation than many non-sex-biased genes. Interestingly, in humans, environmentally responsive genes show also higher variances than genes involved in regulating fundamental cellular processes (Wolf et al., 2023) and a subset of genes shows differential variability depending on which sex expresses them (Khodursky et al., 2022).

The SBI constitutes a cumulative index to compare the variances in sex-biases gene expression in individuals. Interestingly, it shows for many tissues overlapping distributions, especially also in humans. This implies that any individual that falls within the overlapping distributions would not be classifiable to belong to either sex. In addition, there is a general tissue-specificity of sex-biased gene expression (Naqvi et al., 2019; Oliva et al., 2020; Rodríguez-Montes et al., 2023) and we find of a lack of correlation of SBIs between organs of the same individual. This implies that individual males or females are actually composites with respect to displaying more female or more male characters throughout their body and brain, even though they may be easily classifiable based on their primary sexual characteristics. This supports the notion that the use of binary averages without recognizing individual variances in describing sex differences can be rather misleading (Ainsworth, 2015; Maney, 2016).

Conclusions

Our results show that most of the genetic underpinnings of sex-differences show no long-term evolutionary stability, which is in strong contrast to the perceived evolutionary stability of two sexes. This fast evolution is accompanied with high individual variability of sex-biased gene expression with overlapping distributions and different sex-biases in different organs. Individuals constitute therefore a continuous spectrum of sex characteristics that should not be cumulated into a simple binary classification. This is also relevant for the consideration of sex-specific medical treatments. While average differences in disease etiology between sexes are well documented (Mauvais-Jarvis et al., 2020), the decision on where a given individual falls into the spectrum of maleness/femaleness differences becomes more difficult, the more the variances overlap.

Acknowledgements

We thank Christine Pfeifle and Heike Harre for animal care taking and organ preparations, Michaela Schwarz for RNA purification and Wenyu Zhang and Luisa Pallares for discussion and suggestions. The study was funded by institutional resources of the MPG.

Author contributions

C.X. and D.T. designed the experiment, C.X. and S.K. run the experiment and collected the data, C.X. and D.T. did the data analysis and wrote the manuscript.

Methods

Mouse organ samples

Gene expression data were collected from outbred individuals from two Mus musculus subspecies, M. m. domesticus (DOM) and M. m. musculus (MUS) as well as two sister species, M. spretus (SPR) and M. spicilegus (SPI). Age-matched nine adult females and nine adult males were chosen from each of the four taxa, 72 individuals are included in total. As non-gonadal organs we included brain (whole brain), heart, liver (left medial lobe), kidney (right) and mammary gland (fourth, right). Note that the mammary glands in mice have similar sizes in both sexes before lactation and are therefore directly comparable. As gonadal organs we chose ovary (both ovaries), oviduct (both oviducts) and uterus (right uterine horn) from females and testis (right), epididymis (right) and vas deferens (right) from males. All organs were always retrieved from the same individuals, allowing expression comparisons between organs within each individual. In total, 576 samples are included (suppl. Table S1).

RNA sequencing and data analysis

The organs were carefully collected and immediately frozen in liquid nitrogen. Total RNAs were purified using QIAGEN RNeasy 96 Universal Tissue Kit (Catalog no. 74881), and prepared using Illumina TruSeq Stranded mRNA Kit, and sequenced using Illumina NovaSeq S4 (2 x 150 bp) in Kiel Sequencing Center. All procedures were performed in a standardized and parallel way to reduce experimental variance.

Raw sequencing reads were trimmed using Trimmomatic (0.38) (Bolger, Lohse, & Usadel, 2014). Only paired-end reads left were used for following analyses. The trimmed reads were mapped to mouse genome GRCm39 (Martin et al., 2023; Mouse Genome Sequencing et al., 2002) using HISAT2 (2.2.1) (Kim, Langmead, & Salzberg, 2015) with default parameter settings, except for “--score-min” set as “L,0.0,-0.6”, in order to compensate the sequence divergences of individuals from various taxa. Fragments mapped to the genes annotated by Ensembl (Version 104) were counted using featureCounts (2.0.3) (Liao, Smyth, & Shi, 2014).

Assignment of sex-biased gene expression

Various methods have been used to identify genes with sex-biased expression, also depending on whether the focus of the study was on gonadal differences (e.g. (Harrison et al., 2015) or on other organ comparisons (Blekhman, Marioni, Zumbo, Stephens, & Gilad, 2010; Naqvi et al., 2019; Oliva et al., 2020; Reinius et al., 2008; Rodríguez-Montes et al., 2023; Yang et al., 2006). Most of these methods employ some form of parametric statistics, such as a negative binomial distribution, which is necessary when only few samples are available per group. However, this is problematic in two ways. First, transcriptome data are well known for showing extreme outliers, and second, the variances of expression can be different among different genes with similar expression level, and between males and females. Thus, these parametric methods usually lead to exaggerated false positives (Li, Ge, Peng, Li, & Li, 2022). We used here at least nine individuals per sex which allows to perform non-parametric statistics, combined with setting cutoff criteria with a false discovery rate (FDR) correction. For each organ and each taxon, we normalized the fragment counts to TPM values, and added one to all TPM values (“TPM+1”) to avoid dealing with zeros in ratio calculations. Only genes with median of “TPM+1” in at least one sex larger than 2 were kept for analysis. The sex-bias ratio for each gene is calculated as “MEDIAN of females / MEDIAN of males”. We have explored a range of cutoffs and found that a sex-bias ratio of 1.25-fold difference of MEDIAN expression values combined with a Wilcoxon rank sum test and Benjamini-Hochberg FDR correction (using FDR <0.1 as cutoff) (Benjamini & Hochberg, 1995) yields the best compromise between sensitivity and specificity. For controls with non-biased genes, we chose a cutoff of a smaller than 1.05-fold ratio. All steps of the full procedure are listed in suppl. Table S9.

For the non-gonadal organs, we used the comparison between nine females and males each to identify genes with sex-biased expression. For the gonadal organs, we chose ovary versus testis (OvaTes), oviduct versus epididymis (OviEpi), and uterus versus vas deferens (UteVas) for the pairwise comparisons. Except for ovary versus testis, the other two pairs are not really considered to be homologous, but they serve very roughly comparable general functions.

McDonald–Kreitman test

The asymptotic McDonald–Kreitman test (Messer & Petrov, 2013) was performed based on the single nucleotide variants (SNVs) of wild mouse samples derived from the whole genome sequencing data of eight mice of each of the three taxa: M. m. domesticus (DOM), M. m. musculus (MUS), and M. spretus (SPR) (downloaded from “https://wwwuser.gwdg.de/~evolbio/evolgen/wildmouse/vcf”) (Harr et al., 2016). The SNVs were annotated as synonymous and nonsynonymous using snpEff (4.3t) (Cingolani et al., 2012). For each SNV in each taxon, its type, polymorphic or fixed, was determined, and the allele frequency was calculated if it is polymorphic. A fixed SNV in M. spretus was considered to be the ancestral allele, and then the derived allele can be determined in M. m. domesticus and M. m. musculus. Only SNVs on the autosomes were included for the analysis. The numbers of fixed synonymous and nonsynonymous sites, and the numbers of polymorphic synonymous and nonsynonymous sites at each derived allele frequency were used as input for the adapted R script “asymptoticMK_local.R” downloaded from “https://github.com/MesserLab/asymptoticMK”.

The asymptotic McDonald–Kreitman test was performed on eight sets of genes: DOM female-biased genes, DOM male-biased genes, DOM non-biased genes, DOM non-biased genes of which the orthologs are sex-biased in MUS as control, MUS female-biased genes, MUS non-biased genes, MUS male-biased genes, and MUS non-biased genes of which the orthologs are sex-biased in DOM as control.

Data from humans

The human bulk RNA-Seq data were retrieved from the GTEx project (Consortium, 2013, 2020) as TPM files for each organ. To make them best comparable with the mouse data, we chose nine females and nine males for each organ, using the following criteria: maximum age 49 years, death reason category not 0. The full sample list is provided in suppl. Table S4. Note that in contrast to the mouse data, the different organ data come mostly from different individuals.

The human single-nucleus RNA-Seq (snRNA-Seq) data were retrieved from the SEA-AD project (Gabitto et al., 2023). “h5ad” format files, including the normalized read counts using a log transformation of pseudocounts per 10,000 reads, ln(CPTT+1), were downloaded from CZ CELLxGENE. Cells in two brain regions were included: middle temporal gyrus (MTG) and dorsolateral prefrontal cortex (DLPFC). Considering that it is a study of Alzheimer’s disease, we specifically chose “control” samples from individuals without dementia, and with race as “white” for analysis. We only used the data sequenced by the “10x 3’ v3” assay. These criteria lead to seven females and seven males for middle temporal gyrus and six of them each for dorsolateral prefrontal cortex. Only cell types (”subclass”) having at least 100 cells for each of all individuals were kept. The full list of cell types and individuals included is provided in suppl. Table S7 and S8. The expression level of a gene in a cell type of a brain region from an individual was calculated as the mean of the pseudocounts per 10,000 reads of all the cells belonging to the cell type, mean(CPTT+1).

Data availability

The ENA BioProject accession numbers for the wild mouse sequencing data reported in this study are PRJEB36991 and PRJEB50011. The tables with TPM values for all genes and organs, as well as the list of all mouse sex-biased genes are available for download at https://datadryad.org/stash/share/NsE5wg7fQ_b5czVWOlBnREhAfJGSzdT_juSjfghYLlw

Animal ethics and permissions

All mice were kept according to FELASA (Federation of European Laboratory Animal Science Association) guidelines, with the permit from the Veterinäramt Kreis Plön: PLÖ-004697. The Government of Schleswig-Holstein provided permission to sacrifice animals under permit number V312-72241.123-34.