Rapid Evolution of Complete Dosage Compensation in Poecilia.

Dosage compensation balances gene expression between the sexes in systems with diverged heterogametic sex chromosomes. Theory predicts that dosage compensation should rapidly evolve in tandem with the divergence of sex chromosomes to prevent the deleterious effects of dosage imbalances that occur as a result of sex chromosome divergence. Examples of complete dosage compensation, where gene expression of the entire sex chromosome is compensated, are rare, and have only been found in relatively ancient sex chromosome systems. Consequently, very little is known about the evolutionary dynamics of complete dosage compensation systems. Within the family Poeciliidae the subgenus Lebistes share the same sex chromosome system which originated 18.48-26.08 Ma. In Poecilia reticulata and P. wingei, the Y chromosome has been largely maintained, whereas the Y in the closely related species P. picta and P. parae has rapidly degraded. We recently found P. picta to be the first example of complete dosage compensation in a fish. Here, we show that P. parae also has complete dosage compensation, thus complete dosage compensation likely evolved in the short (∼3.7 Myr) interval after the split of the ancestor of these two species from P. reticulata, but before they diverged from each other. These data suggest that novel dosage compensation mechanisms can evolve rapidly, thus supporting the longstanding theoretical prediction that such mechanisms arise in tandem with rapidly diverging sex chromosomes.

Introduction

In organisms with heterogametic sex determination, the Y chromosome diverges from the X when recombination between them is suppressed (Furman et al. 2020). The same process holds for the Z and W chromosomes, but we focus here on male heterogametic systems. Degradation of the Y chromosome can lead to pseudogenization and gene loss resulting in females (XX) having twice as many copies of genes on the sex chromosome compared with males (XY) (Wright et al. 2016). Because genes are normally expressed similarly from both copies of a chromosome, the degradation of the Y means that males only have half the expression of X-linked loci (Ohno 1967; Mank 2013), leading to a dosage imbalance with expression of genes on the autosomes. To resolve this issue, many organisms have evolved mechanisms to equalize expression levels of these sex chromosome genes, known as dosage compensation (Ohno 1967; Mullon et al. 2015). Dosage compensation mechanisms are thought to evolve rapidly in tandem with Y degradation (Ohno 1967), however, the majority of sex chromosomes with complete dosage compensation are relatively old, making it difficult to determine if dosage compensation can evolve in rapidly diverging sex chromosome systems (Mank 2009).

Significance

In species with XY sex chromosomes where the Y has degenerated, females (XX) have twice as many copies of X-linked genes compared with males (XY), leading to unbalanced expression between the sexes. Theory predicts that dosage compensation should evolve rapidly as X and Y chromosomes diverge, but examples of complete dosage compensation in recently diverged sex chromosomes are scarce, making this theory difficult to test. Across Poeciliid species the X and Y chromosomes have recently diversified. Here, we find complete dosage compensation evolved rapidly as the X and Y diverged in the common ancestor of Poecilia parae and P. picta, supporting that novel dosage compensation mechanisms can evolve rapidly in tandem with diverging sex chromosomes. These data confirm longstanding theoretical predictions of sex chromosome evolution.

Dosage compensation can either act by modifying expression on a gene-by-gene basis or globally by modifying expression along the entire chromosome (Graves 2016). Complete dosage compensation is achieved when gene expression is completely restored in the heterogametic sex, whereas incomplete dosage compensation is when expression is only partially restored (Mank 2013). Complete dosage compensation is predicted to arise for sex chromosomes that are rapidly diverging and experiencing extensive gene loss or pseudogenization, and has been more commonly found in male-heterogametic systems (XY) (Wilson Sayres and Makova 2011; Mullon et al. 2015). The most well characterized example for the rapid evolution of complete dosage compensation is in Drosophila where complete dosage compensation followed the emergence and divergence of a new XY sex chromosome system (Marín et al. 1996). The emergence of dosage compensation on neo-sex chromosomes in Drosophila is the result of evolution coopting extant dosage compensation mechanisms that predate the origin of the Drosophila genus (Marín et al. 1996). Although dosage compensation can clearly evolve rapidly, it is unknown if complete dosage compensation can evolve rapidly when it is not present in close relatives.

