Evolution of Young Sex Chromosomes in Two Dioecious Sister Plant Species with Distinct Sex Determination Systems.

In the last decade, progress has been made in methods to identify the sex determination system in plants. This gives the opportunity to study sex chromosomes that arose independently at different phylogenetic scales, and thus allows the discovery and the understanding of early stages of sex chromosome evolution. In the genus Silene, sex chromosomes have evolved independently in at least two clades from a nondioecious ancestor, the Melandrium and Otites sections. In the latter, sex chromosomes could be younger than in the section Melandrium, based on phylogenetic studies and as no heteromorphic sex chromosomes have been detected. This section might also exhibit lability in sex determination, because male heterogamy and female heterogamy have been suggested to occur.In this study, we investigated the sex determination system of two dioecious species in the section Otites (Silene otites and its close relative Silene pseudotites). Applying the new probabilistic method SEX-DETector on RNA-seq data from cross-controlled progenies, we inferred their most likely sex determination system and a list of putative autosomal and sex-linked contigs. We showed that the two phylogenetically close species differed in their sex determination system (XY versus ZW) with sex chromosomes that derived from two different pairs of autosomes. We built a genetic map of the sex chromosomes and showed that both pairs exhibited a large region with lack of recombination. However, the sex-limited chromosomes exhibited no strong degeneration. Finally, using the "ancestral" autosomal expression of sex-linked orthologs of nondioecious S. nutans, we found a slight signature of dosage compensation in the heterogametic females of S. otites.

Introduction

In numerous species with distinct male and female individuals, that is, gonochorism in animals and dioecy in plants, the sexual phenotype is determined by sex chromosomes (Bachtrog et al. 2014). Sex chromosomes evolve from a pair of autosomes that underwent a recombination suppression between sex-determining loci (Bergero and Charlesworth 2009). Consequently, in diploids, recombination is suppressed only in heterogametic individuals producing two types of gametes as opposed to homogametic individuals, producing only one type of gamete and carrying fully recombining sex chromosomes. The Y chromosome in a male heterogametic systems and the W chromosome in female heterogametic systems are therefore sex-limited chromosomes that accumulate deleterious mutations because of the absence of recombination that leads to a reduction of the effective size and an increase of Hill–Robertson effects (Hough et al. 2017) affecting the selection efficacy to purge deleterious mutations (Charlesworth and Charlesworth 2000; Bachtrog, 2006). Gene losses and deficit in gene expression are therefore expected due to indels in regulatory sequences or changes in methylation patterns and constitute one major feature of genomic degeneration. Such genomic changes are mainly deleterious, and mechanisms that counteract these losses are thus expected to be selected for. Notably, dosage compensation has been described in many and diverse organisms, that is, modified expression of alleles on X (or Z) chromosomes or autosomes that has been selected to reach the level of the homogametic sex or to resolve the imbalance between sex chromosomes and autosomal gene expression (Disteche 2016; Gu and Walters 2017). However, not all sex chromosomes evolved such mechanism and its ubiquity has recently been toned down (reviewed in Mank [2013] and Muyle et al. [2017]).

Beyond this theoretical sketch, the study of sex chromosomes on a large survey of the tree of life has revealed that evolutionary pathways can be quite diverse and genera, species, and even population specific (Ogata et al. 2008; Bachtrog et al. 2014). In particular, some sex chromosomes lacking recombination are found homomorphic in lineages that have a long evolutionary history of gonochorism/dioecy, suggesting that mechanisms can prevent degeneration, by resetting the mutation load of the sex-limited chromosome through recurrent turnover of sex chromosomes or occasional recombination in sex-reverted individuals (Perrin 2009). The study of homomorphic sex chromosomes, either young or old, is crucial to assess the tempo and understand the evolutionary pathways leading to differentiated sex chromosomes.

In the last decade, progress has been made to identify the sex determination system in plants: Indeed, sex chromosomes were reported in 48 species across 20 families and at different evolutionary stages (Ming et al. 2011; Harkess and Leebens-Mack 2017; Muyle et al. 2017). This gives the opportunity to study sex chromosomes that arose independently at different phylogenetic scale, and thus allows the discovery and the understanding of early stages of sex chromosome evolution, which is not the case for most model cases in animals. In the genus Silene, sex chromosomes have evolved independently in at least two clades from a nondioecious ancestor, the Melandrium section and Otites section (Desfeux et al. 1996; Mrackova et al. 2008; Marais et al. 2011; Slancarova et al. 2013). The Melandrium section includes the well-studied species Silene latifolia with heteromorphic XY chromosomes (Bernasconi et al. 2009), and Silene diclinis with its multiple sex chromosomes X, Y1, and Y2, resulting from a reciprocal translocation between the ancestral Y and an autosome (Howell et al. 2009). The sex chromosomes of S. latifolia are estimated to be 5–11 Myr old (Rautenberg et al. 2010; Marais et al. 2011; Krasovec et al. 2018), and their origin is complex, involving multiple translocation events (Qiu et al. 2016). Notably, dosage compensation has been reported despite their recent evolution (Muyle et al. 2012, 2017; Papadopulos et al. 2015).

