Genomic Loss and Silencing on the Y Chromosomes of Rumex.

Across many unrelated lineages of plants and animals, Y chromosomes show a recurrent pattern of gene degeneration and loss, but the relative importance of inefficient selection, adaptive gene silencing, and neutral genetic drift in causing degeneration remain poorly understood. Here, we use next-generation genome and transcriptome sequencing to investigate patterns of ongoing Y chromosome degeneration in two annual plant species of Rumex (Polygonaceae) differing in their degree of degeneration and sex chromosome heteromorphism. We find evidence for both gene loss as well as silencing in these young plant sex chromosomes. Our analyses revealed significantly more gene deletion relative to silencing in R. rothschildianus, which has had a larger nonrecombining region for a longer period than R. hastatulus, consistent with this system being at a more advanced stage of degeneration. Intra- and interspecific comparisons of genomic coverage and heterozygosity indicated that loss of expression precedes gene deletion, implying that the final stages of mutation accumulation and gene loss may often occur neutrally. We found no evidence for adaptive silencing of genes that have lost expression. Our results suggest that the initial spread of deleterious regulatory variants and/or epigenetic silencing may be an important driver of early degeneration of Y chromosomes.

Introduction

Sex chromosome evolution is characterized by recurrent patterns of gene degeneration and loss. One of the most extreme examples of this is the nonrecombining human Y chromosome in which only 3% of genes remain functional, despite having ancient homology with the X chromosome (Skaletsky et al. 2003; Wilson Sayres and Makova 2013; Bellott et al. 2014). Even in much younger sex chromosome systems, nonrecombining sex chromosomes consistently fix more nonsynonymous amino acid changes (Bergero and Charlesworth 2011; Chibalina and Filatov 2011; Gschwend et al. 2012; Papadopulos et al. 2015; White et al. 2015; Crowson et al. 2017), exhibit reduced codon bias (Hough et al. 2014; Singh et al. 2014), have higher numbers of insertions and deletions (Bergero et al. 2008, 2015; White et al. 2015), and have high rates of accumulation of transposable elements (reviewed in Kejnovsky et al. 2012). Despite these parallel patterns across unrelated taxa, the relative importance of inefficient selection, positive selection, and/or neutral drift in driving Y chromosome degeneration remains poorly understood.

Linked selection can drive Y degeneration through the recurrent accumulation of deleterious mutations (Charlesworth 1996). When the rate of recombination between sites under selection is low, linked sites can interfere with the population dynamics of neighboring sites, a process known as Hill–Robertson interference (Hill and Robertson 1966). Selection in regions with low recombination rates is less likely to remove slightly deleterious mutations than in regions with higher recombination rates in which the segregation of selected sites is more independent. The intensity of selective interference increases as more sites under selection are locked together in the nonrecombining genomic region (Kaiser and Charlesworth 2009). Even strongly deleterious mutations may accumulate through Muller’s ratchet, the stochastic loss of individuals with the fewest number of deleterious mutations, and this effect should be particularly prominent in species with low effective population size (Ne) (Charlesworth 1996). In addition to interference due to purifying selection, selective interference can be caused by selective sweeps, where sites subject to positive selection can drag along linked deleterious variants to fixation (Maynard Smith and Haigh 1974). Regions of low recombination have been shown to exhibit elevated levels of deleterious mutation and reduced neutral diversity, consistent with model predictions (Betancourt and Presgraves 2002; Cutter and Payseur 2013; Ellegren and Galtier 2016; Elyashiv et al. 2016). Therefore, inefficient selection is considered a major driver of mutation accumulation on the nonrecombining Y chromosome.

Alternatively, differential positive selection could play an important role in Y chromosome evolution. Beneficial alleles will fix at a higher rate in recombining regions, because their fate is uncoupled from nearby deleterious mutations (Peck 1994). As a result, positive selection should proceed at a much higher rate on the recombining X chromosome compared with the Y, due to reduced selective interference (Orr and Kim 1998). Given this differential positive selection, silencing of the less-fit Y homologue may be beneficial. Adaptive silencing is likely to be especially strong on Y homologues of X-linked genes that have recently undergone positive selection (see Crowson et al. 2017).

It is also possible that neutral genetic drift plays a role in sex chromosome evolution, particularly if gene silencing and/or dosage compensation occur early in the process. Genes that lose expression can degenerate neutrally, as selection will no longer act on them. In this manner, loss of expression can amplify the speed of degeneration of the nonrecombining region. One way in which silencing may arise early would be through repetitive element invasion, followed by heterochromatin spreading in cis, and finally degeneration by neutral drift, as has been proposed for neo-Y chromosome systems in Drosophila (Zhou and Bachtrog 2012a; Zhou et al. 2013).

