Introduction

The usual narrative for sex-related differences between individuals is based on binary gametic sex (Goymann, Brumm, & Kappeler, 2023). Such a 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 can evolve fast, when they are subject to sexual selection (Darwin, 1871). The different interests of males and females lead to sexual conflict, characterized by opposing evolutionary constraints on the genes that mediate sex differences. Sexual conflict can have many forms but can be broadly categorized into interlocus and intralocus sexual conflict, with distinct implications for evolutionary biology and species fitness (Schenkel, Pen, Beukeboom, & Billeter, 2018). Each is also expected to have an effect on the evolution of the underlying genes and their regulation (Ingleby, Flis, & Morrow, 2015; Mank, 2017a; Tosto, Beasley, Wong, Mank, & Flanagan, 2023).

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 et al., 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; van der Bijl & Mank, 2021). 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).

The studies of morphological sex-differentiation of the human brain have revealed mosaic patterns where different brain parts of the same individual can show more female or more male characteristics (D Joel et al., 2015). Hence, a binary narrative which is based on gametic sex is not applicable at this level and could be replaced by relative degrees of maleness versus femaleness (D. Joel, Garcia-Falgueras, & Swaab, 2020). 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 degrees of somatic sex-specific variation between individuals even larger.

Sex-biased gene expression is known to evolve fast, supporting the notions of sexual conflict 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 turnover 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 two million years (Figure 1). Given the short reproduction cycles and the fast molecular evolution in rodents, this time corresponds in molecular divergence terms approximately to the equivalent of the evolutionary split between gibbons to humans. Hence, the system allows to study the continuum from microevolutionary to macroevolutionary patterns. 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).

Figure 1:

Most previous studies on sex-biased genes have given special emphasis to the evolution of gene expression in gonadal tissues even when somatic tissues were co-sampled, since the gonads harbor the largest numbers of sex-biased genes (Dean et al., 2017; Harrison et al., 2015; Pointer, Harrison, Wright, & Mank, 2013; Todd et al., 2018). In contrast, our study is mostly focused on somatic tissues that occur in both sexes, but we complement it with comparisons to gonadal tissues that are specific for each sex. 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. We use sets of 18 individuals from each taxon which 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 in the somatic organs, but with marked differences in numbers between the organs in the different evolutionary lineages. As it was previously found, most sex-biased genes are organ-specific (Naqvi et al., 2019; Rodríguez-Montes et al., 2023) and we show here that they are expressed in organ-specific modules. Only a small percentage of genes shows conserved 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 somatic organs, we find variable patterns of distinction and overlap between the sexes. 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. The SBI can also be used to assess mosaicism of sex-biased gene expression between individuals analogous to the brain (D Joel et al., 2015). We show that individuals are characterized by different degrees of mosaicism in different organs.

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 RNA-Seq analysis. We used five somatic organs and three gonadal organ parts for each of the sexes (Figure 1), yielding a total of 576 samples (Supplementary file 1 - Table S1).

To identify sex-biased expression, based on normalized transcript fragment count data (transcripts per million – TPM), we use the ratio of the medians of female (F) and male (M) expression (F/M) with a 1.25-fold ratio cutoff, combined with a Wilcoxon rank sum test and FDR correction < 0.1. These criteria are slightly more stringent than other studies that have focused on somatic sex-expression differences (see Methods for a discussion of these parameters). The fold-differences between male and female expression are continuously distributed, at least for the organs with many sex-biased genes (see distribution histograms in Figure 1 - figure supplement 1).

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 Supplementary file 1 - Table S2, the corresponding gene lists are provided in supplementary data D1. The numbers range between 13 - 5,645 for the somatic organs and 3,444-13,383 for the gonadal pairwise comparisons. Note that the gonadal comparisons are inherently very different to the somatic comparisons, since they are based on comparing morphologically different structures (see Methods for the scheme of pairwise comparisons). Still, we include them as reference comparisons here, but the main focus of our study is on the somatic comparisons.