In the section Otites, only homomorphic sex chromosomes were detected while composed by numerous dioecious species (Oxelman et al. 2013) including S. otites, S. densiflora (Favarger 1946), and S. colpophylla (Mrackova et al. 2008). This suggests that sex chromosomes of the Otites section are younger than those from S. latifolia, assuming the absence of mechanisms resetting mutation load in the sex chromosomes considered. More interestingly, this section exhibits different sexual systems, with a male heterogametic system in S. colpophylla (Mrackova et al. 2008) and S. pseudotites (Sansome 1938), and a female heterogametic system in S. otites (Sansome 1938; Slancarova et al. 2013). However, male heterogamety has also been suggested to be the sexual system of S. otites (Warmke 1942; reviewed by Westergaard 1958) and the taxonomic status of S. pseudotites is ambiguous as it has also been described as a subspecies of S. otites in the Atlas Florae Europaeae (Tutin et al. 1964) or a hybrid between S. otites and S. colpophylla by Wrigley (1986). Therefore, clarification is needed to assess the sex determination system of the two sister species, as they can constitute a new model to decipher the evolutionary mechanisms involved in sex chromosome turnover.

Applying the novel probabilistic method SEX-DETector (Muyle et al. 2016) on RNA-seq data from cross-controlled progenies of S. otites and S. pseudotites, we inferred their most likely sex determination (either XY, ZW, or absence of sex chromosomes) and a list of putative autosomal and sex-linked contigs. We genetically mapped some of these sex-linked contigs and determined their autosomal origin, using the published genetic map of S. latifolia (Papadopulos et al. 2015). We showed that the two sister species differed in their sex determination system, S. pseudotites carrying XY sex chromosomes and S. otites ZW sex chromosomes, each derived from two distinct pairs of ancestral autosomal chromosomes. Both sex chromosomes exhibited a nonrecombining region but exhibited low degeneration, in accordance with their young age. However, we did observe a signature of dosage compensation in heterogametic females of S. otites, when referred to the “ancestral” autosomal expression of nondioecious S. nutans.

Materials and Methods

Plant Material and Sequencing

To identify sex-linked genes, we generated RNA-seq data from a controlled cross for each species. Parents came from seeds sampled in natural populations in France. For S. otites, parents came from a seed bulk from a population located at Valpuisseaux in the Parisian region (latitude N 48.394448 and longitude E 2.301264) collected in 2002 by the Conservatoire Botanique National du Bassin Parisien (provided by Chantal Griveau). For S. pseudotites, parents came from two different out-pollinated maternal progenies collected in 2012 by Jos Kafer (LBBE, Lyon) in a population located in Caussols in Alpes-Maritimes (latitude N 43.743111 and longitude E 6.938917). Vouchers of S. otites and S. pseudotites populations are stored in the herbarium of Conservatoire Botanique National de Bailleul (www.cbnbl.org; last accessed January 23, 2019): BAIL 2019-001 for S. otites and BAIL 2019-002 for S. pseudotites.

The parents and the progenies were sown and grown in the greenhouse with standard conditions (20 °C, 16-h day length). We extracted total RNA from young flower buds of parents, five males and five females of each progeny using the Spectrum Plant Total RNA kit (Sigma. Inc., USA) following the manufacturer’s protocol and treated with a DNAse. Libraries were prepared with the TruSeq RNA Strand Specific Preparation kit (Illumina Inc., USA). Each cDNA library was sequenced using a paired-end protocol on a HiSeq2500 sequencer, producing 100-bp reads (six libraries pooled in equiproportion per lane). Demultiplexing was performed using CASAVA 1.8.1 (Illumina) to produce paired sequence files containing reads for each sample in Illumina FASTQ format. RNA extraction, library preparation, and sequencing were done by the sequencing platform in the AGAP laboratory, Montpellier, France (http://umr-agap.cirad.fr/).

To estimate divergence of the sequences, we used RNA-seq data from flower buds of one hermaphrodite individual of S. nutans, considered as an outgroup species. To test for dosage compensation, we added RNA-seq data from flower buds from five additional hermaphrodite individuals of S. nutans sampled in natural populations and grown in the greenhouse in the same standard conditions, reaching a total of six hermaphrodites.

RNA-seq Cleaning, Assembly, and Genotyping

Raw data of the two families were filtered for sequencing adapters using Cutadapt (Martin 2011), low read quality and poly-A tails with PRINSEQ (Schmieder and Edwards 2011). Reference trancriptomes were assembled de novo mixing all members of a family using default settings of Trinity (Haas et al. 2013) and an additional Cap3 run with default parameters (Huang and Madan 1999) (supplementary table 1, Supplementary Material online). Finally, coding sequences were predicted using TransDecoder (see website: http://transdecoder.github.io/; last accessed January 23, 2019, Haas et al. 2013). Reads from each member of a family were mapped onto their reference transcriptome using Bowtie2 (Langmead and Salzberg 2012). SAM files were compressed and sorted using SAMTOOLS (Li et al. 2009). At each position in individual alignments, diploid genotypes were called according to the method described by Tsagkogeorga et al. (2012) and improved by Gayral et al. (2013), implemented in the reads2snps program. We required a minimum coverage of 10× per position and per individual to call a genotype following Burgarella et al. (2015). Data were not cleaned for paralogous single nucleotide polymorphisms (SNPs). Indeed, XY or ZW SNPs can lead to the identification of paralogous contigs because recombination suppression between X and Y, or Z and W decrease allelic similarity, resulting in highly diverged allele detected as paralogs. Therefore, XY or ZW SNPs tend to be filtered out by the paraclean option implemented in reads2snp that uses a likelihood ratio test based on explicit modeling of paralogy.

Identification of Sex-Linked Genes and Autosomal Origin

The probabilistic methods implemented in SEX-DETector (Muyle et al. 2016) allow to statistically test the presence of sex chromosomes in the data (either an XY or a ZW system). This relies on the Bayesian information criterion (BIC). The lower the BIC, the better the model fits the data. The BIC takes into account the likelihood of the model as well as the number of parameters of the model in order to compare model fit to the data. Contig segregation types were inferred after estimation of parameters using a likelihood-based inference method called EM algorithm (for Expectation Maximization). Posterior segregation type probabilities were filtered to be higher than 0.9. Briefly, the method is based on the genotypes of two parents and ten of their progeny (five males and five females) from which we infer the segregation type of each contig. The model considers that SNPs can be transmitted to the progeny by three segregation modes: 1) autosomal, 2) sex-linked with both X and Y (or Z and W) alleles present, and 3) X (or Z) hemizygous, that is, sex linked with only the X (or Z) allele present (the Y or W allele being inactivated, lost, too weakly expressed or in a different contig due to X/Y or Z/W divergence). The method uses a likelihood-based framework to assess the “posterior” probability of each segregation type for each SNP, given the observed genotype data. The segregation type of each contig is finally inferred by averaging SNP posteriors. Contigs are classified as undefined if no suitable polymorphic sites are available.