These three evolutionary mechanisms acting on the Y chromosome are not mutually exclusive, and disentangling their relative importance remains a challenge. On the neo-Y chromosome of Drosophila miranda there is a general association between loss of gene expression and the rate of molecular evolution (Zhou and Bachtrog 2012b). These results alone cannot distinguish whether silenced genes accumulated missense mutations or whether pseudogenes were actively silenced; however, the observation of a subset of functional genes that have been silenced suggests that heterochromatin formation may sometimes precede gene inactivation (Zhou et al. 2013). Consistent with this finding, in the younger neo-Y chromosome of D. albomicans gene silencing occurs despite little evidence of degeneration (Zhou and Bachtrog 2012a). Collectively, these results from Drosophila suggest that in neo-Y sex chromosome systems neutral evolution following loss of expression may play a significant role in Y chromosome degeneration.

Changes in gene expression are likely to take place early in the evolution of a nonrecombining sex chromosome, making young sex chromosomes valuable for investigating the relative importance of degeneration versus adaptive silencing and changes in gene expression. Furthermore, young plant sex chromosomes provide an interesting comparison to animal systems because the rate of sex chromosome degeneration can be slowed or halted by haploid selection in the pollen (gametophytic) stage of the life cycle (Haldane 1933; Chibalina and Filatov 2011; Scott and Otto 2017). Despite this potential difference, evidence to date from Silene latifolia suggests that gene degeneration may be fairly extensive (Bergero et al. 2015; Papadopulos et al. 2015). In this system, weakly or nonexpressed Y-linked genes were more than twice as likely to be interrupted by premature stop codons than genes with high expression on the Y (Papadopulos et al. 2015). In Silene, gene loss appears to be more prevalent than gene silencing, and there was no evidence that younger evolutionary strata showed proportionally more silenced genes, arguing against extensive adaptive gene silencing or early heterochromatin-based degeneration (Bergero et al. 2015).

Rumex provides excellent opportunities to investigate the early stages of Y chromosome degeneration in plants. The genus includes several clades of dioecious species with sex chromosomes, and dioecy is thought to have evolved ∼15–16 Ma (Navajas-Pérez et al. 2005). The identification of sex-linked genes from controlled crosses in R. hastatulus and R. rothschildianus using RNA sequencing has revealed wide variation in the proportion of genes that have been silenced or lost from the Y chromosome (Hough et al. 2014; Crowson et al. 2017). Rumex hastatulus is of interest because it has a polymorphic sex chromosome system: although timing estimates remain very approximate, males from west of the Mississippi River have an ∼9–16 Ma XY system (Crowson et al. 2017), whereas those from the southeastern United States possess an XYY karyotype postulated to have arisen ∼600,000 years ago from a Robertsonian fusion between an ancestral autosome and the X chromosome shared with western populations (Smith 1964; Quesada del Bosque et al. 2011). Estimates from leaf tissue data found a 24–28% loss in gene expression from the ancestral R. hastatulus Y chromosome, and <8% from the neo-Y (Hough et al. 2014). In contrast, R. rothschildianus has an older 8–13 Ma XYY system, with the second Y chromosome thought to be derived from an X-autosome fusion shared with a large subclade of dioecious Rumex species (Navajas-Pérez et al. 2005; Crowson et al. 2017). Based on cytogenetic studies, Y chromosomes of Rumex species in this subclade are more heterochromatic, suggesting the possibility of a greater degree of degeneration (Degraeve 1976; Navajas-Pérez et al. 2006; Cunado et al. 2007). Consistent with this, recent estimates in R. rothschildianus suggest upward of 90% of genes have lost Y expression, implying much more extensive degeneration, possibly due to greater selective interference caused by the presence of a larger nonrecombining region over a longer period (Crowson et al. 2017).

Recent studies of R. rothschildianus and R. hastatulus suggest that they may have independently evolved sex chromosomes, including neo-sex chromosomes, but there is currently no clear evidence suggesting evolutionary strata (Navajas-Pérez et al. 2005; Hough et al. 2014; Crowson et al. 2017). Thus, neo-sex chromosomes in Rumex have the potential to provide insight into the early events of Y chromosome degeneration because of the occurrence of distinct recombination suppression events between ancestral and neo sex-chromosomes. Furthermore, variation in the number of sex-linked sites offers opportunities to evaluate different amounts of interference selection among Y-linked genes, which may lead to different stages of Y chromosome degeneration (Crowson et al. 2017).