Fish exhibit a high rate of sex chromosome turnover, and although there are some species with incomplete dosage compensation (eg. sticklebacks, flatfish, and rainbow trout) (Shao et al. 2014; White et al. 2015; Hale et al. 2018) complete dosage compensation appears to be rare. We recently identified the first example of complete dosage compensation in a fish; Poeciliapicta. P.picta is a close relative to the common guppy (P.reticulata) and Endler’s guppy (P. wingei) and shares the same XY sex chromosome system as P. reticulata and P. wingei that originated 18.48–26.08 Ma (Rabosky et al. 2018; Darolti et al. 2019) (fig. 1). In P. reticulata and P. wingei, the X and Y have remained largely homomorphic, with little evidence of gene loss on the Y. With the maintenance of gene activity on the Y in P. reticulata and P. wingei, there is no need for dosage compensation, and indeed there is no evidence of dosage compensation in these species (Darolti et al. 2019). However, since the split of P. picta from P. reticulata and P. wingei ∼18.4 Ma (Rabosky et al. 2018) the P. picta Y has diverged substantially from the X across nearly the entire chromosome (∼90%), and a mechanism for complete dosage compensation has evolved (Darolti et al. 2019).

Fig. 1.

A phylogeny of Poecilia species depicting the ∼3.7-Myr interval that is the most parsimonious timeframe for when the Y chromosome degenerated and dosage compensation arose (orange branch) in the common ancestor of P. picta and P. parae. Orange species names indicate X and Y chromosomes are substantially diverged and have complete dosage compensation (P. picta—Darolti et al. 2019, P. parae—this study). Purple species names indicate dosage compensation has been explicitly tested but found to be entirely lacking (Darolti et al. 2019). Gray species names indicate dosage compensation has not yet been directly tested. Divergence times are labeled in gray at the relevant nodes. The phylogeny and divergence times are taken from The Fish Tree of Life (Rabosky et al. 2018), and is redrawn from Sandkam et al. (2021).

Here, we take a comparative approach to narrow the timing of the evolution of complete dosage compensation by testing for dosage compensation in P. parae, a sister taxon to P. picta (fig. 1). We recently characterized the sex chromosomes of P. parae, including five discrete Y haplotypes that control the five male morphs of this species (Sandkam et al. 2021). Importantly the pattern and extent of XY divergence is highly conserved across all five P. parae male morphs and is also conserved between P. parae and P. picta (Sandkam et al. 2021). This conserved pattern of XY divergence in all five P. parae male morphs and between P. parae and P. picta indicates that XY divergence likely occurred in the common ancestor of P. parae and P. picta during the ∼3.7 Myr interval spanning the split of the P. picta–P. parae common ancestor from P. reticulata ∼18.4 Ma, and prior to the split of P. picta and P. parae from each other ∼14.7 Ma (Rabosky et al. 2018). Therefore, if P. parae also has complete dosage compensation, it is most parsimonious that dosage compensation evolved rapidly in tandem with XY chromosome divergence over a period of less than 3.7 Myr (fig. 1).

Results

Characterization of Dosage Compensation in P. parae