The genetic map and reference sequences of S.latifolia from Papadopulos et al. (2015) were used to map contigs of both species using reciprocal best hits on BlastN results between Open Reading Frames of both species and scaffolds of S. latifolia. A hit was considered as valid when alignment length was above 130 bp, sequence similarity above 80%, and e value below e−50.

Genetic Mapping of Putative Sex Chromosomes in S. pseudotites and S. otites

We developed SNP markers for contigs that were inferred as sex linked in S. pseudotites or S. otites (supplementary table 2, Supplementary Material online) with additional autosomal contigs that located on the same S. latifolia linkage group (LG): 10 SNPs on LG6 for S. pseudotites (2 autosomal and 8 XY linked), and 16 SNPs on LG1 (2 autosomal and 14 XY-linked) and 19 ZW-linked SNPs on LG3 for S. otites. The SNPs were chosen in order to be uniformly distributed along the LGs, based on S. latifolia genetic map (Papadopulos et al. 2015).

The SNPs were used to genotype the progenies from the same crosses as the ones used in the SEX-DETector analysis: 37 individuals for S. pseudotites (12 males, 22 females, and 3 with unknown phenotype) and 105 individuals for S. otites (61 males and 44 females).

The Kompetitive Allele-Specific Polymerase chain reaction (KASP) genotyping assay was used to detect the SNPs. This method is based on competitive allele-specific polymerase chain reaction and enables biallelic scoring of single nucleotide polymorphisms. The SNP-specific KASP Assay mix (designed by LGC group, www.lgcgroup.com; last accessed January 23, 2019) and the universal KASP Master mix are added to DNA samples, a thermal cycling reaction is then performed, followed by an end-point fluorescent read. Biallelic discrimination is achieved through the competitive binding of two allele-specific forward primers, each with a unique tail sequence that corresponds with two universal fluorescence resonant energy transfer cassettes; one labeled with FAM dye and the other with HEX dye. The sequences flanking the SNPs that were used to define the primers are given in supplementary table 2, Supplementary Material online.

Genetic maps of the putative sex chromosomes were built using Joinmap3 (Van Ooijen and Voorrips 2001) with a minimum logarithm of the odds-grouping score of 2 and using the Kosambi’s mapping function.

Divergence of Sex-Linked Contigs

We used the putative X-, Y-, Z-, and W-linked sequences inferred by SEX-DETector to test whether the Y- or W-linked copies had incorporated more changes than the X- or Z-linked alleles, as expected if the Y or W is degenerated. We used S. nutans sequences as outgroups. Orthologs of focal species and S. nutans sequences were identified from reciprocal best hits on BlastN results and aligned using MACSE (Ranwez et al. 2011). The number of synonymous and nonsynonymous substitution and the synonymous (dS) and nonsynonymous (dN) substitution per site were calculated using the PopPhyl tool dNdSpiNpiS (Gayral et al. 2013). Although only one sequence of Y (from the father) or W (from the mother) can be reconstructed, up to three alleles can be identified for X- or Z-linked contigs. To avoid bias due to the number of sequences used to calculate these statistics, we randomly chose only one sequence of the putative X or Z. Finally, synonymous divergence max (dSXYmax or dSZWmax) between the X and Y (or Z and W) haplotype contigs were estimated using dNdSpiNpiS, corresponding to the dSXY or dSZW mean of the 5% highest values.