Recent analyses of R. rothschildianus using molecular evolutionary approaches reported that X-linked genes without Y expression show elevated rates of molecular evolution, lower ancestral gene expression, and fewer signals of positive selection (Crowson et al. 2017). These results suggested that inefficient selection may be the primary driver of Y degeneration in Rumex; however, these analyses involved transcriptome data and therefore could not distinguish between genes that have been silenced versus those that have been deleted from the Y chromosome. Furthermore, it is also possible that positive selection for gene silencing has occurred for a subset of genes, and these analyses, which focused on patterns of molecular evolution on a large timescale, did not consider the possibility of recent positive selection on X-linked genes with silenced Y homologues.

Using progeny arrays, Hough et al. (2014) and Crowson et al. (2017) used RNAseq to identify genes with single nucleotide polymorphism (SNP) segregation consistent with hemizygous X-chromosome expression (hereafter hemizygous), and another set with SNP segregation consistent with expression of distinct alleles on the X- and Y-chromosomes (hereafter X&Y). By complementing their SNP segregation analysis with new genomic read data, we can distinguish whether genes with hemizygous expression are silenced or effectively deleted on the Y. First, whereas hemizygous genes which have lost their Y copy will have half the coverage in males compared with females (Vicoso et al. 2013), silenced genes should have equal genomic coverage for each sex. Furthermore, genes in males that have lost their Y copy should have no heterozygosity, whereas genes with an X and a Y copy should have elevated levels of heterozygosity in males as the homologs diverge due to the loss of recombination (Vicoso et al. 2013). Thus, genes that have been silenced should have no heterozygosity in the transcriptome data (where there is no Y allele), but should have elevated amounts of heterozygosity for genomic reads (where the X and Y have diverged).

Here, we compare the results of short-read whole genome sequencing of R. rothschildianus and of XY and XYY karyotypes of R. hastatulus to explore the causes and consequences of changes in gene expression and gene loss in the early stages of degeneration of nonrecombining sex chromosomes. By integrating genome and transcriptome data with population genetic data, our analyses allow the distinction between gene silencing and gene loss for hemizygous genes, as well as for explicit tests of contemporary positive selection on the X-linked homologues of recently silenced Y-linked genes. Our results indicate that gene silencing is likely to precede deletion, although we found no evidence suggesting that gene silencing is adaptive. Our results suggest that the final stages of degeneration may proceed neutrally, although Hill–Robertson interference probably plays an important role in the early degeneration of sex chromosomes.

Materials and Methods

Previous Identification of Sex-Linked Genes

Leaf samples from parents and six male and six female progeny for both subspecies of Rumex hastatulus were previously sequenced using Illumina RNAseq as described in Hough et al. (2014). Hemizygous and X&Y genes were identified on the basis of SNP segregation patterns requiring at least four segregating SNPs, with putative sex-linked genes also validated using polymorphism data from transcriptomes sequenced from a range-wide sample. Similarly, leaf samples from two independent crosses, with parents and six males and six female progeny sequenced for one cross, and four males and six females for the other cross were previously obtained by Crowson et al. (2017). Sex-linked genes were similarly identified in this study, with the requirement of consistent SNP segregation patterns in both families.

DNA Extraction and Sequencing

We extracted DNA from leaf tissue of a male and female plant from R. hastatulus from both Marion, South Carolina and Wesley Chapel, Texas (see Pickup and Barrett 2013 for locality and population information), representative of the XYY and XY karyotypes, respectively. We obtained genomic extractions using the DNeasy Plant Mini Kit (QIAGEN GmbH, Hilden, GERMANY) according to manufacturer’s protocol. Similarly, we obtained DNA from a male and female plant of R. rothschildianus obtained from crosses derived from a sample originally from Tel Aviv Botanical Garden, Israel and described in Crowson et al. (2017).

We obtained Illumina paired-end 250-bp libraries for all samples from the Genome Quebec Innovation Centre, except the Texas male, which was sequenced at The Center for Applied Genomics in Toronto, ON. XYY male and female libraries were multiplexed and run on the same two lanes of Illumina TruSeq HT using the Rapid Run Mode, whereas the female XY and the R. rothschildianus male and female libraries were run separately from the XYY libraries, also on the same two lanes of Illumina TruSeq HT using the Rapid Run Mode. The XY male was run on one lane of Illumina HiSeq X. To assess batch effects caused by runs over different flow cells and sequencing centers, we compared various read and mapping quality metrics (read count, % of reads mapped, median coverage, GC content, Transition/Transversion ratio, SNP count, heterozygous/homozygous SNP ratio, insert size mean, mapping quality [MapQ score] mean), following the approaches of Tom et al. (2017). These results are outlined in supplementary table S1, Supplementary Material online. A Principle Component Analysis (PCA) found differentiation along expected group clusters (species and karyotype), but not between sequencing batches (supplementary fig. S1, Supplementary Material online). We attribute these major group effects to between-species and karyotype differences in genome structure and base composition, rather than technical artifacts. When considering only quality metrics (reads, coverage, MapQ, and % mapped), the PCA explains 1.4% of the variance, suggesting few batch effects.