To test whether complete dosage compensation evolved rapidly (over ∼3.7 Myr) in the common ancestor of P. picta and P. parae, we performed RNA-seq on muscle tissue from males and females of P. parae. There are five discrete Y haplotypes in P. parae that segregate with the five different male morphs (immaculata, yellow melanzona, blue melanzona, red melanzona, and parae morphs). Importantly, these five P. parae Y haplotypes emerged after X–Y recombination was halted in the common ancestor of P. picta and P. parae, ∼18.4 Ma (Lindholm et al. 2004; Sandkam et al. 2021). Therefore, if complete dosage compensation evolved in the common ancestor of P. picta and P. parae, we would expect to see dosage compensation in all male morphs as well. To assess this, we tested for differences in expression from the X and Y chromosome in three of the five male morphs (yellow melanzona, blue melanzona, and parae, hereafter referred to as yellow, blue, and parae males) by aligning reads to the female genome assembly from Sandkam et al. (2021). It is worth noting that all five Y haplotypes show similar patterns of divergence from the X (Sandkam et al. 2021), and so the three morphs we assessed here are indicative of the species as a whole.

For genes that are equally expressed from both sex chromosomes we expect to see a similar proportion of transcripts expressed from each sex chromosome. To test this, we first identified heterozygous transcripts (allowing us to distinguish expression between chromosomes). As expected, we found that heterozygosity did not differ between males and females for transcripts from the autosomes (total 38,986 autosomal transcripts; 17% heterozygous in males and 14% heterozygous in females; Wilcoxon rank sum P value = 0.4523). Nor did the heterozygosity differ between the transcript from the autosomes and those from the sex chromosomes in females (autosomes 17%, sex chromosome 12%; Wilcoxon rank sum P value = 0.0636). Yet in each of the male morphs tested, heterozygosity was substantially lower for transcripts from the sex chromosomes (1% of 1,349) compared with heterozygosity of transcripts from the autosomes (Wilcoxon rank sum P value = 0.0002). Heterozygosity of transcripts from the sex chromosomes was also substantially lower for males (1%) compared with females (12%) (Wilcoxon rank sum P value = 0.0004). This dramatic decrease in transcript heterozygosity on the sex chromosomes in males could be driven by high levels of hemizygosity, as is expected from our previous work indicating highly degraded Y chromosomes in this clade (Darolti et al. 2019; Sandkam et al. 2021).

We then compared the major allele ratios for heterozygous transcripts. We found that the major allele ratio for X-linked genes in females is close to 0.5 indicating roughly equal expression from both copies of the X with some variance due to cis-regulatory variation. In males we found a significant increase (Wilcoxon rank sum P value <0.001) in the major allele ratio compared with females and that the major allele ratio is close to 1 in each of the male morphs. In contrast, autosomal genes are equally expressed from both chromosomes in both sexes (fig. 2, Wilcoxon rank sum P value = 0.1457). These results are consistent with the notion that for sex-linked genes that remain heterozygous in males, gene activity has been greatly reduced from the Y gametolog and expression is primarily produced from the X.

Fig. 2.

Distribution of major allele ratios presented as kernel density plots for expressed genes on the sex chromosome (A) and the autosomes (B) for females and three male morphs of Poecilia parae indicate loss of expression for genes encoded on the Y. A major allele ratio of 0.5 indicates equal expression from both copies of a chromosome, whereas shifts toward 1 indicate expression predominantly comes from just one copy. Vertical dashed lines are median major allele ratio values. (C) Despite loss of expression from the Y, expression levels (log2 counts per million reads mapped [CPM]) of sex chromosome genes (n = 1,242) do not differ from the autosomes (n = 35,839) for any of the male morphs. (D) Male:female expression ratios for genes that exhibit allele-specific expression (ASE) are not different from male:female expression ratios of autosomal genes, demonstrating that a loss of expression from the Y chromosome in males does not result in reduced expression. The horizontal dashed line represents equal expression between males and females. Colors are consistent in all panels and denote sex and/or male morph. Gray, female; blue, blue male morph; yellow, yellow male morph; purple, parae male morph. Data in panes (C) and (D) are presented as box and whisker plots where the horizontal line is the median, the box denotes the 25th and 75th percentile, and “whiskers” are 1.5 times the interquartile range. Numbers above the X axis denote the number of transcripts. (E) Distribution of genes with allele-specific expression (ASE) along the male X chromosome (chromosome 8). Gene locations are demarcated by green lines. The pseudoautosomal region (PAR) is in gray.