Variable numbers of sex-biased genes are shared between pairwise organ comparisons within the different taxa, but only very few sex-biased genes are sex-biased in all organs in a given taxon (Supplementary file 1 - Table S3). These include female biased genes such as the X-chromosomally encoded gene Xist, which codes for a lncRNA required for the dosage compensation of the X-chromosome in females (Loda & Heard, 2019) and Eif2s3x, which encodes the eukaryotic translation initiation factor 2 subunit 3 (Xu, Watkins, & Arnold, 2006).

Most other genes that share a sex-bias between all tissues in each taxon are also encoded on the X-chromosome and most are known to escape the dosage compensation of the X-chromosome in tissue specific patterns (Berletch et al., 2015) (Supplementary file 1 - Table S3). Given that we find different sets of these genes in the different mouse taxa (Supplementary file 1 - Table S3), we can conclude that these genes show also some degree of evolutionary turnover in their regulation.

It has previously been reported that the X-chromosome harbors more female-biased than male-biased genes for the mouse somatic organs (B Reinius et al., 2012). Our data confirm this general observation (numbers and percentages included in Supplementary file 1 - Table S2). When counted across all five somatic tissues, the X-chromosome has a higher percentage of female-biased genes (4.7% female-biased versus 1.8% male biased), but given that the overall fraction of genes on the X-chromosome is also 4.7%, this reflects a deficiency of male-biased genes on the X. The sex-biased genes on the gonadal organs show the opposite trend (3.7% female-biased versus 4.3% male biased). However, there are substantial differences between the different organs and taxa in this respect, implying that there is no single explanation for the evolution of sex-biased genes on the X-chromosome.

Female-biased genes outnumber male-biased genes in all comparisons. Among the somatic 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 somatic 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 SPR and SPI and liver has the highest number of female-biased genes in SPI. For the gonadal organs, the overall numbers remain more constant between the taxa.

Evolutionary turnover of sex-biased gene expression

The evolutionary turnover of genes that become subject to sex-biased expression in any given taxon is very high, even between the closely related taxa that are studied here. Figure 2A 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 somatic 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.

Figure 2:

The fractions of shared genes are higher for the gonadal tissues (Figure 2A), 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 (Supplementary file 1 - Table S2). Hence, although the overall numbers of sex-biased genes do not change much between the taxa for the gonadal organs, the actual composition of the gene sets evolves substantially.

To compare the rates of this turnover of sex-biased gene regulation with non-sex-biased regulatory turnover, we use the same statistical approach for the latter group of genes, but instead of comparing males with females, we use a comparison between the same sex data groups between the taxa. For example, for a turnover comparison between DOM and MUS, we swapped the DOM males with the MUS females in the same statistical algorithm that we use for the sex-bias analysis. We did this also for the SPR-SPI and the MUS-SPR pairs, as well as the reverse sex groups, i.e., a dataset with comparison among females only (“Female Swaps”) and a dataset with comparisons among males only (“Male Swaps”). From the resulting genes lists, we re-sampled 1,000 times 1,000 genes each and calculated averages and standard deviations. The corresponding data and statistics tables are provided as supplementary data D2. The percentages of sex-biased genes that show a regulatory change were then subtracted from the average percentages of the equivalent all-genes comparisons in the same taxa groupings. Hence, when sex-biased gene expression turnover is higher than the all-genes expression turnover, one should find negative values. The results are shown in Figure 2B. Except for the kidney in the DOM-MUS comparison, all somatic organs in all comparisons show negative values, while the gonadal organs in all comparisons mostly show positive values.

This analysis shows that sex-biased gene expression changes occur much more frequently than for the whole gene set. Only the kidney in the DOM-MUS comparison are very different in this comparison, possibly since they have been subject to other rapid evolutionary processes that are currently unknown. The positive value for kidney in the DOM-MUS comparison is mostly due to a relatively high numbers of sex-biased genes that are conserved between these two taxa, for as yet unknown reasons. This pattern changes for the two other taxa comparisons.

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, both at the gene expression level, as well as the protein divergence level.