Read Alignment and Processing

Reads with 50% of bases with a quality score <23, or with <10% of nonambiguous bases, were removed from further analysis using an in-house script (https://github.com/cafeblue/wei_script/blob/master/clean_B_N_HiSeq.pl; last accessed December 6, 2017). We aligned R. hastatulus genomic reads to the female XYY R. hastatulus reference transcriptome assembly (see Hough et al. 2014) and R. rothschildianus to the female R. rothschildianus reference transcriptome assembly (see Crowson et al. 2017) according to GATK best practices (https://software.broadinstitute.org/gatk/best-practices/bp_3step.php?case=GermShortWGS; last accessed December 6, 2017). Briefly, alignments were first conducted using BWA mem, allowing multiple mapping but not discordant mapping, and using a mismatch penalty score of 4, a minimum seed length of 19, a band-width of 100, and a gap penalty of 6 (Li and Durbin 2009). This alignment was followed by further mapping using default STAMPY parameters, with the “–bamkeepgoodreads” option on, set to a substitution rate of 0.001, and disallowing either multiple mapping, or discordant mapping (Lunter and Goodson 2011). The STAMPY step enables mapping of divergent reads, which is important in order to minimize mapping bias due to SNP divergence between X and Y. Separate lane alignment files were processed and merged using PICARD’s AddOrReplaceReadGroups followed by MarkDuplicates, with default parameters (https://broadinstitute.github.io/picard/; last accessed December 6, 2017).

Our estimates of the proportion of genes silenced and deleted from the Y chromosome are robust to a range of MapQ score cutoffs [0, 10, 30, 50] (supplementary table S2, Supplementary Material online). For coverage analysis (see below), we determined coverage from bam files using samtools (Li et al. 2009) and normalized coverage by dividing read depth by locus length. We adjusted read depth by dividing by each sample’s median coverage. We excluded loci with <5 reads or >2 SD coverage in either males or females, to retain only single copy genes. For heterozygosity analysis (see below), we called SNPs using GATK Haplotypecaller (McKenna et al. 2010). SNPs were filtered using GATK VariantFiltration with a cutoff site quality score of >50.

Identification of Silenced or Deleted Genes

Previous work in Rumex identified assembled transcripts as autosomal, X&Y expressed, or hemizygous expressed based on RNAseq SNPs segregating in a family pedigree (Hough et al. 2014; Crowson et al. 2017). We used these subsets of assembled transcripts for genomic status using the coverage and heterozygosity techniques developed by Vicoso et al. (2013). We performed univariate Gaussian mixed model analyses on male to female genomic coverage ratios using the R package mclust (Fraley and Raftery 2012). Genes in the cluster nearest to half genomic coverage in male as compared with female (log2=−1) were defined as deleted, whereas genes in the other cluster, near equal coverage in males and females, were defined as silenced. As further validation of the presence of both gene loss and silencing, we assessed heterozygous SNP density in X&Y and hemizygous genes in our genomic data. A subset of the previously identified hemizygous genes exhibited heterozygous SNP calls in the male transcriptome, possibly due to either genotyping error and/or chimeric transcripts; to simplify our analyses and increase stringency, transcripts with heterozygosity in the male RNAseq data were removed from this analysis. In our heterozygosity analysis, we compared genes with no heterozygous SNP calls in male genomic reads to those with heterozygous SNPs, as an alternative means of estimating the extent of deletion versus gene silencing.

Statistical Comparisons

To ascertain the timing of silencing and deletion in R. hastatulus, we quantified genes as “old hemizygous,” if they were hemizygous in both karyotypes—suggestive of gene loss in the common ancestor—or as “new hemizygous” if they were hemizygous in one karyotype but X&Y in the other. Genes with inconclusive segregation in either karyotype were removed from these analyses. New hemizygous genes from each karyotype could be grouped together because their genomic coverage ratio is with respect to the karyotype from which they were lost; however, old hemizygous genes can have distinct genomic coverage ratios in each of the two karyotypes. For old hemizygous genes, the results from each karyotype are included in table 1. In further analyses, we include only the results for the XYY karyotype for old hemizygous genes.

Table 1

Counts and Ratio of Silenced and Deleted Hemizygous Genes in Rumex hastatulus, Split by the Segregation Pattern of Their Homologs in the Other Karyotype as well as by Coverage Respective to Each Karyotype

	Hemizygous in Both			X&Y
	Coverage in XY	Coverage in XYY		Coverage in XY	Coverage in XYY
Silenced	1	2		31	27
Deleted	7	6		8	11
Silenced/Deleted	0.14	0.33		3.88	2.45

Polymorphism data of females from XY R. hastatulus populations used for π and Tajima’s D estimates are from Hough et al, and dn/ds ratios and PAML (Yang 1997) results for XYY R. hastatulus karyotype are from Crowson et al. (2017). A gene enrichment test was performed using BLAST2GO (Gotz et al. 2008).

Results

Genomic Coverage

At the genome-wide scale, R. rothschildianus and R. hastatulus showed no significant signals of gene loss (supplementary fig. S2, Supplementary Material online); in particular, we failed to observe a clear genome-wide bimodal distribution for log2 male/female coverage ratios. This pattern suggested that either there was no degeneration of the sex chromosomes in R. hastatulus or R. rothschildianus, or that patterns of degeneration were too restricted in the transcriptome as a whole to be clearly observable at the genome-wide scale. We therefore investigated the latter possibility by isolating sex-linked transcripts.

X&Y and Hemizygous Coverage

By using previously identified sex-linked genes based on SNPs from RNAseq data in pedigrees, we investigated whether our identified sex-linked transcripts showed evidence of effective deletion that was not observable using genome-wide data. As expected, X&Y genes showed a log2 male/female coverage ratio distribution centered around zero for both species (supplementary fig. S2, Supplementary Material online), indicating equivalent coverage in males and females. In contrast, hemizygous genes in R. hastatulus were significantly different from the X&Y coverage (Wilcoxon rank sum test, XY-karyotype: W = 7436.5, P < 2.2e-16; XYY-karyotype: W = 52,410, P < 2.2e-16) and exhibited a bimodal distribution of genomic coverage, suggesting the presence of both silencing and effective deletion on the Y chromosome for both sex chromosome karyotypes. Univariate Gaussian mixed model analyses estimated that 55% (54 of 98) and 70% (81 of 115) of the Y homologs of hemizygous genes were deleted in leaf tissue for the XY and XYY karyotypes, respectively (fig. 1). A χ2 contingency test indicated significantly more silencing relative to deletion in the XY system as compared with the XYY system (χ2= 4.72, P = 0.03). We also observed evidence for a greater fraction of gene deletion in the hemizygous gene set of the more degenerated R. rothschildianus (fig. 1). We estimate that for 92% (215 of 233) of hemizygous genes in R. rothschildianus, the Y homolog has been deleted as opposed to silenced, which is significantly greater than in either R. hastatulus karyotype (for XYY, χ2 =24.63, P = 6.96e-07).

. 1.

—Distribution of log2 of the male to female genomic coverage ratio across hemizygous genes in (a) Rumex hastatulus XY karyotype, (b) R. hastatulus XYY karyotype, and (c) R. rothschildianus. Sex linkage was determined using SNP segregation patterns. Bars represent histograms, whereas curves in read ratio plots represent kernel density estimates (d–f). Gene classification predicted using a univariate Gaussian mixed model.

Heterozygosity

Further evidence for silencing of a subset of hemizygous genes was obtained by comparing heterozygosity in R. hastatulus XY and XYY karyotypes, and in R. rothschildianus. As expected, X&Y genes showed elevated heterozygosity in the male for DNA reads (supplementary fig. S3, Supplementary Material online), due to sequence divergence between the sex chromosomes. In contrast, homologs of 50% of hemizygous genes were found to have heterozygous calls in the male genomic sequence (putatively silenced), whereas the remainder had no heterozygous calls in the male (putatively deleted) for the XY karyotype (fig. 2). This proportion increases to 56% with no genomic heterozygosity (putatively deleted) in the XYY karyotype (fig. 2). These proportions are not significantly different from the proportions predicted using coverage ratio (for the XY karyotype: χ2 = 0.16, P = 0.69). Consistent with this, genes in R. hastatulus with no male genomic heterozygosity also had a female biased coverage ratio, whereas genes with high male genomic heterozygosity had no coverage ratio bias (color in fig. 2). Overall, the patterns of heterozygosity we observed are consistent with the conclusion that hemizygous genes are a mixture of silenced and deleted genes. These results further support the conclusion that there is reduced expression and extensive silencing of Y-linked genes, in addition to gene loss.

. 2.

—Number of heterozygous DNA read SNPs in males versus females for hemizygous genes in (a) Rumex hastatulus XY karyotype, (b) R. hastatulus XYY karyotype, and (c) R. rothschildianus. Shapes colored according to log2 of male to female genomic coverage ratio, where yellow is sex-unbiased genomic coverage and blue is female biased.

Timing of Silencing versus Loss

Our interspecific comparisons indicated that hemizygous genes in the more degenerated system of R. rothschildianus have a greater proportion of gene deletion relative to gene silencing than in R. hastatulus. This difference suggests that the earlier stages of Y degeneration involve gene silencing before gene loss. To investigate this further, we examined cases in R. hastatulus where genes on the ancestral Y have recently lost expression in one of the two karyotypes, giving us a set of genes with very recent expression loss (hereafter young hemizygous). We compared our set of young hemizygous genes to the set of genes that could be identified as being hemizygous in both karyotypes (hereafter old hemizygous). By comparing whether these genes are more likely to be silenced or lost, we can establish whether early loss of expression is more likely to be due to silencing or gene loss. We hypothesized that the early stages of gene loss would involve silencing, followed later by deletion from the Y chromosome. Indeed, homologs of 80% of old hemizygous genes (hemizygous in both karyoptypes) are effectively deleted, whereas homologs of only 19% of young hemizygous genes (hemizygous in only one of the two karyotypes) show evidence of deletion (χ2 =4.2, P = 0.04) (fig. 3 and table 1).

. 3.

—Distribution of log2 of the male to female genomic coverage ratio across genes for (a) old and (b) young hemizygous genes of Rumex hastatulus. Bars represent histograms and curves represent kernel density estimates. (c and d) Gene classification predicted using a univariate Gaussian mixed model.

No Evidence of Active Silencing

To investigate evidence of adaptive silencing due to differential positive selection (Orr and Kim 1998), we asked whether genes with signals of positive selection on the X are more likely to be silenced on the Y. Recent analysis of R. rothschildianus using patterns of molecular evolution found no evidence for elevated signals of positive selection for hemizygous genes (Crowson et al. 2017). However, this study could not distinguish between genes that were silenced versus those that were deleted, and could not identify the possibility of recent lineage-specific positive selection. In our analysis, a Mann–Whitney U-test found no significant difference for π at synonymous sites between X&Y to hemizygous genes in the XY karyotype of R. hastatulus (U = 28,077, P = 0.52) or between silenced or deleted genes (U = 1031.5, P = 0.96). Similarly, we found no significant difference in Tajima’s D at synonymous sites between X&Y to hemizygous genes (U = 12,070, P = 0.68) or between silenced to deleted genes (U = 452, P = 0.63). Using the hermaphroditic R. bucephalophorus as an outgroup, we found no significant difference for dn/ds between silenced to deleted genes in the XYY karyotype of R. hastatulus (U = 66, P = 0.59). A search for enrichment in gene ontology families found no significant differences in functional enrichment for effectively deleted compared with silenced genes.

Finally, to determine whether dosage compensated genes were more likely to be silenced, we tested whether deleted and silenced genes showed significant differences in the ratio of male to female expression. A Mann–Whitney U-test found no significant difference between silenced and deleted genes for male to female expression in the XY karyotype (U = 1,040, P = 0.36) or XYY karyotype (U = 1418.5, P = 0.80) of R. hastatulus, or in R. rothschildianus (U = 2,400, P = 0.09) (supplementary fig. S4, Supplementary Material online).

Discussion

Sex chromosomes in plants and animals are derived from an autosomal pair that diverged and degenerated following the loss of recombination (Charlesworth 1996). In many ancient animal sex chromosome systems such as those of birds, mammals, and flies, the degeneration of the Y chromosome is extensive (Handley et al. 2004; Bellott et al. 2014; Vicoso and Bachtrog 2015). In flowering plants, however, heteromorphic sex chromosomes are relatively rare (Bachtrog et al. 2014). This observation has led some to propose that flowering plant sex chromosomes do not lose genes to the same extent as animal sex chromosomes because of substantial selection in the haploid gametophytic phase of the lifecycle (Haldane 1933; Nei 1970; Chibalina and Filatov 2011). However, the hypothesis that widespread degeneration of plant Y chromosomes is prevented by haploid selection is inconsistent with recent work from Silene latifolia, where the Y chromosome has experienced significant pseudogenization (Papadopulos et al. 2015) and ∼85% gene deletion (Bergero et al. 2015; Papadopulos et al. 2015). Our results are in line with these recent conclusions as we found extensive gene deletion in R. rothschildianus, with the majority of sex-linked hemizygous genes showing evidence for deletion from the Y. These results suggest that haploid selection is not sufficient to completely halt gene loss on plant Y chromosomes. Nevertheless, haploid gametophytic selection may still contribute to the evolution of plant genomes in general and of their sex chromosomes in particular (Chibalina and Filatov 2011; Arunkumar et al. 2013; Crowson et al. 2017; Scott and Otto 2017).

On the other hand, our results for R. hastatulus, particularly the XY karyotype, contrast with those from S. latifolia, despite apparently comparable ages of sex chromosomes. This discrepancy may be explained by differences in the intensity of selective interference in the two systems. Factors such as the strength of selection, gene density on the sex chromosomes, size of the nonrecombining region, and time since the loss of recombination are important factors for consideration (Kaiser and Charlesworth 2009; Crowson et al. 2017). Although precise values for most of these metrics are difficult to ascertain, the difference in the distribution of changes at synonymous sites (ds) between R. hastatulus (Hough et al. 2014) and S. latifolia (Bergero et al. 2007; Papadopulos et al. 2015) suggests varying patterns and histories of the loss of recombination, and could explain the contrasting amount of gene loss observed between these two species. However, our evidence for a greater proportion of deleted relative to silenced hemizygous genes in the XYY karyotype of R. hastatulus compared with the XY karyotype suggests that the extent of selective interference may well be an important determinant of the amount and rate of gene loss.

The use of a transcriptome assembly may have limited the number and types of genes that we found and bias our study to genes expressed in leaf tissue. In particular, if sexually (Rice 1987) or ploidally (Scott and Otto 2017) antagonistic genes are enriched on the sex chromosomes but are depleted in leaf tissue (where we expect relatively little sexual antagonism), we may be limiting the breadth of genes identified. However, tissue specific expression seems to be rare in plants; estimates suggest that approximately ¾ of genes expressed in one tissue are expressed in all tissues (Schmid et al. 2005; Zhang et al. 2017). Furthermore, we found that 70–80% of our genomic reads mapped to our assembly (supplementary table S1, Supplementary Material online). Thus, we expect our conclusions to be broadly robust, although the patterns of gene deletion and silencing may be different for tissue-specific genes.

Y Chromosome Degeneration Is a Progressive Process

Despite extensive study of sex chromosome systems, the relative contributions of adaptive silencing, Hill–Robertson interference, and neutral genetic drift to Y chromosome degeneration have yet to be fully disentangled. Because genes on nonrecombining Y chromosomes should be under less efficient selection than their X counterparts, it may be beneficial for Y chromosome copies to be adaptively silenced (Peck 1994; Orr and Kim 1998). This pattern is expected to be strongest for genes that have experienced positive selection on the X (Crowson et al. 2017). Previous work in R. rothschildianus by Crowson et al. (2017) found no evidence of elevated rates of positive selection for hemizygous genes. However, because the test used in that study relied on recurrent positive selection over long evolutionary timescales, we investigated whether signs of recent within-species selective sweeps could remain.

In agreement with Crowson et al. (2017), our population genetic analyses found no evidence for greater positive selection on X-homologs of Y-silenced genes. Thus, we find no evidence for adaptive silencing of Y chromosome genes due to differential positive selection. However, our results alone cannot exclude the possibility that degenerate genes on the Y are actively silenced. In particular, Y copies may also be adaptively silenced if they have accumulated many deleterious mutations and their products interfere with proper cell function (Charlesworth 1978). As Y-linked copies accumulate deleterious mutations, selection may act to silence them even in the absence of positive selection on the X. This type of adaptive silencing is especially likely for highly constrained or highly expressed genes. But Crowson et al. (2017) found that genes that are less constrained and that have lower expression are more likely to be lost from the Y. This finding suggests less functionally important genes or dosage insensitive genes are likely to be lost first, in contradiction to the prediction that important genes would be silenced as they degenerate. Therefore, there appears to be little evidence in support of adaptive silencing of Y-linked genes in Rumex.

Two lines of evidence suggest that degeneration may be a progressive process, where loss of expression may often occur prior to gene loss. First, we find a much higher proportion of gene loss relative to silencing in hemizygous genes of R. rothchildianus than in R. hastatulus. Given that the former system shows a greater proportion of hemizygous genes and is thought to have had a much larger nonrecombining region for a longer period, R. hastatulus may be at an earlier stage of degeneration and shows a concomitantly higher proportion of silenced genes. Second, genes that have recently become hemizygous within one karyotype of R. hastatulus are mostly silenced, and show little evidence for gene loss. These patterns are consistent with the hypothesis that gene loss may often follow loss of expression.

Our results suggesting that gene silencing precedes deletion imply that gene loss may often proceed neutrally following reduction of gene expression on the Y. This result is in line with findings on the very young neo-Y chromosomes of Drosophila albomicans (Zhou and Bachtrog 2012a). Zhou et al. (2012a, 2012b) concluded that early heterochromatin accumulation following transposable element (TE) invasion on the sex chromosomes is likely to cause the observed degeneration in Drosophila. TE accumulation is likely in regions of low recombination (reviewed in Kent et al. 2017) and this accumulation followed by silencing, which can spread to adjacent genes, probably contributes to Y-chromosome degeneration. Further genomic work in Rumex examining the possible role of TE invasion in heterochromatin formation and gene expression loss on the Y chromosome would be profitable.

Our results suggest that a significant amount of gene degeneration occurs neutrally, however identifying the relative roles of degeneration by neutral drift following early silencing, from that of the accumulation of deleterious mutations by selective interference, remains unclear. Recent findings in Rumex suggest interference selection is likely to still play a major role in Y chromosome evolution. In particular, a recent molecular population genetic analysis of R. hastatulus reported Y-linked polymorphism to be dramatically reduced compared with X and autosomal polymorphism (Hough et al. 2017), in line with expectations of interference selection on Y chromosomes (see Charlesworth 2012). This analysis also suggested that the observed elevated rate of molecular evolution of genes remaining on the Y chromosome was explainable by interference selection, without needing to invoke additional neutral processes caused by gene silencing. Therefore, although our evidence suggests that the final stages of gene loss may be neutral following silencing, this does not preclude a major role for Hill–Roberson interference during gene degeneration.

Y Silencing by Selective Interference

One way to unite the role of linked selection and gene silencing is a model in which change in expression is itself caused by Hill–Robertson interference. This model is not mutually exclusive with the model of early TE invasion and epigenetic silencing, and both are likely to contribute to a varying extent in different systems. Hill–Robertson interference lowers effective population size (Ne) in the low recombination region (reviewed in Comeron et al. 2008; Good et al. 2014). This reduced Ne broadens the range of polymorphism that fall into the “nearly neutral zone,” where the effects of drift overwhelm those of selection (Ohta and Tachida, 1990). This effect has been well-described for coding regions and is thought to explain elevated rates of amino acid evolution on the Y chromosome as compared with the X (Charlesworth 1996; Hough et al. 2017). However, the consequences of Hill–Robertson interference can also be extended to regulatory regions such as promoters and enhancers. The successive fixation of mutations that each slightly decrease expression would cause lower expression of genes on the Y and could culminate in complete gene silencing. The lower expression of Y alleles than X alleles has been observed in Rumex, Silene, and many animal systems (Bergero and Charlesworth 2011; Chibalina and Filatov 2011; Meisel et al. 2013; Vicoso et al. 2013; Hough et al. 2014; Singh et al. 2014; Papadopulos et al. 2015). Such mutations causing slight decreases in gene expression may be more likely to fix than large-effect indels so that gene silencing would appear more likely than deletion. Following silencing, genes would be free to degenerate neutrally.

This expression-drift model is like that proposed for the process of loss of gene duplicates following whole genome duplication, whereby mildly deleterious mutations affecting gene expression lead to eventual loss of expression of one of two duplicate genes (Gout and Lynch 2015; Thompson et al. 2016). Although the dosage compensatory drift models assume a symmetrical distribution of mutations affecting expression, the general loss of expression on the Y implies that new mutations in regulatory regions generally decrease expression, or mutations that increase expression are often more deleterious. Some studies do in fact suggest that new mutations are more likely to decrease than increase expression, including empirical findings from de novo mutations and polymorphism in enhancers in mice (Kwasnieski et al. 2012; Patwardhan et al. 2012), Drosophila (Arnold et al. 2014), and yeast (Metzger et al. 2015). Therefore, mutation pressure and inefficient selection may drive the loss of expression of Y-linked genes until they become fully silenced.

Several genomic studies have reported that regulatory regions are under weaker purifying selection than coding regions, making them more likely than protein-coding changes to be susceptible to mutational decay on the Y chromosome. In Caenorhabditis elegans, Hill–Robertson interference is a predominant force acting on expression variation (Rockman et al. 2010). In Capsella grandiflora, most expression QTLs (eQTLs) experience weak purifying selection (Josephs et al. 2015; Steige et al. 2017), and the distribution of fitness effects is lower in regulatory regions of the Tobacco etch virus than in regions of nonsynonymous change (Bernet and Elena, 2015). Furthermore, conserved noncoding sites (CNS) are thought to be enriched for regulatory regions (Josephs et al. 2015) and, in Capsella, CNS experience intermediate or weak purifying selection (Williamson et al. 2014). Similarly, in mice 35% of CNS mutations, as compared with 77% of coding sites, have a fitness effect more than NeS = 10, and 44% of CNS mutations have fitness effect less than NeS =1 (Halligan et al. 2013). Work in Drosophila suggests selection on the order of Ns = 1 for deviations from optimal expression of 1 SD (Bedford and Hartl 2009). Overall, we therefore expect a large class of nearly neutral mutations to lie within regulatory regions, leaving open the possibility of significant mutation accumulation on Y chromosomes. Future work investigating the rate of mutational decay in regulatory sequences compared with protein-coding regions will help assess the possible role of regulatory decay in the earliest stages of gene loss on the Y chromosome.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.