Delineating the W Sex Chromosome in the Clam Shrimp, Eulimnadia texana

In many animals and plants, a specific region of the sex chromosomes (termed the “sex-determining region”) contains genes that code for an organism’s sex (male or a female) and sexual development [1, 2]. Sex chromosomal systems commonly occur in two forms: XY and ZW systems [3]. In XY systems, females are homogametic (XX) and males are heterogametic (XY). In ZW systems, males are homogametic (ZZ) and females are heterogametic (ZW).

Sex chromosomes have evolved independently numerous times in various animal and plant lineages [3‒6]. Of those lineages, a common characteristic of sex chromosomes is that the sex-limited (Y or W) chromosome is degraded. This represents an example of evolutionary convergence of sex chromosomes across distant lineages, which allows a unique opportunity to study the evolutionary forces that act on sex chromosomes [7].

Sex chromosomes begin their evolutionary trajectory as “proto-sex chromosomes” that have reduced recombination around a small sex-determining region to reduce the likelihood of crossing over between “sexually antagonistic” genes (genes that benefit one sex but are detrimental in the opposite sex), and by extension, reduce the likelihood of producing low-viability “intersexes” [7‒9]. Such crossover suppression should create a linkage unit near the sex-determining region of the proto-sex chromosomes [7, 10]. However, a side effect of recombination suppression is that transposable elements (TEs), a type of repeat sequence that moves from one location in the genome to another, are predicted to accumulate in the non-recombining region of the permanently heterozygous chromosome (Y or W) [11‒13]. This accumulation of TEs can cause gene mutations, deleterious effects in gene expression, and mutations in regulatory sequences due to their insertion into coding regions and replication within the genome [14‒16]. As a result, the non-recombining sex-determining region of sex chromosomes commonly contains high TE accumulation and low gene content compared to rest of the genome [17].

There are two major classes of TEs. Class I are retrotransposons that can move through copy-paste mechanisms using intermediate RNA and can be subdivided into long terminal repeat retrotransposons (termed “LTR”; e.g., gypsy) and non-LTR retrotransposons (e.g., L1 and ALU) [18, 19]. Class II retrotransposons are DNA transposons that can move via cut-and-paste mechanisms (e.g., hATs and Helitrons) [20]. Distribution patterns of TE types can be used as markers to identify the boundaries of the recombining and non-recombining regions of a sex chromosome. L1 and gypsy elements accumulate in heterochromatic regions of suppressed recombination, whereas Alu and DNA transposons elements accumulate in regions of high recombination [21‒24]. The accumulation of TEs can disrupt the functions of the genes present in the non-recombining region and as a result, the gene functions can be lost leading to functional and physical degradation of the Y(W) chromosome [13, 25‒27].

Studies of sex chromosomes in mammals, Drosophila, and Caenorhabditis elegans support theoretical predictions that the later stages in the evolution of permanently heterozygous sex chromosome (Y or W) accumulate mutations, have low GC content (indicative of low recombination rates), lose genetic information, and physically degrade [7, 28‒30]. However, because the sex chromosomes of these model organisms are already highly diverged, it is important to also explore new study systems with “young” sex chromosomes (i.e., at early evolutionary stages of sex chromosome evolution) [9, 31] to allow the delineation of the evolutionary processes involved in the earlier stages of sex chromosome divergence.

The freshwater crustacean Eulimnadia texana is an excellent candidate to serve as a new study system, as these “clam shrimps” have an androdioecious (males and hermaphrodites) breeding system [32, 33]. E. texana is a branchiopod crustacean known for its unique sex-determining system [34] having three common mating types formed by combining two “proto-sex chromosomes”: Z and W [34, 35]. Males are homozygous for the Z chromosome (ZZ), while hermaphrodites are ZW or WW. Unlike in most animals, E. texana expresses both the Z and W chromosomes in homogametic form, allowing the purging of deleterious recessive alleles from both chromosomal types, which should retard physical degradation of both the Z and W [35]. This unique feature of E. texana allows for the study of proto-sex chromosomes [36, 37], making it an ideal species for exploring early-stage sex chromosome evolution.