As a measure of variance for the gene expression, we used the interquartile range IQR (difference between the 75th and 25th percentiles of the data) divided by the median of the data range (representing the non-parametric version of the coefficient of variation). This measure buffers against outlier values that can frequently be observed in the data (Xie et al., 2024). For these comparisons we used also a set of non-sex-biased genes, defined as genes having a F/M ratio of < 1.05-fold. Sex-biased genes show in most cases higher variances than non-biased genes, but to different extents in the two sexes and in the different organs (Figure 3A). The pattern is much more distinct for somatic genes than for gonadal genes. To assess whether sex-biased genes are recruited from genes with higher general expression variances, we created reciprocal gene sets with the orthologous sex-biased and non-sex-biased genes in pairwise taxa comparisons (DOM - MUS and SPR -SPI). These show variance levels more comparable to those found among the sex-biased genes (Figure 3A), implying that their generally higher variation might help them to become more easily sex-biased, as it has also been suggested for plant sex-biased genes (Scharmann, Rebelo, & Pannell, 2021).

Figure 3:

Protein divergence patterns have also been studied for sex-biased genes. 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). Sex-biased genes in sex-dimorphic plant leaves, on the other hand, did not show increased divergence compared to non-sex-biased genes (Scharmann et al., 2021). However, all the previous tests on this question were a somewhat limited by not taking possibly confounding demographic factors into account.

A frequently used test for elevated positive selection in population comparisons, is the McDonald-Kreitman (MK) test. It compares the level of polymorphism versus fixation of coding and non-coding 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 somatic 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 to correctly assign coding and non-coding sites. The polymorphism data were taken from previous whole-genome sequencing studies of these populations (Harr et al., 2016). Classification of fixed versus polymorphic variants was done by comparing to SPR as outgroup. Statistical estimates are based on 1,000 bootstrap replications.

The results are shown in Figure 3B. Both, female-biased and male-biased genes show higher alpha values compared to non-biased genes. Interestingly, female biased rates are higher than male-biased rates, corroborating the observation in brain expressed genes in birds (Mank et al., 2007). Similar as for the expression variance comparisons, we generated also reciprocal gene sets, i.e. sex-biased in either DOM or MUS, but not sex-biased in the respective other taxon. The alpha values for the non-sex-biased reciprocal orthologs are significantly lower (Figure 3B) than those for the sex-biased genes. This comparison suggests that the evolution to a sex-biased expression state is often accompanied with an enhanced adaptive substitution rate at coding positions of the respective genes.

Module analysis of sex-biased genes

We were interested to assess whether sex-biased genes represent a subset of all genes expressed in a given organ, or whether they are derived from particular co-regulated modules. We used weighted gene co-expression network analysis (WGCNA) (Langfelder & Horvath, 2008) to determine gene expression modules for each organ, based on a set of transcriptomes from 48 DOM females that were prepared in parallel to the remainder of the samples, because WGCNA requires such large sample sets to generate reasonably robust results. The analysis of the scale independence and mean connectivity plots for the different organs in this data set show that saturation is reached at levels of soft power β between 4 to 8 (see Methods). Since a similarly large set is currently not available for males, we restrict the further analysis to the female-biased genes, which constitute actually the larger fraction of sex-biased genes between the sexes.

The module assignments for each of the five somatic organs are listed in Supplementary file 1 - Table S7. We found between 21 to 64 modules in the different organs. The numbers of genes in each module for each taxon and tissue are provided in Supplementary file 1 - Table S8 (note that module numbers are not comparable between the organs - they are always only determined for a given organ). The comparisons show that the sex-biased genes are not simply a proportional reflection of all modules, but are enriched for some subsets of modules. To directly visualize this, we subtracted the normalized fraction of gene numbers in each module from the normalized fraction of sex-biased genes for each taxon. The full data are provided in Supplementary file 1 - Table S8, a subset of the data is plotted in Figure 4. This subset is restricted to four somatic organs, since brain has too few sex-biased genes overall. Further, we did not include higher-numbered modules with few genes only, unless the differences to the fractions of all genes is more than 3% in the sex-biased genes.