Dosage Balance and Compensation

Dosage compensation in S. pseudotites and S. otites was studied comparing expression level of sex-unbiased sex-linked contigs of the focal species with the expression of the orthologs in S. nutans. If dosage compensation occurs, the expression in the heterogametic sex should be equivalent to the expression in the homogametic one, whereas the allele expression in the heterogametic sex should be different, the Y or the W being lower than the X or the Z. To be able to study separately the X, Y, Z, and W allele expression, expression was studied at the SNP level.

To identify sex-unbiased contigs, we summed the number of reads at each SNP within a contig, sum that was then divided by the number of SNPs within a contig to obtain an average expression over SNPs for each contig. The average expression was then normalized using the TMM normalization method implemented in edgeR (Robinson et al. 2010). A QL F-test was done to test for differential expression between sexes, using edgeR. A contig is considered as biased if its expression is significantly at least doubled in one sex compared with the other one (fold changes significantly >2).

We focused on the average expression over SNPs of contigs of the six females of a family (the mother and its five daughters), the six males of a family (the father and its five sons), and six hermaphroditic individuals of S. nutans. For each contig, we computed the following log 2 ratios: (Xm + Ym):XXf, which should be close to 0 as we worked on sex-unbiased contigs; Xm:XXf, which should be greater than −1 if dosage compensation occurs; and Ym:Xm as an estimate of Y degeneration. To have a more accurate estimation of Y degeneration and possible subsequent dosage compensation in males, we also used S. nutans expression of sex-linked orthologs as a proxy of ancestral autosomal expression (Mank 2013), and thus computed the following log 2 ratios: XXf:AAnut, (Xm + Ym):AAnut, Xm:AAnut, and Ym:AAnut. In the same way, for Z/W contigs, we computed the following log 2 ratios: (Zf + Wf):ZZm, Zf:ZZm, Wf:ZZm, Wf:Zf, ZZm:AAnut, (Zf + Wf):AAnut, Zf:AAnut, and Wf:AAnut.

Results

Sex Determination System of S. pseudotites and S. otites

The analysis by SEX-DETector revealed that, for both species, the model without sex chromosomes was the less likely (with the highest BIC, see Materials and Methods for details), suggesting the presence of sex chromosomes (table 1).

Table 1

Model Comparison Using SEX-DETector on Silene pseudotites and Silene otites

Species	Sex Determination System	Chromosomal Category	Number of Contigs	BIC
Silene pseudotites	XY	Sex linked	174	3,422,710.589
		Autosomal	7,233
		Undefined	37,351
	ZW	Sex linked	0	3,430,155.856
		Autosomal	7,608
		Undefined	37,150
	No sex chromosome	Sex linked	0	3,438,114.512
		Autosomal	8,638
		Undefined	36,120
Silene otites	XY	Sex linked	222	2,844,016.544
		Autosomal	4,827
		Undefined	50,005
	ZW	Sex linked	329	2,844,920.705
		Autosomal	5,232
		Undefined	49,493
	No sex chromosome	Sex linked	0	2,862,178.838
		Autosomal	7,008
		Undefined	48,046

For S. pseudotites, the XY system was the most likely with the lowest BIC and 174 contigs were inferred as sex linked, whereas the ZW system had a higher BIC and no contigs inferred as sex linked (table 1). Hereafter, we therefore considered that S. pseudotites had an XY system. We did not detect any X-hemizygous loci, that is, only found on the X chromosome and not detected on the Y, due to either their divergence or loss on the Y. Over the 174 contigs inferred as sex linked, 167 contigs had an ortholog on S. latifolia scaffolds from Papadopulos et al. (2015) and 44 of them were mapped on their genetic map linked (supplementary table 2, Supplementary Material online). In the same way, over 7,233 contigs inferred as autosomes in S. pseudotites, 6,668 contigs had an ortholog on S. latifolia scaffolds and 1,154 of them were mapped. Although inferred autosomal contigs of S. pseudotites were spread over all LGs of S. latifolia (including the X chromosome of S. latifolia), inferred sex-linked contigs of S. pseudotites were grouped together on the LG6 of S. latifolia (fig. 1).

. 1.

—Genetic map of (a, b) Silene pseudotites and (c, d) Silene otites. (a, c) Contigs with chromosomal category inferred by SEX-DETector were mapped on the genetic map of Silene latifolia (Papadopulos et al. 2015): contigs inferred as autosomal are in gray, contigs inferred as sex linked and among which some were confirmed using a genetic map are in orange and in S. otites, contigs inferred as XY linked but subsequently found to be associated with a male-specific lethal locus (Msl1) are in blue. (b, d) Comparative genetic mapping of LG6 from S. pseudotites and LG1 and LG3 from S. otites. Chromosomal regions containing loci in complete linkage with sex are indicated by an orange star. Loci mapped on LG1 are represented by a circle and the colors indicate distorted segregation in S. otites males: The gradient color goes from white (40% of males are homozygous, i.e., close to 1:1 segregation as expected), to blue (up to 80% of males are heterozygous, i.e., there is segregation distorsion).