In a recent study, the E. texana WW genome was assembled producing a genome size of 120 Mb and 108 scaffolds [38]. Baldwin-Brown et al. [38] suggested that the putative sex chromosome was located in the largest scaffold (scaffold 1) of the genome. However, they were unable to identify the sex-determining region on scaffold 1. In a separate study, a low genome-wide recombination rate was observed in ZW hermaphrodites implying that recombination on the sex chromosome is suppressed in these hermaphrodites [39]. Baldwin-Brown and Long [39] proposed that only Z chromosomes undergo recombination in males [39], and by extension, that the W chromosome never experiences crossing over. However, a recent study of TE content in E. texana showed no signs of TE accumulation throughout the entire W chromosome relative to the autosomes [40], indicating that the W chromosome does not have crossover suppression throughout its length. Instead, Lang et al. [40] proposed that a smaller portion of the W had an accumulation of TEs, suggestive of a smaller region of crossover suppression.

The objective of our study is to assess whether E. texana’s scaffold 1 contains a smaller region of crossover suppression by testing for a TE “hotspot,” low gene density, and low GC content. All three indicators do indeed suggest that a smaller region (close to previously identified sex-linked markers) has low-to-no crossing over, indicating that E. texana’s W chromosome is indeed in a proto-sex chromosome stage of early evolutionary development.

The genomic data for WW E. texana were collected from NCBI GenBank, under BioProject “PRJNA352082” and genome accession number “NKDA00000000.” Microsatellite sequences of sex-linked markers were collected from NCBI, accession number AY515024.1 for microsatellite CS8 sequence, accession number “AY515026.1” for microsatellite CS11 sequence, and accession number “AY515028.1” for microsatellite CS15 sequence.