Figure 4:

The results suggest that the module classes for the sex-biased genes tend to be retained between the taxa, despite the turnover of genes. This is further detailed for each organ below.

For heart, there is an excess of genes in module 1 for DOM, MUS and SPR. Module 1 includes 3010 genes associated with multiple physiological processes in a GO analysis. DOM includes 1247 of these among the total of 2405 sex-biased genes (= 52%). MUS includes 171 out of 337 (=51 %) and SPR 30 out of 76 (= 39%). In MUS, 54 out of the 171 sex-biased genes in module 1 are not sex-biased in DOM, and in SPR there are 6 sex-biased genes in module 1, that are neither present in DOM nor in MUS. Hence, there is a substantial turnover of sex-biased genes within module 1 between the taxa. In SPI, on the other hand, there is mostly a relative loss of sex-biased genes in module 1 (it retains only 6 such genes), but a relative major gain in a new module. This is module 26 with a total of 20 genes, 7 of which occur in SPI as sex-biased genes (four are transcription factors involved in circadian clock entrainment). The profiles for the taxon-specific sex-biased genes (blue) are very similar to all sex-biased genes (red) in each comparison, implying that both are parts of the same modules in each taxon.

For kidney, the patterns are more heterogeneous, with little commonality between the taxa. Also, the patterns between all sex-biased genes and that taxon specific sex-biased genes differ for some of the modules.

For liver, it is mostly module 2 that is enriched for sex-biased genes. The module has a total of 389 genes, which are enriched for innate immunity regulation GO terms. Here, MUS and SPI have with 47% and 52% the largest fractions of sex-biased genes in this module. The patterns of the recruitment of taxon-specific sex-biased genes is somewhat heterogeneous and module 2 is not specifically notable.

For mammary, it is module 6 that is most strongly enriched for sex-biased genes. The module has 540 genes, which relate mostly to epithelial cell morphology and tubule formation. MUS and DOM have the largest fractions of sex-biased genes in this class, with 83% and 73% respectively. Again, the patterns for the taxon-specific sex-biased genes are somewhat heterogeneous, module 6 is not specifically notable. But the changed overall pattern in SPI is evidently driven by genes newly recruited to become sex-biased in a subset of modules in this taxon.

We conclude from this analysis that sex-biased genes tend to be expressed in a subset of modules and that evolutionary turnover may sometimes occur by recruiting previously not sex-biased genes from the same module to become sex-biased. However, this pattern is clearly heterogeneous between organs and taxa. It applies to some comparisons, but not to all.

A sex-biased gene expression index

Given that the sex-biased expression of genes is thought to be 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. 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 (Xie et al., 2024).

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. While this provides a measure of the overall expression levels of sex-biased genes in the organs, this may in itself be of less relevance. For example, in plants it was shown that the sex-biased morphological variation of the leaves does not correlate with the numbers and expression of sex-biased genes in these leaves (Scharmann et al., 2021). But independent of considering the numbers and absolute expression levels of the genes, one can look at the variance of expression between the individuals for each organ and how this variance evolves between the taxa. Such a comparison has not been done before in any dataset and the SBI is uniquely suitable to visualize such comparisons (see also Discussion).

To focus on the variance aspect of the SBI, we normalize the values to center the distributions around 0. In this display, 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 most somatic phenotypes, such as height measures or brain structures are not completely sex-specific (D Joel, 2021; Maney, 2016).