To determine whether Y degeneration has been coupled with X chromosome dosage compensation, we compared average expression for genes from the X chromosome to the autosomal gene expression levels in both sexes. We found expression of sex chromosome genes was not different from autosomal genes in any of the male morphs (Wilcoxon rank sum yellow P value = 0.3544, blue P value = 0.745, parae P value = 0.121), or between genes on the X chromosome or autosomes in females (Wilcoxon rank sum P value = 0.4568) (fig. 2). There was no systemic difference in the distribution of the M:F ratio for genes on the sex chromosome (supplementary fig. S1, Supplementary Material online). The distribution of the M:F ratio was not standard in the yellow or blue male morph (yellow mean = −0.14±0.51 P < 0.001, t-test P value <0.001; blue mean −0.13±0.54, t-test P value <0.001) but did not significantly differ from zero in the parae male morph (parae mean = −0.005±0.75, t-test P value = 0.84) (supplementary fig. S1, Supplementary Material online).

Moreover, we observe no significant differences in the M:F expression ratio for sex-linked genes with an allele-specific expression (ASE) pattern compared with the M:F ratio for autosomal genes in any of the three male morphs (Wilcoxon rank sum yellow P value = 0.11, blue P value = 0.105, parae P value = 0.107) (fig. 2). We also found no consistent differences in the distribution of the M:F expression ratios for genes with ASE on the sex chromosome (yellow mean = −0.16±0.39, P value = 0.00; blue mean = −0.16±0.44, P value = 0.00; parae mean = 0.11±0.62, P value = 0.21) (supplementary fig. S2, Supplementary Material online). ASE genes were distributed along the entire X chromosome providing further support for global dosage compensation for genes along the entire sex chromosome and not regional dosage compensation or dosage compensation on a gene-by-gene basis (fig. 2). Taken together, these data indicate complete dosage compensation in P. parae.

Discussion

Evidence for Upregulation of the X Chromosome in Males

In this study, we identify a complete dosage compensation system in P. parae. In general, there are two ways in which complete dosage compensation has been observed in XY systems. In eutherian and marsupial mammals, one of the two X chromosomes is silenced in females. Although this balances sex chromosome gene expression between males and females, it does not address expression differences between X-linked and autosomes genes. In fact, X inactivation in females means that both sexes on average express X-linked genes less than the autosomal average, and only dosage sensitive genes on the X are upregulated in both sexes to counter this (Pessia et al. 2012). Alternatively, in Drosophila (Marín et al. 2000), and Anolis (Marin et al. 2017), dosage compensation is achieved by doubling the expression of genes on the X chromosome in males.

In P. parae, we found very few sex chromosome genes with heterozygous expression (1%) which strongly suggests that most of the sex chromosome genes in males are hemizygous. Despite evidence for hemizygous expression for sex chromosome genes in males, genes on the sex chromosome are expressed at similar levels in males and females, and between the sex chromosome and autosomes in males. Moreover, heterozygous genes exhibit ASE of a single allele in males, and this pattern of ASE is not observed in females or for autosomal genes in males. Furthermore, ASE genes in males are distributed along the entire X chromosome instead of being clustered in a specific region (fig. 2). Taken together, these data strongly support widespread gene loss on the Y chromosome and an upregulation of genes along the entire X chromosome in males indicating the evolution of a complete dosage compensation mechanisms in P. parae. Although it is possible that some of these ASE genes are Y gametologs expressed at higher levels than the X, we would expect the opposite for most loci given the widespread loss of gene activity on the Y that occurs with XY divergence. Furthermore, Y upregulation would result in male-biased expression when combined with X chromosome dosage compensation, and we instead observe similar male and female expression levels for ASE genes. Taken together, these data suggest that complete dosage compensation in P. parae is more similar to dosage compensation in Drosophila and Anolis, where genes on the X are hyper expressed in males. This provides an excellent avenue to explore the mechanisms controlling expression across entire chromosomes.