To identify species-specific repeat sequences, we used RepeatModeler (https://www.repeatmasker.org/RepeatModeler/), a de novo TE family identification and modeling package, with default parameters. Using the LTRStruct option in RepeatModeler, LTRs present in the repeat sequences were verified. To annotate the TEs in the genome, RepeatMasker (https://www.repeatmasker.org/RepeatMasker/) was used against RMBlast using the species-specific repeat library.

To annotate the TEs in the genome, RepeatMasker has built-in, statistically optimal scoring matrices to align the sequences and identify the hits with biological significance. For the analysis of the distribution of different types of TEs across different regions of the genome, a randomization test was used [41]. A window of 10,000 bp was used to break the genome into “segments.” Scaffold 1 contained 4,369 segments, whereas rest of the genome contained 7,839 segments. Here we considered scaffold 1 as sample 1 and rest of the genome as sample 2. For each segment in the two samples, the number of TEs present was identified, and the mean number of TEs present for each sample was calculated. The mean difference (“observed” mean difference) between sample 1 and sample 2 was calculated. Then we combined all the segments from the two samples into one pool and randomly drew 4,369 segments and assigned them to sample 1 and 7,839 segments to sample 2. The mean number of TEs presents for each randomized sample, and the mean difference (“expected” mean difference) between the two samples was then calculated. We repeated this expected difference calculation 10,000 times. The p value was then the percentage of times the observed mean difference value was equal to or greater than the expected mean difference value. If the p value is less than 0.05 (p < 0.05) then the TE accumulation on scaffold 1 was deemed significantly higher than the rest of the genome.

The above-mentioned randomization method was used a second time to compare regions within scaffold 1. Scaffold 1 was split into 20 equally sized “bins” to determine if the TEs were uniformly distributed across the scaffold or were instead concentrated into one or more “hotspot” region(s). Herein, each bin of scaffold 1 was treated as sample 1 and rest of the bins of scaffold 1 plus the rest of the genome was treated as sample 2. The likelihood that a bin had a higher observed mean difference than expected was calculated as noted above.

Microsatellite sex-linked marker sequences were blast against the scaffold 1 sequence to identify the exact genomic positions of the sex-linked markers on scaffold 1. These sex-linked markers were then mapped onto scaffold 1 using R scripts.

Known TEs ≥ 500 bp were mapped onto the genome. The information on TE type and their locations were obtained from the RepeatMasker output files and were then mapped onto the genome using R scripts.

To generate information on gene density and GC content of scaffold 1, we used the genome FASTA and the gene annotation GTF files of the WW genome, which were available from WW genome assembly project [38]. First, the genome was indexed using SAMtools, and then we extracted the sequence lengths of each scaffold from the output.fai file. We created a 2,200 kb sliding window file using BEDTools in order to break scaffold 1 into 20 bins which allowed for the comparison of gene density and GC content in those bins. Second, we generated a BED file storing coordinates of coding regions in using a GTF file. Third, we calculated the fraction of coding sequences in the genomic regions using the coverage BEDTool. Finally, using the same 2,200 kb genome sliding window file and BEDTools “nuc” option, we calculated the GC content of the genomic regions. Then gene density and GC content values for scaffold 1 were extracted from files containing the fraction of coding sequences in the genomic regions and the GC content file, respectively. To visualize gene density and GC content present in the 20 bins of scaffold 1, R scripts were used.

To visualize the GC content of scaffold 1 in a boxplot and conduct pairwise Wilcoxon tests, we generated a GC file with a 50 kb sliding window using the above-mentioned method. From the 50 kb window GC content file, we extracted the information on window coordinates, percentage values of AT and GC and added a new column called “location.” Here, the coordinates which fall within the putative sex-determining region were labeled as “WSDR” (W sex-determining region), and rest of the coordinates were labeled “PAR” (pseudoautosomal region), under the location column. The boxplot generation and pairwise comparisons for genomic location versus GC content using Wilcoxon rank sum test were conducted using an R script.

To identify gene models on scaffold 1, we utilized the genome annotation file produced by the Maker pipeline [42]. The protein sequences derived from the MAKER pipeline for WW individuals were subjected to a protein BLAST search [43] against the Drosophila proteome database, which was obtained from UniProt database [44]. The BLAST results were filtered based on the following criteria: e-value <1e−10, bit score >50, and percentage identity ≥80%. The genes associated with sex determination, spermatogenesis, and oogenesis were extracted from UniProt and then identified on scaffold 1. Finally, we calculated the number of genes found in each bin of scaffold 1.

The RepeatMasker scan of the WW genome of E. texana found that 30.8% of the genome was composed of TEs (online suppl. Fig. 1; for all online suppl. material, see https://doi.org/10.1159/000542284). The majority of the repeat content was composed of unknown/unclassified TEs, whereas DNA transposons and retrotransposons comprised ∼4% and ∼10% of the genome, respectively. LTRs were the dominant TE family compared to other known TE families (online suppl. Fig. 1).

To localize the putative sex-determining region of the WW genome, the known TE accumulation was compared on scaffold 1 to the rest of the genome. Scaffold 1 was not significantly higher in TEs than the rest of the genome (observed mean difference = −0.418, p value = 1). Therefore, we split scaffold 1 into 20 equal sections (“bins”) and carried out another set of randomization tests, treating each bin as an individual sample to compare to the rest of the genome. The genomic region of 6.6 Mb through 8.8 Mb of scaffold 1 (fourth bin, Table 1) showed significantly higher accumulation of TEs (observed mean difference = 0.317, p value = 0.013), while other bins of the scaffold 1 showed no significant accumulations of TEs compared to the rest of the genome. Therefore, from this point onwards, we will refer to E. texana’s bin 4 as the putative WSDR.

All three sex-linked marker microsatellites were found on scaffold 1 (Fig. 1a): CS8 was located at ∼2.1 Mb, CS11 was located at ∼7.4 Mb (within the above-mentioned WSDR), and CS15 was located at the center of the scaffold (∼22 Mb).

On scaffold 1, TE classes were enriched (above 20% of TEs) in the first quarter portion, while TEs along the remainder of the scaffold were scattered (Fig. 1b). At the last quarter portion of the scaffold, DNA TEs had high peaks (above 20%) and LINE had medium peaks (10–20% of TEs) (Fig. 1b). Interestingly, there were two peaks of LTR (indicated by dense purple peaks above 10% of TEs) in the region of 5–10 Mb (Fig. 1b). A similar enrichment pattern was observed in LTR and LINE repeats when the TEs above 500 bp were mapped separately to scaffold 1 (online suppl. Fig. 2a). DNA repeats showed random distribution along scaffold 1; they were enriched even in the regions where LTR and LINE repeats were poor (online suppl. Fig. 2a). LINE/I and LTR/Gypsy repeats, the two most abundant TE superfamilies on scaffold 1, tended to occur in the same genomic regions displaying similar enrichment patterns (online suppl. Fig. 2b).

The gene density along the length of scaffold 1 was found to be nonuniform, ranging from a low of 0.308 to a high of 0.452 (Fig. 1c). Notably, bin 4, which corresponds to the putative WSDR, displayed the lowest gene density (0.308), while bin 19 exhibited the highest gene density (0.452) (Fig. 1c). Analysis of sex determination-related genes across scaffold 1 revealed that the first four bins showed a higher concentration of genes associated with reproductive processes (drawn from Drosophila melanogaster database) compared to the remaining bins (online suppl. Table 1).

Similar to gene density patterns, the levels of GC content along scaffold 1 were found to be nonuniform, ranging from 0.382 to 0.412 (Fig. 1c), with bin 7 having the lowest GC content (0.382), bin 4 having second lowest GC content (0.384), and bin 20 having the highest GC content (0.412) (Fig. 1c). When comparing GC content of PAR (pseudoautosomal region) and WSDR (bin 4) genomic locations on scaffold 1, WSDR had a lower GC content than PAR region, and the Wilcoxon rank sum test revealed that the difference was significant (Fig. 2, p value = 0.00001).

Many studies provide evidence that TEs play a pivotal role in sex chromosome evolution [45, 46]. Moreover, the non-recombining sex-determining regions are found to be rich in TEs along with having low gene density and low GC content [17, 47]. In the current study, we used information on TE distribution, sex-linked markers, gene density, and GC content to delineate the likely WSDR of E. texana.

In agreement with an earlier study of TE distribution in E. texana [40], three sex-linked markers were found on scaffold 1 reinforcing the notion that this is indeed the sex chromosome. Contrary to the expectation of a uniform TE enrichment across scaffold 1 [39], we found that scaffold 1 has sparse accumulation of TEs throughout most of its length (Fig. 1b, online suppl. Fig. 2a, b). Scaffold 1 does not have a significantly higher TE content throughout its length relative to the rest of the genome, suggesting that the entirety of scaffold 1 is unlikely to be non-recombining. These results agree with the findings of Lang et al. [40].

In addition to the TE predictions, if scaffold 1 is non-recombining across its entirety, we would expect a low gene density [17] and a low GC content across the entire length of scaffold 1 [47, 48]. Our results revealed that the gene density and GC content levels varied across the length of scaffold 1. While some regions of scaffold 1 had high gene densities, others had low gene densities (Fig. 1c). Similarly, some regions of scaffold 1 had high GC content, whereas other regions had low GC content (Fig. 1c). Both results indicate that the entire length of scaffold 1 is not subjected to recombination suppression, again reinforcing previous findings [40]. Therefore, although scaffold 1 harbors sex-linked markers, it is not likely to be a fully non-recombining chromosome, as previously predicted [39].

Lang et al. [40] suggested that a smaller region of E. texana’s scaffold 1 might harbor the sex-determining region, which has been found in other organisms [49‒51]. However, the number of TEs masked in the Lang et al. [40] study was too small to test this hypothesis. In the current study, the inclusion of species-specific TEs allowed us to test for evidence of a smaller region of scaffold 1 harboring a non-recombining region. Scaffold 1 was divided into 20 equal-sized section (“bins”) to investigate TE accumulation in each section against the rest of the genome. Our results showed that only bin 4 (∼5% of scaffold 1 ranging from 6.6 to 8.8 Mb) was found to be accumulating higher TE content compared to the rest of the genome. None of the remaining bins had higher TE accumulation compared to the rest of the genome (Table 1). Therefore, higher TE accumulation was limited to a small region of scaffold 1, implying that the non-recombining region of the W is concomitantly small, as predicted by Lang et al. [40]. Two of the three previously found [35] sex-linked markers also mapped close to this likely non-recombining region (blue in Fig. 1a). Additionally, the gene density analysis (Fig. 1c) revealed that bin 4 has the lowest gene density among the 20 bins of E. texana’s scaffold 1, indicating this region has a high level of gene degradation. However, the WSDR (Bin 4) retained several genes associated with reproductive processes in Drosophila melanogaster (online suppl. Table 1), including the mts gene (involved in oogenesis [52]), the Rab6 gene (plays a role in germarium-derived egg chamber formation, oocyte microtubule cytoskeleton polarization, and pole plasm oskar mRNA localization [53]), and the dac gene (associated with genital disc sexually dimorphic development [54]). Similar to our results, the smaller sex-determining regions in young sex chromosomes of spinach, papaya, and Mercurialis annua were also found to exhibit patterns of high repeat content and lower gene density [55‒57]. Therefore, the observed patterns of TE distribution, linked genes, and gene content in E. texana all suggest that a smaller region of scaffold 1 is the likely sex-determining region.

GC content serves as a useful proxy for identifying non-recombining regions [29, 48]. Many species (e.g., mammals, Drosophila melanogaster, Caenorhabditis elegans, and yeast) have low GC content in non-recombining regions [29, 48]. Researchers have observed significantly lower GC content on degenerated sex-limited W(Y) chromosomes of chickens, zebra finches, humans, and mice when compared to autosomes and Z(X) chromosomes [58]. Our results for GC content revealed that bins 7 and 4 have the lowest GC content among the 20 bins of scaffold 1 (Fig. 1c), suggesting that these two bins might be non-recombining regions. Although the GC content is low in bin 7 (Fig. 1c), its gene density was relatively high compared to bin 4 (Fig. 1c) and its TE content was relatively low compared to the entire genome (Table 1), both of which are inconsistent with bin 7 being the sex-determining region of the W. However, the association between GC content and gene density is complex and not well understood [30]. A study of the rattlesnake W chromosome found that the youngest evolutionary strata were relatively gene-rich despite recombination suppression [58]. Therefore, bin 7 might be a young evolving stratum with recombination suppression and thus exhibits a lesser degree of degradation than the potentially older evolutionary stratum, bin 4. Moreover, our results on comparison of WSDR (bin 4) with PAR (the rest of the bins in scaffold 1) revealed that WSDR has significantly low GC content compared to PAR (p value = 0.00001, Fig. 2). Overall, our findings from analyses of TE accumulation, gene density, and GC content suggest that bin 4 likely harbors the putative sex-determining region of the proto-W sex chromosome.

As mentioned earlier, many studies on sex chromosome evolution have been done on model organisms with advanced sex chromosomal degradation, which provide much information on the later stages of degeneration of non-recombining chromosomes [7, 28]. However, there is little understanding of the evolutionary forces driving recombination suppression and proto-sex chromosome differentiation during the early stages of sex chromosome evolution [59]. To overcome these challenges, it is essential to explore species (e.g., E. texana) that have “young” sex chromosomes that are at the early stages of sex chromosome differentiation [9].

In conclusion, this study confirms previous findings [40] that challenge the prediction that Scaffold 1 represents a wholly non-recombining sex chromosome [39]. Instead, we have evidence of a smaller sex-determining region embedded within the W chromosome that exhibits low recombination, in agreement with the idea that E. texana’s W is a proto-sex chromosome [35, 40]. Hence, E. texana is an ideal candidate for studying evolutionary processes that drive early sex chromosome evolution. Future studies of E. texana should explore the genomic composition of the proto-sex-determining region and the evolution of recombination suppression, as well as whether LTR/Gypsy TEs contribute to structural changes in these sex chromosomes [60]. To explore this role, TE composition and sex-specific TE abundance on ZZ and WW sex chromosomes would need to be compared [60]. In the near future, we plan to assemble the genome of ZZ males so that we can compare the TE composition of the Z to the W chromosomes to understand how TEs may influence sex chromosome differentiation in proto-sex chromosomes.

We thank following people for their continuous support and guidance: Agnieszka Lipinska, Octavio Manuel Palacios Giménez, and Sean Mitchuson.

All work was conducted on nonregulated invertebrate animals; thus, neither IACUC reviews nor a review by the Ethical Review Board were necessary. An ethics statement was not required for this animal study as the species used does not require ethics approval according to local/national guidelines.

The authors have no conflicts of interest to declare.

No external funding was used in the conduct of these studies.

C.E. designed the study, performed data collection, analysis, and interpretation, and wrote the manuscript. S.C.W. helped conceive the study and edited the manuscript.

All data generated or analyzed during this study are included in this article and its online supplementary material. Further inquiries can be directed to the corresponding author.