The sex determination system of S.otites was less clear at first sight. Although the XY model had a slightly lower BIC than the ZW model, it also inferred a hundred sex-linked contigs less than the ZW model (table 1). We hereafter considered as autosomal contigs that were inferred as autosomal under both XY and ZW models (i.e., 2,160 contigs), XY sex-linked contigs as contigs that were inferred as sex linked under the XY model and undefined under the ZW model (i.e., 155 contigs) and ZW sex-linked contigs when they were inferred as sex linked under the ZW model and undefined under the XY model (i.e., 188 contigs). Once again, as it was the case for S. pseudotites, no hemizygous sex-linked loci were detected. Over the 2,160 inferred autosomal contigs, 2,008 contigs had an ortholog on a S. latifolia scaffolds and 326 of them were mapped on the genetic map of S. latifolia. In the same way, over the 155 inferred XY sex-linked contigs, 143 contigs found an ortholog on a S. latifolia scaffolds and 37 of them were mapped. And over the 188 inferred ZW sex-linked contigs, 75 contigs found an ortholog on S. latifolia scaffolds and 44 of them were mapped (see supplementary table 2, Supplementary Material online, for sex-linked contig description). Although inferred autosomal contigs of S. otites were spread over all LGs of S. latifolia (including the X chromosome of S. latifolia), inferred XY sex-linked contigs were grouped on LG1 of S. latifolia and inferred ZW sex-linked contigs were grouped on LG3 (fig. 1).

Genetic Linkage Map of Sex-Linked Contigs of S. pseudotites and S. otites

To validate the sex-linked status of some contigs and investigate the level of recombination among them, we constructed a genetic linkage map for each species using segregating crosses for sex. Using the same cross as above, we genotyped 37 individuals for S. pseudotites (12 males, 22 females, and 3 with unknown phenotype) and 105 individuals for S. otites (61 males and 44 females).

In S. pseudotites, the 8 XY sex-linked and 2 autosomal SNPs, as inferred by SEX-DETector, formed a single LG with a total map length of 32 cM, compared with 61 cM in S. latifolia. This reduction of genetic length was due to the complete linkage of the 8 XY-linked SNPs (indicating recombination cessation), whereas the most distant loci were 29 cM apart in S. latifolia LG6. They were also completely linked to the sexual phenotype confirming their status of sex-linked loci: all females were homozygous and all males were heterozygous at these markers, as expected in a XY system (fig. 1).

We then tried to disentangle the seemingly complex sex determination system of S. otites with potentially two pairs of sex chromosomes. The 19 SEX-DETector inferred ZW-linked SNPs exhibited a very strong linkage (14 in complete linkage) with the sex determination region (females being heterozygous and male homozygous), generating a linkage map with a reduced length when compared with S. latifolia LG3: the LG of S. otites had an estimated length of 4 cM, whereas the corresponding LG in S. latifolia was 57 cM long (fig. 1 and supplementary table 2, Supplementary Material online).

Using 2 SEX-DETector inferred autosomal and 14 SEX-DETector inferred XY-linked SNPs (supplementary table 2, Supplementary Material online), we were able to build a map that had a length of 52 cM in comparison with the 69 cM length of S. latifolia LG1 (fig. 1). Although XY-linked SNP segregation was independent with sex (being either homozygous or heterozygous in more or less equal proportion), some loci exhibited a strong segregation distortion in males, male being preferentially heterozygous at these loci (supplementary table 2, Supplementary Material online). A putative gene acting epistatically in males must be located in the vicinity of loci exhibiting distorted segregation in males. We observed an increase of the ratio of heterozygous males/total of males in a window around 50 cM in our LG1 genetic map (fig. 1).

Overall, we confirmed the existence of two different sex determination systems, XY in S. pseudotites and ZW in S. otites, with sex chromosomes that derived from two different autosomes. The additional XY system inferred by SEX-DETector and corresponding to the autosomal LG1 in S. latifolia seems to rather be the signature of a locus acting negatively in males when homozygous for a recessive allele. Therefore, we did not consider this LG to be a genuine sex chromosome and thus only assessed the XY system of S. pseudotites and ZW system of S. otites in the rest of the study.

Divergence of Sex Chromosomes

Sex-limited chromosomes (Y and W) should accumulate more deleterious mutations than X and Z chromosomes due to their lower effective size. We calculated the number of synonymous and nonsynonymous substitutions and the synonymous (dS) and nonsynonymous (dN) divergence using S. nutans as outgroup. In S. pseudotites, both the number of synonymous and the number of nonsynonymous substitutions were slightly but significantly higher in the Y haplotype than in the X haplotype (G-test, G = 4.491, P value = 0.034, and G = 9.692, P value = 0.002, respectively, table 2). Therefore, dS and dN were slightly higher in the Y haplotype (table 2). On the contrary, in S. otites, we only found a small significant excess of nonsynonymous substitution in the W-linked contigs (G-test, G = 3.952, P value = 0.046), otherwise differences were not significant (table 2). The divergence between the X and Y haplotypes of S. pseudotites (average dSXY = 0.0192 and dSXYmax = 0.0561) was higher than the divergence between the Z and W haplotypes of S. otites (average dSZW = 0.0111 and dSZWmax = 0.0352).