The normalized distributions of individual sex-bias indices are plotted as density functions across all mouse organs and taxa in Figure 5. We find indeed overlapping patterns in most comparisons of the somatic 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 and livers was previously described and is linked to hormonal signaling pathways, especially in males (Clodfelter et al., 2006; 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 more 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 5:

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) (Supplementary file 1 - Table S9). 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.

Mosaic patterns of sex-biased gene expression

The comparative analysis of sex-biased brain structures based on MRI datasets of more than 1,400 human brains has suggested that brains are usually composed of a mosaic of sex-biased structures (D Joel et al., 2015; D. Joel et al., 2020). The above correlation analysis of SBI values between mouse organs has also suggested that such a mosaic pattern could apply (see above). To visualize this more directly, we use the heatmap approach developed by (D Joel et al., 2015) to plot the SBI values for each individual and organ via a normalized color scale (Figure 6). These plots show that the SBIs differentiate the individuals such that each has a more or less unique combination. While the patterns of males and females are distinct, they show also overlaps for individual values in the middle ranges.

Figure 6:

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 after a long disease phase. 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 (27 female-male comparison sets; Supplementary file 1 - Table S10), 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 (Supplementary file 1 - Table S11). Among 27 tissues included, only 10 have at least five sex-biased genes in either sex (Figure 7A). The tissues with the largest numbers of sex-biased genes are “Adipose Subcutaneous” and “Breast Mammary Tissue” (Figure 7B - note that the Y-axis scale of this figure is 15-fold expanded compared to Figure 7A). 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 7:

Calculation of the SBI 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 7C). 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 ovary 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 7D). This is due to the inclusion of many additional genes that fall below our threshold of 1.25-fold change (compare numbers in Supplementary file 1 - Table S11).

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 single cells from dorsolateral prefrontal cortex (DLPFC) samples and seven individuals of each sex for single cells from middle temporal gyrus (MTG) samples (Gabitto et al., 2023). There are 15 different cell types with sufficient coverage for an analysis in the DLPFC data and 13 in MTG data (Supplementary file 1 - Tables S13 and S14). A subset of 12 cell types overlaps and can therefore be compared between the tissues with the results shown in Table 1 in and detailed in Supplementary file 1 - Table S15. Sex-biased expression is still at a low level in these data, similar as found for brain tissues in general and too low to generate meaningful SBI distribution plots.

Table 1:

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.

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 two 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. This has also been suggested in a phylogenetic study on sex-biased genes in plants (Scharmann et al., 2021).

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 (Supplementary file 1 - Table S16). We then focused on those genes that occur in at least two of the analyzed tissues. This filters for genes that can be expected to have a more general function in regulating sex-biased expression. Based on this set, we asked which orthologues of the genes show also sex-biased expression in any of the human tissues. Interestingly, only very few genes are left when applying these filters (Figure 8). 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.

Figure 8:

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 (Figure 8), implying that they depend also on sex-bias controlling expression networks.

Discussion

The present study is the first comparative study on sex-biased gene expression in animals that includes sufficient data to allow within-population comparisons in parallel to a phylogenetic approach from closely related taxa. Further, its focus is on somatic sex-biased expression of genes, which is of special relevance for understanding the variances of sex-related phenotypes of individuals in populations. The data show that sex- biased gene expression evolves fast, even between the closely related taxa studied here. Many genes can 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) and in sex-biased genes expressed in sexually dimorphic leaves of plants (Scharmann et al., 2021). Variances of sex-biased gene expression between individuals are often overlapping between the sexes, which can be shown by plotting sex- bias index (SBI) values.

Neutral or adaptive?

Patterns of fast evolution of molecular characters such as gene expression or substitution rates raise always the question of which fraction could be due to neutral divergence versus adaptive processes, i.e., positive selection for new variants. This question has previously been studied in the same system of mouse populations and taxa that are used here. These studies showed that the overall patterns of gene expression changes between these populations are mostly compatible with neutral divergence models, based on tests of intra-group variability with between-group divergence (Bryk, Somel, Lorenc, & Teschke, 2013; Fabian Staubach, Teschke, Voolstra, Wolf, & Tautz, 2010; Voolstra, Tautz, Farbrother, Eichinger, & Harr, 2007).