Rapid Evolution of Dosage Compensation

Although theory suggests complete dosage compensation should evolve rapidly in tandem with Y degradation (Ohno 1967), gene expression studies in nonmodel systems with heteromorphic sex chromosomes have demonstrated that complete dosage compensation is actually quite rare and is not a guaranteed outcome of sex chromosome evolution (Mank et al. 2011). Until the recent characterization of complete dosage compensation in P. picta, complete dosage compensation in vertebrates has been observed in a limited number of lineages with relatively ancient (>160 Myr) sex chromosomes (Marin et al. 2017). The age of these systems makes it difficult to refine estimates for the speed at which complete dosage compensation can arise.

Within the family Poeciliidae the subgenus Lebistes is particularly well suited to address this question as it contains several species with characterized sex chromosomes including P. reticulata, P. wingei, P. picta, and P. parae which have recently diversified (Darolti et al. 2019). Although it is possible that the sex chromosomes in each of these lineages has evolved on the same linkage group independently, there is strong evidence that all Lebistes share the same sex chromosome system which originated in a common ancestor 18.48–26.08 Ma (Rabosky et al. 2018; Darolti et al. 2019). Despite sharing the same XY system, the extent of Y degradation differs dramatically, and although the Y is largely intact in P. reticulata and P. wingei, it is highly degraded in P. picta and P. parae (Darolti et al. 2019; Sandkam et al. 2021). Without gene loss, there is no selective pressure to evolve dosage compensation, thus it is not surprising that dosage compensation was not found in either P. reticulata and P. wingei (Darolti et al. 2019), where there is little evidence of decreased gene activity from the Y chromosome.

In contrast to P. reticulata and P. wingei, the Y chromosomes in P. picta and P. parae exhibit substantial divergence along the entire chromosome (Darolti et al. 2019; Sandkam et al. 2021). Similarities in the boundary of the pseudoautosomal region, the extent of Y divergence, and shared Y-specific K-mers present in both P. picta and P. parae (Sandkam et al. 2021) suggest a shared history of X and Y divergence in a common ancestor that predates the divergence of these two species. Here, we present evidence for complete dosage compensation in multiple morphs of P. parae. It is possible that degeneration and complete dosage compensation independently evolved twice in these two closely related species, or even within each of the morphs separately and independently of P. picta. However, the most parsimonious explanation for the similarities in the patterns of Y degeneration between these two species is that dosage compensation evolved once over a period of less than 4 Myr in their common ancestor (fig. 1).

The rapid evolution of complete dosage compensation is typically associated with the recruitment of extant dosage compensation machinery, such as ancestral or pre-existing dosage compensation mechanism (Marín et al. 1996; Marin et al. 2017), or through the recruitment of sex-specific autosomal transregulatory mechanisms to compensate for the degeneration of cis-regulatory elements (Lenormand et al. 2020). In fishes, complete dosage compensation is rare, which may be the result of frequent sex chromosome turnover and a paucity of heteromorphic sex chromosomes (Ashman et al. 2014) that makes complete dosage compensation evolution unlikely (Vicoso 2019). As such dosage compensation in fish is frequently accomplished on a gene-by-gene basis and remains overall incomplete (Shao et al. 2014; White et al. 2015; Darolti et al. 2019) with the exception of P. picta (Darolti et al. 2019) and now also P. parae. Further work elucidating the mechanism of X chromosome dosage compensation in P. picta and P. parae will provide novel insights in the evolution of dosage compensation mechanisms.

Materials and Methods

RNA Isolation and Sequencing