Table 2

X-Y and Z-W Divergence

	Silene pseudotites		Silene otites
	X Linked	Y Linked	Z Linked	W Linked
Total number of contigs	125		127
Number of synonymous substitution	1,821	1,951	1,798	1,860
dS	0.060	0.0647	0.0525	0.0543
	[0.0548; 0.0667]	[0.0586; 0.0715]	[0.0447; 0.0621]	[0.0462; 0.0640]
Number of nonsynonymous substitution	782	910	980	1,070
dN	0.0087	0.0101	0.0096	0.0105
	[0.0075; 0.0100]	[0.0088; 0.0115]	[0.0077; 0.0122]	[0.0086; 0.0131]
dN/dS	0.1439	0.1563	0.1839	0.1941
	[0.1242; 0.1660]	[0.1355; 0.1798]	[0.1558; 0.2124]	[0.1664; 0.2219]
dSXY or dSZW	0.0192		0.0111
	[0.0166; 0.0222]		[0.0093; 0.0131]

Note.—Synonymous dS and nonsynonymous dN substitution per site using S. nutans as outgroup and the synonymous divergence between the X (or Z) and Y (or W) haplotype (dSXY or dSZW).

Degeneration and Dosage Compensation of the Sex-Limited Chromosome

The expression of sex-linked contigs in both species was highly correlated between males and females (supplementary fig. 1, Supplementary Material online). In S. pseudotites out of 174 XY-linked contigs, 151 were sex unbiased (after removing five contigs with missing data), whereas in S. otites, 160 were sex unbiased out of 188 ZW-linked contigs (after removing 10 contigs with missing data). The fact that most sex-linked genes were sex unbiased could suggest that the Y or W chromosomes have a limited degeneration or that dosage compensation is effective on the X or Z chromosomes in the heterogametic sex.

Focusing on sex-unbiased expressed genes that were heterozygous in the heterogametic sex, we first compared the levels of expression of the sex chromosomes in the heterogametic sex with the homogametic sex: we expected the Y or W chromosome to be less expressed than X or Z in the homogametic sex in case of degeneration, and a higher expression of X or Z in the heterogametic sex than in the homogametic sex in case of dosage compensation.

In S. pseudotites, (Xm + Ym):XXf ratio was centered to zero as expected for sex-unbiased-expressed contigs. The Xm:XXf median was significantly higher than the Ym:XXf median (Wilcoxon signed-rank paired test, V = 6,500, P value < 10−5). Consequently, Ym:Xm had a distribution centered slightly below 0 (median = −0.1906, supplementary fig. 2a, Supplementary Material online). Similarly, in S. otites, on contigs sex unbiased in their expression ([Zf + Wf]:ZZm centered on zero), the W allele had a slightly lower expression than the Z allele of the heterogametic sex when compared with the genotypic expression of the homogametic sex (median of Zf:ZZm = −0.8317 and median of Wf:ZZm = −0.9842, Wilcoxon signed-rank paired test, V = 8,228, P value < 10−4, supplementary fig. 2b, Supplementary Material online).

Second, as suggested by Mank (2013), we measured the expression of sex-linked genes in reference to the autosomal and ancestral-like expression of nondioecious S. nutans, a sister species that belongs to the same subgenus.

Overall, expression of S. pseudotites was slightly higher than S. nutans (supplementary fig. 2a, Supplementary Material online). Xm:AAnut was higher than Ym:AAnut (median = −0.8229 and −1.0244, respectively, Wilcoxon signed-rank paired test, V = 3,482, P value < 10−3, supplementary fig. 2a, Supplementary Material online) suggesting a mild Y degeneration, whereas Xm:AA was not significantly higher than Xf:AA (median of Xf:AA = −0.8524, Wilcoxon signed-rank paired test, V = 5,027, P value = 0.572, supplementary fig. 2a, Supplementary Material online). When Xm:AAnut was plotted against Ym:AAnut (a proxy of Y degeneration), there were no outliers to the regression line for low values of Ym:AAnut, and therefore no clear signature of dosage compensation (fig. 2).

. 2.

—Correlation between the expression of X (or Z) and Y (or W) in the heterogametic sex for (a) S. pseudotites and (b) S. otites. Xm = male X allele expression, Ym = male Y allele expression, Zf = female Z allele expression, Wf = female W allele expression, and AAnut = genotype expression of orthologs in the nondioecious S. nutans that do not have sex chromosomes (therefore, orthologs are on autosomes). Diagonal lines across the plot show no change in expression.

In S. otites, the sex-limited chromosome had a reduced expression in the heterogametic sex (medians of Zf:AAnut and Wf:AAnut were −0.8093 and −1.0327, respectively, Wilcoxon signed-rank paired test, V = 4,041, P value = 0.003, supplementary fig. 2b, Supplementary Material online). The Z expression level in the heterogametic sex was significantly higher than the one in the homogametic sex (median of [ZZm/2]:AAnut = −0.9326, V = 1,848, P value < 10−3, supplementary fig. 2b, Supplementary Material online) with a higher level of Zf:AAnut relatively to Wf:AAnut when Wf:AAnut was low (i.e., W chromosome was degenerated, fig. 2). This suggests that some ZW genes are dosage compensated in S. otites.