Given that sex-biased genes appear to be recruited from a set of genes that show generally a higher variance (compare Figure 3A), it is thus possible that the sex-bias status is due to neutral evolution mechanisms, in which segregating variants in promotors that respond to sex-specific regulatory mechanisms become randomly fixed over time. However, this assumption is not compatible with the finding that the turnover of sex-biased gene expression is faster than non-sex-biased gene expression (Figure 2B).

For a neutral mechanism acting at the promotor levels, one would expect that different organs and different taxa show similar overall numbers of sex-biased genes and that these should be randomly recruited from the different gene expression modules. Both of these ptterns are clearly not the case (see Figures 1 and 4). For example, with a neutral mechanism it would be difficult to explain why DOM has such a large number of sex- biased genes in the heart, most of which belong to only one module of gene expression. Similarly, also the other module-based patterns in other organs could not be explained by such a mechanism.

The alternative for explaining the module patterns would be that transcription factors that regulate these modules become sex-biased through neutral mechanisms, with the consequence that they drag with them a number of genes from within the respective modules to make them sex-biased. However, this would not explain why we see a faster adaptive protein evolution among the sex-biased genes. If their sex-bias status is passively driven by neutrally evolving transcription factors, one would not expect a differential selection pressure acting on their protein sequences.

In the combination of these considerations, we favor an interpretation where positive selection drives sex- biased gene expression, probably including also transacting regulatory genes to explain the module patterns. It is generally known from selective sweep studies (Ihle, Ravaoarimanana, Thomas, & Tautz, 2006; F. Staubach et al., 2012; Teschke, Mukabayire, Wiehe, & Tautz, 2008) and amino acid substitution comparisons (Halligan, Oliver, Eyre-Walker, Harr, & Keightley, 2010) that the mouse populations are subject to massive positive selection. The high propensity for positive selection in mice is likely due to their short generation times and high effective populations sizes. Both are very different in humans, which might explain the generally smaller numbers of sex-biased genes in humans.

Role of sexual selection

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 adaptive evolution of sex-biased gene expression due to the existence of continuous sexual conflict (Cox & Calsbeek, 2009) or sex-specific selection, e.g. (Oliver & Monteiro, 2011).

Many previous studies on sex-specific expression have focused in differences between gonads. But gonads are composed of very different cell types and tissues in males versus females. Hence many of the sex-biased genes in gonads may simply be sex-biased because of the different cell-types that generate these different tissues.

Yet, gonads produce also the hormones that are required for the sex-specific differentiation of the tissues, including the gonads themselves. Hence, in populations where individuals assume intermediate reproductive tactics between males and females, it was also possible to correlate this intermediacy with corresponding patterns of sex-biased gene expression in the gonads (Dean et al., 2017; Pointer et al., 2013; Todd et al., 2018). In a study in fish that show individuals with female mimicry, it was possible to find also differences in brain expression patterns and modules according to reproductive tactics (Cardoso, Gonçalves, Goesmann, Canário, & Oliveira, 2018).

These studies support the notion that sexual selection can be a strong driver of adaptive sex-biased gene expression patterns. This could include patterns that are evolving in only one of the sexes, or antagonistic sexual selection to resolve evolutionary conflicts between the sexes. Given the still relative lack of studies that focus on somatic sex-biased genes, it is currently not possible to say which of these mechanisms is more prevalent for somatic tissues.

Intragroup variation and the SBI

