Genomic Loss and Silencing on the Y Chromosomes of Rumex

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 (log₂=−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.

Polymorphism data of females from XY R. hastatulus populations used for π and Tajima’s D estimates are from Hough et al. (2017), 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).

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 log₂ 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 log₂ 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 χ² contingency test indicated significantly more silencing relative to deletion in the XY system as compared with the XYY system (χ²= 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, χ² =24.63, P =6.96e-07).

Fig. 1.

—Distribution of log₂ 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. 2a). This proportion increases to 56% with no genomic heterozygosity (putatively deleted) in the XYY karyotype (fig. 2b). These proportions are not significantly different from the proportions predicted using coverage ratio (for the XY karyotype: χ² = 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.

Fig. 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 log₂ 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 (χ² =4.2, P =0.04) (fig. 3 and table 1).

Fig. 3.

—Distribution of log₂ 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).

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 (N_e) in the low recombination region (reviewed in Comeron et al. 2008; Good et al. 2014). This reduced N_e 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 N_eS = 10, and 44% of CNS mutations have fitness effect less than N_eS =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.