Discussion

We confirmed that S. pseudotites is male heterogametic as suggested by Sansome (1938). Its sex-linked genes are not orthologous to S. latifolia X-linked genes but to LG6-linked genes (fig. 1). These results support the idea that the sex chromosomes of S. latifolia and S. pseudotites have evolved from different pairs of autosomes, in agreement with an independent origin of sex chromosomes, as suggested by the phylogeny of the genus (Desfeux et al. 1996; Mrackova et al. 2008; Marais et al. 2011). Mrackova et al. (2008) found that S. colpophylla, a close relative to S. pseudotites in the section Otites, is also male heterogametic and did not share sex-linked genes with S. latifolia either. Recent results indicate that S. colpophylla sex chromosomes originated also from autosomal LG6 suggesting that S. colpophylla and S. pseudotites might share the same pair of sex chromosomes (Balounova et al. 2018, https://www.biorxiv.org/content/early/2018/05/28/325068).

Our results suggest that female heterogamety occurs in. S otites also in accordance with the seminal study of Sansome (1938). More recently, Slancarova et al. (2013) found genetic evidence for female heterogamety in S. otites using amplified fragment length polymorphism. Our RNA-seq data from one family confirm the presence of a LG that confers ZW determination. We were able to identify that the ZW sex chromosomes of S. otites are homologous to S. latifolia autosomal LG3.

Both LGs associated with sex in S. pseudotites and S. otites exhibited a reduction of recombination, a classical feature of sex chromosomes (Wright et al. 2016; Kirkpatrick 2017). As observed in another young sex chromosome pair, the loss of recombination might be too recent to impact the genomic structure of the sex-limited chromosome leading ultimately to heteromorphy (Pucholt et al. 2017).

In addition to a ZW system, we detected an additional locus on the chromosome homologous to S. latifolia LG1 that was first inferred as an XY pair by the SEX-DETector analysis of S. otites. Although females could be found at equal frequency at the heterozygous or recessive homozygous state at the locus, males were preferentially found at the heterozygous state. This suggests that one genotype was missing: Male homozygous recessive at the locus that should represent half of the males at the zygote stage.

What are the possible selective forces at play? First, gametic drive in males must be discarded as daughters inherit in equal frequency the two alleles from their father. Therefore, distortion segregation in males must be due to factors acting at the diploid stage of males, most likely at the seed stage because we did not see any drastic juvenile mortality in the progeny. A putative epistatic negative interaction between this male-specific lethal locus on LG1 (Msl1) and one Z-linked locus (or Z-linked loci) when both are homozygous for the recessive allele could cause male lethality, leading to the loss of half of the males. More generally, this negative interaction could involve any factor specifically expressed in males. Such male-specific lethal genes have been described in Drosophila leading to abnormal sex ratio (Fukunaga et al. 1975). These Msl genes can code for proteins involved in the hypertranscription of the single X chromosome in males (Marín and Baker 2000). However, in our case, the male is homogametic.

In addition, it must be noted that we did not observe an unbalanced sex ratio in the studied progeny at the adult stage despite the putative loss of half of the males. Further studies are needed to estimate the occurrence of such an allele in natural populations, and the associated fitness cost in homozygous males.

Molecular Evolution of Young Sex Chromosomes

As a result of selective events occurring in a nonrecombining region, the effective population size of Y-linked genes is expected to be reduced, weakening the efficacy of selection (Charlesworth and Charlesworth 2000). Although we found a higher number of synonymous substitutions in the Y compare to the X haplotype, which can be explained by a higher germline cell-division rate in male than in female, provoking more mutations events on the Y chromosome only present in male (Filatov and Charlesworth 2002; Goetting-Minesky and Makova 2006), dN/dS tended to be higher for Y-linked genes consistent with a reduced efficacy of purifying selection, which could ultimately lead to the degeneration of the Y chromosome (Charlesworth and Charlesworth 2000). Nevertheless, the observed pattern is not as strong as the one observed in S. latifolia (Papadopulos et al. 2015; Krasovec et al. 2018), suggesting a much younger system in S. pseudotites. Indeed, the dSXYmax of S. pseudotites (using the dSXYmax as a proxy of the age of the sex chromosomes when recombination started to stop) is four times as low as S. latifolia’s (0.06 vs. 0.25; Nicolas et al. 2005; Papadopulos et al. 2015). In addition, contrarily to S. latifolia, we did not detect any hemizygous loci in S. pseudotites XY chromosome pair. The absence of hemizygous loci is also observed in S. otites suggesting that its sex chromosomes are younger than S. latifolia ones. The significant lower divergence of Z–W would suggest that S. otites sex chromosomes are even younger than S. pseudotites XY chromosomes. However, the recent phylogenetic study of Balounova et al. (http://dx.doi.org/10.1101/325068) suggests that ZW sex determination system is ancestral. In favor of this hypothesis is the lower number of sex-linked contigs in S. pseudotites than in S. otites, and a less drastic reduction of the recombination rate of the sex chromosome. The lower ZW divergence of S. otites could be explained by a slower evolution than S. pseudotites XY (Ellegren 2011).

Dosage Compensation in Young Systems

When focusing on sex-unbiased sex-linked genes in both species, Y and W alleles were slightly less expressed than X or Z alleles in the heterogametic sex when compared with their autosomal orthologs in S. nutans, used as a proxy of ancestral expression (supplementary fig. 2a, Supplementary Material online). Although a dosage compensation signature was not statistically significant in S. pseudotites, it was significant in S. otites, the increase of expression of Z in females being stronger with the level of degeneration of the expression of the sex-specific chromosome (supplementary fig. 2b, Supplementary Material online and fig. 2). Note that because some of the sex-linked contigs might be located in the pseudoautosomal region, we might underestimate dosage compensation in both species.

Evidence for dosage compensation in young system is rare in plants as well as in animals. However, the 5–10 Myr of S. latifolia sex chromosomes enabled de novo dosage compensation to be selected for some genes as soon as Y expression started to decline (Muyle et al. 2012; Papadopulos et al. 2015). In Rumex hastatulus, sex chromosomes evolved within the past 15–16 Myr and a neo-Y sex chromosome system recently derived (Navajas-Pérez et al. 2005). Both old and young sex-linked genes showed an overall trend of reduced Y expression relative to X-linked alleles while male and female expressions were not different (Hough et al. 2014).

Some female heterogametic species seem to lack global dosage compensation when compared with male heterogametic species (e.g., in chicken, Ellegren et al. 2007; or in snake, Vicoso et al. 2013). Such inequality can be explained by sexual conflict together with stronger selection and greater reproductive variance in males that slow down the selection for dosage compensation in Z compared with X (Mullon et al. 2015).

In S. otites, we found trends of dosage compensation in ZW system that was clearer than in S. pseudotites, while it is expected to be more efficiently selected in a XY system. This pattern is in accordance with the hypothesis that the ZW system is older than the XY system.

Sex chromosome dosage compensation seems not to be universal, and when it occurs in a given species, it might concern only a subset of sex-linked genes (reviewed in Gu and Walters 2017), which are dosage sensitive (e.g., Pessia et al. 2012; Zimmer et al. 2016; Naqvi et al. 2018). Therefore, the difference between S. otites and S. pseudotites in dosage compensation pattern could also be the signature of more dosage-sensitive genes on the S. otites Z than on the S. pseudotites X.

Turnover of Sex Chromosomes in Section Otites

We found evidence of two different sex chromosome pairs between S. pseudotites and S. otites. With S. latifolia and S. diclinis, this rises the number of distinct sex chromosomes to at least 3 in the Silene genus (Howell et al. 2009). In fish (Ross et al. 2009), amphibians (Dufresnes et al. 2015), or diptera (Vicoso and Bachtrog 2015), such diversity due to sex chromosome turnover is well known. In plants, this phenomenon has been described in Salicaceae (Hou et al. 2015) and Fragaria (Goldberg et al. 2010; Wei et al. 2017; Tennessen et al. 2018). Several mechanisms have been suggested to explain the lability of genetic sex determination, among them, sex ratio bias (Ogata et al. 2003; Vuilleumier et al. 2007), X–Y recombination (Dufresnes et al. 2015), or sex-antagonistic selection (van Doorn and Kirkpatrick 2007, 2010) (reviewed in Muyle et al. 2017). Interestingly, on theoretical grounds, turnover is facilitated when the sex-limited chromosome is not too degenerated (explaining that turnover is observed in species with homomorphic sex chromosomes), which is the case of ZW system, considered as the ancestral state. In addition, the occurrence in S. otites of a locus that seems to be male-specific lethal raises the question of its maintenance in the species, and its possible effect on the evolution of its bearing chromosome. As the locus exhibits an antagonistic selective effect among sexes, could it favor the emergence of a new sex chromosome that ultimately would displace the ZW chromosome pair?

Given this rapid turnover that seems to occur in this subsection, composed of more than ten species that are dioecious (Oxelman et al. 2013), the investigation of sex determination system at the scale of the subsection will be necessary to have a better picture of sex chromosome origin and evolution, using the methodology used in the present study as a very efficient first step. Additional comparative genetic mapping of the dioecious species and a nondioecious sister species (e.g., S. nutans) will allow to identify not only the autosomal ancestor of sex chromosomes but also the possible causes of the loss of recombination between sex chromosomes, such as inversion or translocation events. For example, in S. latifolia, Qiu et al. (2016) were able to establish an evolutionary scenario of the sex chromosome involving three different LGs of nondioecious S. vulgaris. The study of this subsection offers the opportunity to compare the evolution of X/Z and Y/W chromosomes between dioecious species that are phylogenetically close (>4 Myr, Slancarova et al. 2013) in regards to theoretical expectations that predict different features such as the impact of mating systems on X/Z genetic diversity, the level of degeneration of Y/W, and the level of dosage compensation (Bachtrog et al. 2011; Ellegren 2011; Mullon et al. 2015).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.