Sex-biased genes do not only evolve fast between lineages, but they represent also the fraction of genes that are more variable within lineages. This implies that they contribute more to individual variation than 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-biased gene expression in individuals. It has some intrinsic features that make it particularly suitable to reflect the individual variances. In contrast to measures that would focus on either male-biased or female-biased genes, it puts both into relation to each other, by subtracting the respective median expressions of male sex-biased genes from those of female sex- biased genes. The use of the medians, rather than the averages reduces the impact of genes with extremely high sex-bias because of specific molecular functions, e.g., the gene Xist that is required in every cell for the dosage compensation of the X-chromosomes in females. However, note that even this gene shows a rather variable expression between females (values included in supplementary Data D1)

The use of the cumulative median across all genes in either sex-biased class is chosen following the general concepts of “genomic prediction” (Meuwissen, Hayes, & Goddard, 2001) which are based on summarizing values derived from the polygenic view of the generation of the phenotype (Tautz, Reeves, & Pallares, 2023). Such a cumulative view across the effects of all expressed genes in a tissue is nowadays also developed for understanding human genetic diseases, summarized as a polygenic risk score (Crouch & Bodmer, 2020; Gibson, 2019). The SBI represents a similar approach, but specifically for studying sex-specific differences.

The use of a subtraction rather than a ratio between each of the sex-biased medians has the advantage of implicitly retaining the information on the relative expression levels, as well as allowing a simple comparison on a linear scale. Linear scales are also usually used to compare distributions of morphological measures especially in the context of visualizing overlapping sex-specific patterns (Maney, 2016). If one wants to focus on the comparison of variance patterns only, one can easily scale the SBI values to center them around zero and use a smoothing function, as we have done in Figures 5 and 7. Normalized SBI values can also be used in a heatmap visualization to represent a mosaic pattern of sex-biased differences (D Joel et al., 2015), as we have done it in Figure 6.

The SBI distribution comparisons in Figures 5 and 7 show for many tissues overlapping distributions, especially in humans. This implies that any individual that falls within the overlapping distributions would not be classifiable to belong to either sex. The gonadal distributions are always very distinct and non-overlapping, but this can be explained by the different tissue and cell-types, as discussed above. Interestingly however, we find also distinct distribution of SBI variances for kidney and liver in DOM and MUS mouse populations. This shows that even somatic organs that do not differ in their cell types can show very distinct SBI variances. Interestingly, the two other taxa, SPR and SPI show overlapping distributions for these organs, implying that the very distinct patterns in DOM and MUS are based on adaptive differences between the sexes that are maintained by some form of sexual selection in these species. Note that it is known that the gene expression differences in these organs are directly regulated by hormone action (Clodfelter et al., 2006; Laulhe et al., 2021), but this does not explain why the strong differences exist in the first place, especially when a similarly strong differentiation is not seen in the other mouse taxa or humans.

There is a general tissue-specificity of sex-biased gene expression in mammals (Naqvi et al., 2019; Oliva et al., 2020; Rodríguez-Montes et al., 2023; Yang et al., 2006) 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. Hence, while hormones are undoubtedly the drivers of dimorphic sex-differentiation, they act via different sets of transcription factors in the different organs (Williams & Carroll, 2009). This allows to generate independent variances with respect to sex-biased gene expression in different cell types and different organs, which are reflected in the SBI distributions.

The generally broader overlaps in SBI distribution in humans compared to mice can be ascribed to two factors. First, there are fewer sex-biased genes in humans, possibly because the relatively lower population size in humans cannot maintain so many adaptive differences as in mice (see above). In addition, the human samples are much more heterogeneous with respect to age and environmental exposure, compared to the mice that were collected at very similar ages and grew up under controlled environmental conditions. However, this implies that humans reflect the more natural conditions, while the mouse data represent a somewhat more artificial accentuation of the genetic differences. Evidently, both, for the action of adaptive processes, as well as for our own thinking about sex-differences, the natural conditions are more relevant - and they seem to favor larger overlapping variances.

Conclusions