Animals used in this study were collected in Spring 2019 from natural populations in Suriname and brought to the University of British Columbia (Vancouver, BC, Canada) aquatics facility, where they were kept in 20 L glass aquaria on a 12:12 day:night cycle at 26 °C and 13ppt salinity (Instant Ocean Sea Salt), and fed Hikari Fancy Guppy pellets supplemented with live brine shrimp daily. Individuals were euthanized using a lethal overdose of MS-222 and muscular tail tissue was taken from the anal pore to the base of the pectoral fin. RNA was immediately isolated using RNeasy spin columns with on-column DNase treatment (Qiagen) following the manufacturer’s recommended protocol. RNA was isolated from three females and from three of each of the male morphs (12 samples in total). Library preparation and 100-bp paired-end sequencing was performed on an Illumina NovaSeq 6000 at McGill University and the Génome Québec Innovation Centre. Adaptor sequences were removed and reads were quality filtered and trimmed using trimmomatic (v0.36) using a sliding window of four bases and a minimum Phred score of 15. Reads with leading and trailing bases with a Phred score <3 were also removed. Sequencing libraries consisted of ∼88 million reads (supplementary table S1, Supplementary Material online).

Transcript Alignment and Filtering

Reads were aligned to a female P. parae genome assembly (Sandkam et al. 2021) using the two-pass method for STAR align (v2.7.2) (Dobin et al. 2013) with default settings and the maximum number of alignments per read (–outFilterMultimapNmax) set to one. Alignments were sorted by coordinate and converted to BAM format using SAMtools (v1.9) (Li et al. 2009) using default settings. To find the full list of nonredundant P. parae transcripts we generated GTF files for each individual using StringTie (v1.3.6) (Pertea et al. 2015) using default settings, then merged all GTF files. To remove noncoding RNA (ncRNA) we first compiled a database of all ncRNAs in reference genomes of close relatives on Ensemble: Poecilia formosa (PoeFor_5.1.2), Oryzias latipes (ASM223467v1), Gasterosteus aculeatus (BROAD S1), and Danio rerio (GRCz11). We then used BlastN v2.8.1 (Altschul et al. 1990) with -evalue 10e-10 to identify putative ncRNAs. Using this approach, we identified 514 putative ncRNAs (0.8% of all transcripts) that we removed from our P. parae transcript database.

Allele-Specific Expression

To ensure our results are comparable to our previous results in P. picta we followed the same pipeline to identify ASE (Darolti et al. 2019). In short, for each sex and morph, we identified SNPs separately using SAMtools mpileup (v1.9) and VarScan (v2.4.3) (Koboldt et al. 2012) with parameters –min-coverage 2, –min-ave-qual 20, –min-freq-for-hom 0.90, and excluding triallelic SNPs. We then filtered SNPs for a minimum site coverage of 15 to account for sequencing errors, and used a variable coverage filter to account for potential effects of sequencing errors due to variable coverage levels (an error rate of one in 100 and a maximum coverage for a given site of 100,000) (Quinn et al. 2014). We then removed SNP clusters of more than five SNPs in 100 bp window to limit potential biases from read assignments to a single reference sequence (Stevenson et al. 2013). We then used a Wilcoxon rank sum test to test for differences in the median major allele frequency between males and females for genes on the sex chromosome or genes on the autosomes.

Expression Level

We extracted read counts using the featureCounts program from the Subread package (Liao et al. 2014) and the ncRNA filtered GTF file described above. We then used the RLE normalization method in edgeR v3.24.3 (Robinson et al. 2010) to control for library size variation. To account for biological and technical variation in sequencing and read-mapping, genes with low expression (less than 5% of the mean) and the top 0.01% of expressed genes were removed from the data set. We then used a Wilcoxon rank sum test to compare expression levels, in counts per million, between groups using the wilcox.test() function in R (P < 0.05). For the analysis of M:F expression ratios for genes with ASE on the sex chromosome, we subsampled the M:F ratios of autosomal genes based on the number of ASE genes in each morph 1,000 times. We then performed a Wilcoxon rank sum test for each subsampled group and calculated the mean P value for each morph.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.