Our results show that most of the genetic underpinnings of sex-differences show no long-term evolutionary stability. 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 (Ainsworth, 2015; Maney, 2016; Sharpe et al., 2023). 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. The fast evolution of sex-biased expression is also a warning sign that mouse models may not be suitable for developing gender specific medicinal treatments for humans, especially in view of the fact that humans have by far fewer sex- biased genes than mice.

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). Nine age-matched adult females and adult males each were chosen from each of the four taxa, 72 individuals are included in total in the overall analysis. As somatic 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 were included (Supplementary file 1 - Table S1).

RNA sequencing and data analysis

The organs were carefully prepared and immediately frozen in liquid nitrogen. Total RNA was 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 or on other organ comparisons (Blekhman, Marioni, Zumbo, Stephens, & Gilad, 2010; Harrison et al., 2015; Naqvi et al., 2019; Oliva et al., 2020; B. 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 > 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. This was assessed by comparing sex-randomized datasets with the actual data. In these pre-tests we found that the inclusion of the Wilcoxon test with FDR correction was most effective in increasing the contrast between randomized scores and actual scores to an at least 20-fold difference. For controls with non-biased genes, we chose a cutoff of a < 1.05-fold ratio without Wilcoxon test. All steps of the full procedure are listed in Supplementary file 1 - Table S17.

For the somatic 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.

To obtain the distribution of alpha in each McDonald–Kreitman test, we used 1000 rounds of resampling via bootstrapping from variants. The results are included in Supplementary file 1 - Table S6.

WGCNA analysis

In order to assess co-regulatory modules among the sex-biased genes, we used the weighted gene co- expression network analysis (WGCNA) approach (Langfelder & Horvath, 2008) to generate module assignments. According to its manual, our sample size, nine, is too low to generate biologically meaningful gene networks. We made therefore use of data from a parallel study (Xie et al., 2024) for which we obtained 39 additional DOM female samples for each of the somatic organs in the present analysis (brain, heart, kidney, liver, and mammary gland). Note that all the samples were obtained in parallel at the same time by the same people and with the same procedure. The total sample size, 48, is sufficient for a reliable WGCNA analysis.

Following the standard procedure, we first removed the outlier samples: none for brain and mammary, and one for each of heart (H35232), kidney (H35331), and liver (H35302). Then we chose the parameter “power” with the help of its “pickSoftThreshold” function as 4 for brain, 8 for heart, kidney, and mammary, and 6 for liver. Finally, we ran its main function “blockwiseModules” to get module assignments of genes in each organ with parameters as “corType = ’bicor’, maxPOutliers = 0.1, networkType = ’signed hybrid’, TOMType = ’signed’”. The module assignments for all genes are included in the supplemental data D1.

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 Supplementary file 1 - Table S10. Note that in contrast to the mouse data, the different organ data come from different sets of 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 (Abdulla et al., 2023; Megill et al., 2021). 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 Supplementary file 1 - 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 supplemental data folders D1 - D3 including the tables with TPM values for all genes and organs are available for download at the Edmond repository of the MPG: https://edmond.mpg.de/privateurl.xhtml?token=d97a1f75-75e0-4177-a048-b186d23f74ae

This will be converted into a published dataset after final publication of the paper.

Code availability

The code for sex-biased gene detection and analysis is at “https://github.com/cxiepku/sex-biased”

Acknowledgements

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

Additional information

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.

Animal ethics and permissions

Maintenance and handling of the animals were conducted in accordance with German Animal Welfare Act and FELASA guidelines. The project was approved with the number 1158 by the Animal Welfare Officers of the University of Kiel according to the German Animal Welfare Act §4 “Killing animals and organ withdrawal for scientific purpose”. Permits for keeping mice were obtained from the veterinary office “Veterinäramt Kreis Plön” under permit number: PLÖ-000 4697 (08.04.2014).

Additional files

Supplemental File 1 with all supplementary tables

Figure supplement 1. Histograms of distributions of female/male sex-bias ratios for all mouse and human organ datasets.

Figure 4 - Figure supplement 1. Scale independence and mean connectivity plots for the different organs for determining the soft threshold parameter β for the WGCNA module assignments.