https://orcid.org/0000-0002-4085-3804
Abstract
Phenotypic invariance—the outcome of purifying selection—is a hallmark of biological importance. However, invariant phenotypes might be controlled by diverged genetic systems in different species. Here, we explore how an important and invariant phenotype—the development of sexually differentiated individuals—is controlled in over two dozen species in the frog family Pipidae. We uncovered evidence in different species for 1) an ancestral W chromosome that is not found in many females and is found in some males, 2) independent losses and 3) autosomal segregation of this W chromosome, 4) changes in male versus female heterogamy, and 5) substantial variation among species in recombination suppression on sex chromosomes. We further provide evidence of, and evolutionary context for, the origins of at least seven distinct systems for regulating sex determination among three closely related genera. These systems are distinct in their genomic locations, evolutionary origins, and/or male versus female heterogamy. Our findings demonstrate that the developmental control of sexual differentiation changed via loss, sidelining, and empowerment of a mechanistically influential gene, and offer insights into novel factors that impinge on the diverse evolutionary fates of sex chromosomes.
Introduction
Developmental Systems Drift and Sex Determination
An important discovery is that major developmental cascades are governed by conserved suites of genes, suggesting the existence of a “genetic toolkit” whose components orchestrate core developmental processes across diverse organisms (e.g., Pax6/eyeless in eye development Xu et al. 1997; distal-less in limb development Cohn and Tickle 1999). This conservation allows us to understand fundamental biological mechanisms by studying a handful of model organisms. However, conserved phenotypes may be controlled by genetic systems that differ among species. For example, there is variation among closely related species in transcription factor binding positions despite conservation of regulatory function (Villar et al. 2014). This divergence in genetic systems that underpin conserved phenotypes is called developmental systems drift (DSD) (Weiss and Fullerton 2000; True and Haag 2001). How DSD occurs remains an open question, and could involve pathway switching, convergence, and rapid evolution (Haag and Doty 2005), changed pleiotropic interactions (Pavlicev and Wagner 2012), and neutral variation could play an important role (Lynch and Hagner 2015).
Because it is fundamentally linked to reproduction and therefore fitness, it is surprising that the genetic control of two differentiated sex phenotypes (male and female) is so frequently subject to DSD. For example, differences in which sex is the heterogametic sex (male for XY systems and female for ZW systems) indicates DSD of the regulation of sexual differentiation in amphibians, fishes, and reptiles (Evans et al. 2012; Pennell et al. 2018). This is also true at finer phylogenetic scales, such as within gecko lizards (Gamble et al. 2015), stickleback fishes (Ross et al. 2009), and ranid frogs (Jeffries et al. 2018). Even in mammals, several species have independently lost the male determining gene SRY along with the rest of their Y chromosomes (Sutou et al. 2001; Matveevsky et al. 2017). More broadly, a diversity of environmental cues, which impinge on and thus are coupled to variable genetic systems, can govern sexual differentiation (McCabe and Dunn 1997; Refsnider and Janzen 2016).
Diverse sex-determination systems differ in the ways that genes act to trigger sexual differentiation (such as differences between the sexes in gene dosage, splice variants, gene content, and the ratio of X to autosomal chromosomes) and also in what genes actually trigger sexual differentiation (Graves 2008). Although triggers of sexual differentiation can evolve rapidly, several genes involved with the developmental process of sexual differentiation have conserved roles (Angelopoulou et al. 2012). An example of a highly conserved sex-related gene is dmrt1; this gene is a transcription factor that contains a DNA-binding DM-domain and operates in a downstream capacity in sexual differentiation of eutherian mammals (Matson et al. 2011), but is the sex-determining locus in birds (Smith et al. 2009; Lambeth et al. 2014), and (as a paralog) in various fish and amphibians (Matsuda et al. 2002; Yoshimoto et al. 2008). Conservation of these sex-related genes suggests that DSD often involves one central system with distinctive inputs in different species (i.e., variation on a theme), as opposed to wholesale replacement of major systems (Graves and Peichel 2010).
How Does DSD Influence Sex Chromosomes?
Sex chromosomes are distinguished from each other and from autosomal pairs by the presence of genetic variation that triggers (and thus is associated with) sexual differentiation. Some sex chromosome pairs are nearly identical in sequence and structure (homomorphic), whereas others evolve additional differences in gene content and repetitive elements, that causes them to be cytogenetically distinctive from each other (heteromorphic). Sex chromosome heteromorphy is caused by Hill-Robertson effects (Hill and Robertson 1966; Felsenstein 1974) that stem from recombination suppression over portions of the sex chromosomes for an extended period of evolutionary time (Charlesworth and Charlesworth 2000). In homomorphic sex chromosomes, recombination suppression is presumably restricted to a small region of the sex chromosomes and/or for a short period of time. Whether sex chromosomes evolve to be heteromorphic or stay homomorphic (the presumptive ancestral-state; Bachtrog et al. 2014) may depend on several factors. For example, de novo evolution or translocation of a sex-determining gene could convert a pair of autosomal chromosomes into homomorphic sex chromosomes (Volff et al. 2007). The “rapid turnover” explanation for homomorphic sex chromosomes predicts that diverged species with homomorphic sex chromosomes should have nonhomologous genomic locations of their triggers for sex determination. Rapid turnover could be favored or facilitated by temporal or geographic variation in the selective advantage of differing offspring sex ratios (West et al. 2000), neutral variation in genetic triggers for sex determination (Bull and Charnov 1977), or a sex-specific fitness advantage of a new sex-determining system (Bull and Charnov 1977).
Another possibility is that sex chromosome homomorphy could exist if recombination is common enough between the sex chromosomes with sufficient frequency to prevent divergence between them (Perrin 2009), but not so common as to frequently disrupt the sex-determining locus. The advantages of a high recombination rate might still be enjoyed in non-sex linked portions of the genome if the rate of recombination is sex biased (heterochiasmy or achiasmy) with a higher rate in the homogametic sex (females for XY systems, males for ZW systems). Under this mechanism, the “fountain of youth” explanation for homomorphic sex chromosomes, does not require variation among species in the genomic location of the trigger for sex determination, and predicts a lower rate of recombination rates in the heterogametic sex (males for XY systems, females for ZW systems).
Intra-locus sexual antagonism—when a mutation has opposite relative fitness effects on each sex—could also contribute to sex chromosome divergence in several ways. For instance, recombination suppression should be favored by natural selection when it links alleles with sexually antagonistic function to the phenotypic sex that these alleles benefit (Jordan and Charlesworth 2012). Thus, homomorphic sex chromosomes may exist in species in which sexual antagonism is rare or resolved by mechanisms other than physical linkage to a sex-determining locus, such as by the evolution of sex-biased expression. Under this scenario, different species with homomorphic sex chromosomes may not vary in the genomic location of their trigger for sex determination, but if this is the case, sexual antagonism should either be rare, or resolved in other ways.
Conversely, it is also possible that the origin of a mutation with sexually antagonistic effects on an autosome could drive turnover of sex chromosomes (van Doorn and Kirkpatrick 2010), in which case expectations would match those of the rapid turnover hypothesis. The “hot potato” model considers the distinctive responses of natural selection to mutations with sexually antagonistic effects, and to deleterious mutations that accumulate in nonrecombining regions of a sex chromosome that contain the trigger for sex determination (Blaser et al. 2014). Natural selection on the former type of mutation is expected to favor the expansion of recombination suppression surrounding sex-linked regions, whereas selection on the latter type of mutation is expected to favor the establishment of new sex chromosomes with small nonrecombining regions. In this way, the hot potato model predicts that transitions in sex chromosomes will occur continuously through evolution (Blaser et al. 2014). This cycle could be broken by “evolutionary traps” that disfavor sex chromosome transitions, such as when genes that are essential for development of one sex become sex-linked (Bachtrog et al. 2014). Divergence of sex chromosomes may also be limited if an inability to cope with sex-differences in gene dosage (i.e., hemizygosity in one sex) make the selective cost of sex chromosome degeneration too high (Adolfsson and Ellegren 2013).
Pipid Frogs and DSD of Sex Determination and Sex Chromosomes
Species in the frog family Pipidae are a compelling group with which to explore possible links between genes that trigger sex determination and the evolutionary fate of the sex chromosomes on which they reside. Pipids have considerable variation in the genomic locations of genes that trigger sexual differentiation, and also have variation among and within species in male versus female heterogamy, and in the presence or absence of at least one trigger for sex determination (dm-w; Yoshimoto et al. 2008; Olmstead et al. 2010; Bewick et al. 2011; Roco et al. 2015; Furman and Evans 2016). There is also variation among pipid species in the extent of recombination suppression in genomic regions linked to different sex-determining loci (Bewick et al. 2013; Furman and Evans 2018), even though morphologically diverged sex chromosomes have not been observed (Wickbom 1950; Tymowska 1991).
Using genomic approaches, we identify here new variation in the developmental systems that underpin sexual differentiation in pipid frogs by demonstrating that dm-w has variable effects on female differentiation in different species. We also find new variation in the degree to which recombination is suppressed in regions linked to sex-determining genes. By studying sex chromosomes in this group, we elucidate three distinct ways in which DSD can happen, including loss, sidelining, or empowerment of a mechanistically influential gene. Through interactions with genes with sexually antagonistic functions, we speculate that variation in the efficiency of genetic triggers for sex determination is a key determinant of the evolutionary paths taken by sex chromosomes.
Results
Loss, Sidelining, and Empowerment of an Imperfect Regulator of Sexual Differentiation
In the African clawed frog Xenopus laevis, the gene dm-w exists as a single allele in females and triggers their sexual differentiation, and is a partial duplicate of the gene dmrt1, which is involved in the sex-determining cascades of many species (Yoshimoto et al. 2008). Dm-w is also known from several species that are closely related to X. laevis (Bewick et al. 2011), and is distinguished from dmrt1 by nucleotide divergence and several insertion/deletions (see Supplementary Material online for details). These features allowed us to use targeted next-generation sequencing of genomic DNA to identify dm-w in assembled reads from a panel of Xenopus species. Our probes targeted regions of dm-w that are homologous to dmrt1 (i.e., dm-w exons 2 and 3) but not regions of dm-w that are not homologous to dmrt1 (i.e., dm-w exons 1 and 4).
In these targeted next-generation sequences, we detected dm-w exons 2 and 3 in a paraphyletic assemblage of female individuals from seven tetraploids, two octoploids, and one dodecaploid, and failed to detect dm-w in females of one diploid, seven tetraploids, three octoploids, and three dodecaploids (fig. 1). There was perfect concordance in the identification of independently assembled dm-w exons 2 and 3 (i.e., if dm-w was detected at all in an individual, both exons were independently detected in that individual). We did not detect missense or nonsense mutations in any of the dm-w sequences we recovered from exon 2 or 3.
The sex-determining gene dm-w was detected in at least one individual (red, pink, light blue with dagger) using targeted next-generation sequencing, but was not detected in at least one individual in other species (light blue, black, white). In some species, dm-w was detected, but was not female-specific based on PCR assays (daggers). Data are plotted on a Bayesian consensus phylogeny estimated from complete mitochondrial genome sequences, as described in Supplementary Material online. All nodes have 100% posterior probability except where labeled. The most recent common ancestor of all species in which dm-w was detected, indicated with an asterisk, has a 94% likelihood of having this gene not fixed only in females and a 5% likelihood of it being fixed only in females. The scale bar indicates the time in millions of years ago.
As expected because dm-w originated in an ancestor of subgenus Xenopus after divergence of an ancestor of subgenus Silurana (Bewick et al. 2011), we did not identify dm-w in the targeted next-generation sequences from any of the species in subgenus Silurana. Included in our capture experiment were several male individuals of multiple species (supplementary table S1, Supplementary Material online). As expected, males from two species in subgenus Xenopus (X. boumbaensis and X. lenduensis) did not have either dm-w exon. Surprisingly, however, males of two other species (X. itombwensis and X. pygmaeus) carried both exons of dm-w that were targeted by our probes.
Our results are congruent with findings based on PCR assays by Bewick et al. (2011), and extend them by 1) ruling out a role of primer mismatch in the failure to detect dm-w, 2) independently validating all detections with two exons, and 3) surveying several species that were discovered recently, including X. kobeli—the only dodecaploid species with dm-w. Thus, dm-w is present in one or more tetraploid, octoploid, and dodecaploid species. Bewick et al. (2011) detected dm-w exon 2 in a female X. pygmaeus, and here, we detected exons 2 and 3 in a male individual of this species.
Several lines of evidence in addition to the perfect concordance with the PCR assays of Bewick et al. (2011) allow us to conclude that the targeted next-generation sequence results were minimally influenced by sequencing depth, cross hybridization with nontarget genomic regions, or rapid divergence of dm-w in some species. First, the sequencing depth of capture libraries of samples in which dm-w was detected was not substantially different from (and on an average lower than) that of samples that where dm-w was not detected (3.4–5.6 million or 3.4–23.6 million sequence reads were captured from libraries with and without dm-w, respectively, supplementary table S1, Supplementary Material online). Second, we were able to detect both expected dmrt1 ohnologs in all tetraploid species for both dmrt1 exons, including allelic variation for several loci, and an additional paralog of dmrt1 in X. borealis that arose after allotetraploidization and that was confirmed independently using Sanger sequencing (data not shown). This indicates that our approach efficiently captured sequences in each genome that were similar to our probes. We also identified most expected paralogs in all octoploid and dodecaploid species, although in most of these species, less than four or six paralogs (respectively) are expected due to gene loss (Bewick et al. 2011). Third, we compared the assembled captured sequences to other DM-domain containing genes from X. tropicalis (dmrt1, dmrt2, dmrt3, and dmrt5) and found them to be highly diverged in exonic and flanking intronic regions from the dmrt1 and dm-w sequences that we assembled from the captured sequences (data not shown). Fourth, for essentially all dm-w and dmrt1 sequences, we also captured and assembled flanking intronic sequences that were clearly homologous to the flanking intronic sequences from the genome sequences of X. laevis and X. tropicalis (Hellsten et al. 2010; Session et al. 2016), even though our capture probes did not target these intronic sequences. This indicates that, as expected, these intronic sequences are contiguous with exonic sequences of dmrt1 and dm-w that matched our probes. And fifth, previous efforts provide unambiguous evidence that at least one species—X. borealis—does not use dm-w to trigger sex determination: whole-genome sequences from one female and one male individual do not contain this gene, and the newly identified sex chromosomes of X. borealis (chromosome 8 L) are not the same as the chromosome on which dm-w resides in X. laevis (chromosome 2 L) (Furman and Evans 2016, 2018). That at least one Xenopus species appears to have lost dm-w and evolved new sex chromosomes increases our prior expectations that this could happen in other species as well.
The region of dm-w that we sequenced using targeted next-generation sequencing is relatively rapidly evolving. For example, uncorrected pairwise divergence between dm-w exon 2 from X. laevis and X. clivii is 8% at the nucleotide level (16 out of 209 ungapped bp) and 13.4% at the amino acid level (9 out of 67 ungapped amino acids). Uncorrected pairwise divergence of dm-w exon 3 from these two species is 11.0% at the nucleotide level (16 out of 145 ungapped bp) and 12.5% at the amino acid level (6 out of 48 ungapped amino acids). Thus, the failure to detect dm-w in some species (apart from X. borealis) using targeted next-generation sequencing could conceivably be false negatives due to rapid evolution of dm-w. However, if this were the case in some other species (e.g., X. longipes), then it would have to affect the entire DM domain (which is encompassed by our data from exons 2 and 3). To the extent that well supported relationships among mitochondrial genomes are congruent with evolutionary relationships among the primarily maternally inherited dm-w, false negatives due to rapid evolution of dm-w also would have to have occurred independently several times (i.e., separately in X. parafraseri, X. allofraseri, X. wittei, and more) and in some cases over rather short periods of time (e.g., X. ruwenzoriensis; fig. 1). As a conservative measure that accommodates this possibility of false negatives, the itemization of new sex-determination systems in the Discussion does not consider the sex-determining systems of species with undetected dm-w (apart from X. borealis) as necessarily distinctive (although they may be).
These results prompted us to turn to museum specimens of wild caught animals, and laboratory animals from six species, to test, using PCR and Sanger sequencing, whether dm-w was restricted to animals that were phenotypically female. We found that dm-w was female-specific in two species, including X. laevis, but not in four others, including X. clivii, which is the most divergent species from X. laevis that carries dm-w (table 1 and fig. 1). For three species (X. clivii, supplementary fig. S7, Supplementary Material online; X. pygmaeus, supplementary fig. S6, Supplementary Material online; and X. victorianus), dm-w was most commonly observed in phenotypic females, but it was also occasionally found in phenotypic males. For one of these species (X. victorianus), dm-w was found in all females (n = 20), but also found in one-third of the males (5 out of 15 tested). For two of these three species (X. clivii, X. pygmaeus), dm-w was not detected in roughly one-third of the females we surveyed (5 out of 16 female X. clivii and 3 out of 9 female X. pygmaeus did not have detectable dm-w) and dm-w was detected in several males (table 1). In one species (X. itombwensis), dm-w was observed in all individuals of both phenotypic sexes (present in 25 individuals including 20 males, supplementary fig. S3, Supplementary Material online). Furthermore, in X. itombwensis dm-w, heterozygous positions were observed in a few individuals of both phenotypic sexes (supplementary fig. S4, Supplementary Material online). Heterozygous positions were also observed in dm-w carried by five male and one female X. victorianus (supplementary table S4 and fig. S5, Supplementary Material online).
Dm-w Is and Is Not Female Specific (Y or N, respectively) in Several Xenopus Species.
| Species | No. of Females (no. with dm-w) | No. of Males (no. with dm-w) | Female-Specificity of dm-w | Notes |
|---|---|---|---|---|
| X. itombwensis | 5 (5) | 20 (20) | N | a |
| X. pygmaeus | 9 (6) | 11 (2) | N | |
| X. clivii | 16 (11) | 29 (5) | N | b |
| X. victorianus | 20 (20) | 15 (5) | N | c |
| X. laevis | 24 (24) | 12 (0) | Y | |
| X. gilli | 13 (13) | 7 (0) | Y |
| Species | No. of Females (no. with dm-w) | No. of Males (no. with dm-w) | Female-Specificity of dm-w | Notes |
|---|---|---|---|---|
| X. itombwensis | 5 (5) | 20 (20) | N | a |
| X. pygmaeus | 9 (6) | 11 (2) | N | |
| X. clivii | 16 (11) | 29 (5) | N | b |
| X. victorianus | 20 (20) | 15 (5) | N | c |
| X. laevis | 24 (24) | 12 (0) | Y | |
| X. gilli | 13 (13) | 7 (0) | Y |
Note.—Targeted next-generation sequencing, PCR, and Sanger sequencing were used to assess how many phenotypic females and males (no. of females and no. of males, respectively) carry dm-w in parentheses for each phenotypic sex.
dm-w is fixed or almost fixed; heterozygotes observed in both sexes suggesting autosomal segregation.
Includes 24 wild samples analyzed in Furman and Evans (2016).
Heterozygotes observed in one female and five males from eastern DRC.
Dm-w Is and Is Not Female Specific (Y or N, respectively) in Several Xenopus Species.
| Species | No. of Females (no. with dm-w) | No. of Males (no. with dm-w) | Female-Specificity of dm-w | Notes |
|---|---|---|---|---|
| X. itombwensis | 5 (5) | 20 (20) | N | a |
| X. pygmaeus | 9 (6) | 11 (2) | N | |
| X. clivii | 16 (11) | 29 (5) | N | b |
| X. victorianus | 20 (20) | 15 (5) | N | c |
| X. laevis | 24 (24) | 12 (0) | Y | |
| X. gilli | 13 (13) | 7 (0) | Y |
| Species | No. of Females (no. with dm-w) | No. of Males (no. with dm-w) | Female-Specificity of dm-w | Notes |
|---|---|---|---|---|
| X. itombwensis | 5 (5) | 20 (20) | N | a |
| X. pygmaeus | 9 (6) | 11 (2) | N | |
| X. clivii | 16 (11) | 29 (5) | N | b |
| X. victorianus | 20 (20) | 15 (5) | N | c |
| X. laevis | 24 (24) | 12 (0) | Y | |
| X. gilli | 13 (13) | 7 (0) | Y |
Note.—Targeted next-generation sequencing, PCR, and Sanger sequencing were used to assess how many phenotypic females and males (no. of females and no. of males, respectively) carry dm-w in parentheses for each phenotypic sex.
dm-w is fixed or almost fixed; heterozygotes observed in both sexes suggesting autosomal segregation.
Includes 24 wild samples analyzed in Furman and Evans (2016).
Heterozygotes observed in one female and five males from eastern DRC.
A caveat to a lack of positive PCR amplification of dm-w in females of some species is that there could be polymorphism at primer sites that prevented amplification in some individuals. However, in females lacking dm-w, the negative amplifications were independently replicated with two separate amplicons that targeted different and nonoverlapping portions of dm-w using different pairs of primers. We only considered samples for which a positive control (i.e., the amplification of a portion of the mitochondrial 16S ribosomal RNA gene) was successful. This caveat notwithstanding, we observed positive independent amplifications of multiple independent regions of dm-w in several males from multiple species (table 1). This provides unambiguous evidence that dm-w is not completely linked to female differentiation in some species, even though it is usually found in females in these same species. Ancestral-state reconstruction of these observations suggests that dm-w was not fixed only in females in the most recent common ancestor of species that carry this gene (fig. 1). Further details about these methods and results are discussed in Supplementary Material online.
Additional Variation in Sex-Linkage, Recombination Suppression, and Heterogamy in Other Pipid Species
Using reduced representation-genome sequencing (RRGS) and shotgun whole-genome sequencing (WGS) from captive-bred families, we extended our analysis to the sex chromosomes of three other pipid species whose ancestors never carried dm-w. This effort also uncovered variation in sex chromosomes, and definitively demonstrates that males are heterogametic in Hymenochirus boettgeri and Pipa parva, and suggests that females are heterogametic in X. mellotropicalis. The sex chromosomes of H. boettgeri and P. parva correspond to X. tropicalis chromosomes 4 and 6, respectively, but the sex chromosomes of X. mellotropicalis correspond to X. tropicalis chromosome 7. Additionally, based on the genomic locations where scaffolds containing sex-linked SNPs in the RRGS data mapped to the X. tropicalis genome assembly, the extent of recombination suppression on the sex chromosomes varied substantially among these species, with a large region in H. boettgeri (∼79.8 Mb, i.e., ∼60% of the assembled chromosome, containing 11 sex-linked SNPs out of 29 SNPs total in the sex-linked region), a medium-sized region in P. parva (at least ∼8.6 Mb, containing 3 sex-linked SNPs out of 5 SNPs total in the sex-linked region), and a small region in X. mellotropicalis (only 1 scaffold <1 Mb contained one sex-linked SNP out of 5 SNPs total on this scaffold) (fig. 2). Because these sex-linked SNPs in each family are based on RRGS data, they represent a small fraction of the total number of sex-linked SNPs in the genomes of each species. Embedded within these inferred regions of recombination suppression, we did observe several SNPs in the RRGS genotypes that did not have a completely sex-linked pattern of inheritance. We suspect that these are genotype errors stemming from low coverage because they are flanked by completely sex-linked SNPs that were confirmed by Sanger sequencing described below.
The sex chromosomes of Hymenochirus boettgeri (Chr.04), Pipa parva (Chr.06), and Xenopus mellotropicalis (Chr.07) are not homologous. For each species, orange or black dots represent the homologous genomic locations of sex-linked or not sex-linked SNPs, respectively, on the ten chromosomes of X. tropicalis. The density of SNPs is highest for X. mellotropicalis because it is most closely related to X. tropicalis and because scaffolds from both subgenomes of this allotetraploid species map to only one region of the diploid X. tropicalis reference genome. Numbers refer to genomic regions that were validated by Sanger sequencing, and are in blue, or orange font based on whether the Sanger sequences had completely, or partially sex-linked SNPs, respectively, for the following loci: (1: sall1, 2: dmrt5, 3: hmcn1, 4: kctd1, 5: ncoa2, 6: mmp16, 7: or8h1). Putative false negatives for sex-linkage for H. boettgeri and P. parva (some black dots on chromosomes 4 and 6, respectively), and a false positive for sex-linkage for X. mellotropicalis (an orange dot on chromosome 8) are discussed in Supplementary Material online.
Our inferences concerning the extent of recombination suppression in each focal species are based on the locations of genomic regions in X. tropicalis that are homologous to dozens of sex-linked scaffolds on our three focal species (supplementary results, Supplementary Material online and fig. 2). For each focal species, we performed simulations to evaluate the probability of observing the number of sex-linked SNPs by chance. For H. boettgeri and P. parva, but not X. mellotropicalis, the simulations indicate that a number of identified sex-linked SNPs are highly unlikely to occur by chance (supplementary fig. S1, Supplementary Material online). Thus, even though the genome assemblies for each of our focal species were highly fragmented and the family sizes of each was relatively small, the generally high synteny among frog genomes (Sun et al. 2015; Session et al. 2016) and close phylogenetic affinities of pipid genomes allowed us to demarcate large regions of sex-linked recombination suppression using RRGS and WGS data in two pipid species. For X. mellotropicalis, our data indicate with high confidence that the sex-specific region of the sex chromosomes of this species is small; a consequence of this finding is that the support for female heterogamy is relatively weak in terms of the number of sex-linked SNPs in the RRGS data (one completely sex-linked SNP on an unplaced scaffold and one almost sex-linked SNP on chromosome 7; supplementary information, table S2, and fig. S1, Supplementary Material online).
To further validate these inferences based on RRGS data mapped to WGS assemblies, we Sanger sequenced genes from our laboratory families that we inferred to be within the sex-linked regions of each species based on the RRGS and WGS data. These amplifications targeted regions that included sex-linked variation that was spanned by the RRGS data, and also regions that were not covered by the RRGS data but that were nonetheless putatively sex-linked in each focal species under the assumption of synteny with the X. tropicalis genome. The Sanger sequences verified complete sex-linkage for H. boettgeri (four independent amplifications of portions of three genes) and P. parva (three independent amplifications of portions of three genes). Sanger sequencing verification was unsuccessful for a sex-linked scaffold identified by the RRGS data for X. mellotropicalis due to repetitive sequences in the very small sex-linked region that we identified. However, a paralog of one gene (or8h1) was identified based on its genomic position in the X. tropicalis genome that had an almost female-specific sex-specific amplification. This suggests that or8h1 is in close genomic proximity to the trigger for sex determination in X. mellotropicalis, but that recombination does occasionally occur between these two genes (supplementary fig. S2 and results, Supplementary Material online). This locus and another sex-linked SNP from X. mellotropicalis map to the short arm of chromosome 7, which is also where the master sex-determining gene of X. tropicalis is thought to reside (fig. 1; Olmstead et al. 2010; Roco et al. 2015). The homologous genomic locations of the trigger for sex determination in the closely related species X. tropicalis and X. mellotropicalis provides further validation for our inference of the identity of the sex chromosomes of X. mellotropicalis, even though this region is very small in both species.
Discussion
Developmental Systems Drift of Sex Determination and Recombination Suppression
Developmental processes are orchestrated by networks of genetic interactions. When compensatory changes occur in network components (i.e., genes) (Johnson and Porter 2006; Lynch and Hagner 2015), these networks can diverge in different species with minimal phenotypic consequence; this phenomenon is called DSD (True and Haag 2001). Several lines of evidence—some newly reported here—demonstrate that DSD of sex determination occurred in pipid frogs. A gene called dm-w that reliably triggers female differentiation in X. laevis is not present in some species. In other species, this gene is genetically associated with female differentiation with imperfect efficiency (i.e., usually—but not always—because it is sometimes present in males). And, in some species, we failed to detect dm-w in some females, suggesting that this gene is not required for female differentiation in these species. Variation among pipid species in male versus female heterogamy also evidences DSD because it indicates variation in the dosages of sex-specific factors that trigger sexual differentiation.
Our findings concerning the presence, absence, and sex-specificity of a sex-determination gene (dm-w) identify several mechanisms by which DSD can occur. In some species this gene is found in the majority of phenotypic females and also in a minority of phenotypic males. Interpretation of these results in a phylogenetic context (fig. 1) argues that the ancestral capacity of dm-w to trigger female differentiation was strong but incomplete. By extension, this also argues that when dm-w initially began to contribute to the genetic control of sex determination, this gene was usually present as a single allele in females, because most of their fathers would have lacked dm-w. Dm-w is now strictly coupled to female differentiation in two closely related species (X. laevis and X. gilli), which indicates DSD occurred via an empowerment of the sex-determination capacity of dm-w with respect to its ancestral condition.
Although dm-w originated in an ancestor of some extant Xenopus, in most (four of six) of the descendant species of this ancestor that we surveyed, this gene was not female-specific. In two of these species (X. itombwensis and X. victorianus), two alleles of dm-w were detected in both sexes (heterozygous sites are present in both sexes). Furthermore, the genomic location of dm-w in X. itombwensis may be on an autosome or pseudoautosomal region (because dm-w was present in all individuals of both sexes; table 1 and supplementary table S4, Supplementary Material online). That X. itombwensis is derived from an ancestor in which dm-w was usually segregating as a single allele in females, suggests that DSD also occurred via mechanistic sidelining of dm-w (i.e., this gene is still present and possibly functional, but not as in a sex-specific or sex-biased capacity as a trigger for sex determination). We did not find evidence in the coding regions of exons 2 and 3 that this sidelining rendered dm-w nonfunctional in X. itombwensis. Rather, its role in sexual differentiation may instead operate in a downstream capacity relative to some other, as yet unidentified, trigger for sexual differentiation. In several species, we did not observe dm-w in the targeted next-generation sequencing of one female. Whether and how many times dm-w was lost depends on how many of these species also do not carry this allele in other individuals. However, several species in which one female lacks dm-w are closely related; this is consistent with the possibility that this gene was lost in their collective most recent common ancestor (fig. 1)—a scenario which is consistent with DSD by gene deletion. It is possible some of these observations are a consequence of intraspecific regulatory variation in dm-w (e.g., dm-w alleles with high or low expression may tend to occur in females or males, respectively), and/or intraspecific variation in molecular variation in other portions of dm-w not surveyed here (e.g., coding regions of exon 4) or other genes (e.g., scan-w in X. laevis; Mawaribuchi, Takahashi, et al. 2017).
Overall, this information combined with other reports cited below evidence at least seven distinct systems that trigger sex determination in pipid frogs, with each being distinguished by genomic location (expressed below in terms of orthology to chromosomes of X. tropicalis), independent origin, and heterogamy. These systems include:
The male heterogametic system on chromosome 4 of H. boettgeri (this study).
The male heterogametic system on chromosome 6 of P. parva (this study).
The putative female heterogametic system on chromosome 7 in subgenus Silurana of X. mellotropicalis (this study), some populations of X. tropicalis (Roco et al. 2015), and possibly close relatives.
The male heterogametic system on chromosome 7 of X. tropicalis (Roco et al. 2015), which the principle of maximum parsimony suggests evolved from a female heterogametic ancestor (this study).
The female heterogametic system on chromosome 2 in subgenus Xenopus of species where dm-w is female-specific (e.g., X. laevis; Yoshimoto et al. 2008; Session et al. 2016) or usually in females (e.g., X. clivii; this study).
The female heterogametic system on chromosome 8 of X. borealis, and possibly close relatives whose ancestor lost dm-w (Furman and Evans 2016).
The non-dm-w based systems on unknown chromosome(s) of other species in subgenus Xenopus that might have lost dm-w completely (e.g., X. longipes), or where it segregates autosomally (X. itombwensis) (this study; table 1).
In species that do not vary in the heterogametic sex but do differ in the genomic location of the trigger for sex determination (e.g., P. parva and H. boettgeri), DSD may have occurred if there exists a novel trigger for sex determination, or alternatively translocation of a trigger may have occurred without DSD. We did not attempt to collect information about whether translocations account for the variation in the genomic locations of triggers for sex determination because, apart from dm-w, these triggers have not been identified.
Within pipid frogs, a roughly bimodal pattern is evident in which recombination suppression on the sex chromosomes is either restricted to a very small region (less than a few Mb) or alternatively spans a large region (greater than ∼10 Mb). Previous work on the allotetraploid species X. laevis and X. borealis, for example, illustrates that recombination suppression affects only a small region linked to dm-w on chromosome 2 L in the former species, but a large region spanning ∼50 Mb that carries an unidentified sex-determining factor on chromosome 8 L in the latter (Furman and Evans 2018). Similar to X. laevis, X. tropicalis also has a very small region of suppressed recombination (Bewick et al. 2013; Mitros et al. 2019). Here, we demonstrate that the pattern in the X. mellotropicalis resembles that of X. laevis and X. tropicalis, but the patterns in H. boettgeri and P. parva resemble X. borealis. In P. parva, we estimated the size of the region of suppressed recombination to minimally span multiple Mb (fig. 2); there was a dearth of high-confidence genotypes on one side, so the size of this region may be even larger. The sex chromosomes of pipid frogs thus demonstrate that homomorphic sex chromosomes might have large regions of suppressed recombination; a similar situation is found in ratite birds (Vicoso et al. 2013; Yazdi and Ellegren 2014).
One speculation that could account for these species-level differences is that there were (or are) more loci with sexually antagonistic effects on the chromosomes with more recombination suppression (chr 4, 6, 8) than those with less (chr 2, 7) (Rice 1987). However, another plausible explanation for these differences in the extent of recombination suppression instead relates to the efficacy of the trigger for sex determination in each species; we discuss this possibility next.
Inefficient Sex Determination as a Mechanism for Nondivergence of Sex Chromosomes
These findings of DSD of triggers for sex determination have the potential to affect genome evolution in many ways, including the evolution and recombination patterns of entire sex chromosomes. For example, there are several ways that the function and age of the trigger for sex determination could impinge on the genetic fate of entire sex chromosomes—including whether they are heteromorphic or homomorphic and whether a large or small region of recombination suppression is present. The fountain of youth hypothesis, for example, predicts that nondiverged sex chromosomes could be associated with phenotypic sex reversal if the heterogametic sex has low rates of recombination (Perrin 2009). However, in X. laevis and X. borealis, recombination occurs more frequently in oogenesis than spermatogenesis, including on the pseudoautosomal regions of the sex chromosomes, even though females are heterogametic in both of these species (Furman and Evans 2018). In Xenopus, recombination on the ends of chromosomes is lower in females than males (Furman and Evans 2018). Thus, because dm-w resides on the tip of chromosome 2 L, recombination is less likely to disrupt this gene when it is found in females.
The rapid turnover hypothesis posits that new sex-determination loci arise in different locations before recombination suppression results in sex chromosome heteromorphy (Volff et al. 2007). Our findings of at least seven distinct sex-determination systems in pipid frogs are generally consistent with the expectations of the rapid turnover hypothesis, but suggest other mechanisms are at play as well. For example, the sex chromosomes of X. clivii and X. laevis—species whose ancestors diverged millions of years ago (fig. 1; Evans et al. 2015)—are homomorphic (Tymowska 1991), even though they both have dm-w and presumably orthologous sex chromosomes. Another example is the sex chromosomes of X. borealis, which have a recently derived sex-determining system, but widespread recombination suppression (Furman and Evans 2018), indicating that recent turnovers are not necessarily associated with small regions of recombination suppression.
Another possibility is that the efficiency of the trigger for sex determination could influence the evolutionary fate of entire sex chromosomes. The potential to resolve genomic conflict associated with sexual antagonism via linkage to an inefficient trigger for sex determination (i.e., that is not sex-specific) is presumably modest compared with linkage to an efficient trigger (that is sex-specific). If resolution of sexual antagonism favors recombination suppression in sex-linked regions, it follows that a small region of suppressed recombination would be expected near an inefficient trigger for sex determination, even in the absence of rapid turnover. Indeed, our findings suggest that the ancestral sex-determining locus of X. laevis used to be inefficient, and the region of recombination suppression in this species is known to be small (Furman and Evans 2018). Likewise, the sex-determining gene in X. mellotropicalis is possibly homologous to that of some strains of X. tropicalis, and the most recent common ancestor of these two species might have had a polymorphic and/or inefficient trigger of sexual differentiation as X. tropicalis does now (Roco et al. 2015; Mitros et al. 2019). This inefficient/polymorphic sex-determining locus could explain why both of these species have small nonrecombining regions on their sex chromosomes (fig. 2; Bewick et al. 2013). Conversely, species with large regions of suppressed recombination, such as H. boettgeri, P. parva, and X. borealis (Furman and Evans 2018) might have triggers of sex determination that operate with high fidelity. In these species, expansion of recombination suppression in genomic regions that flank efficient sex-determining loci may have been favored by natural selection because it resolved genomic conflict associated with mutations with sexually antagonistic fitness effects.
Other examples of inefficient triggers for sex determination have been reported. In the frog Rana temporaria—a species with homomorphic sex chromosomes—dmrt1 is a candidate trigger for sex determination, but shows incomplete linkage to the sex phenotype (Rodrigues et al. 2017). This also could arise if this trigger were inefficient in some contexts (or populations), and/or if there is a polygenic basis of sex determination. Although the mechanism of sex determination differs from frogs (Hediger et al. 2004), a similar genomic situation is present in the housefly Musca domestica, where an inefficient sex-determining locus on a Y chromosome is linked to male differentiation, with the extent of linkage depending on input from multiple sex-determining alleles on other chromosome pairs (Denholm et al. 1986; Tomita and Wada 1989). These examples of inefficient triggers for sex determination, including dm-w in some Xenopus, allow for a novel explanation for why sex chromosomes may stay homomorphic; this explanation does not require or preclude rapid turnover of sex chromosomes, or sex differences in the rate of recombination.
What Factors Govern the Efficiency of dm-w?
In X. laevis, dm-w is thought to drive female differentiation via competitive inhibition of the male-related gene dmrt1 (Yoshimoto et al. 2008, 2010) and expression of one dmrt1 paralog (dmrt1-S) is higher in male than female tadpole gonads (Mawaribuchi, Musashijima, et al. 2017). Related to this, the mechanism by which allopolyploidization is thought to occur in Xenopus does not involve duplication of the W-chromosome (or dm-w), but does involve duplication of autosomal genes such as dmrt1 (Kobel and Du Pasquier 1986). Thus, a single dose of dm-w may compete with variable dosages of dmrt1 in different individuals or species, depending on ploidy level and pseudogenization. Interestingly, pseudogenization of dmrt1 homologs occurred independently several times in Xenopus polyploids (Bewick et al. 2011), and it is conceivable that population-level variation in pseudogenization if dmrt1 and/or other loci affects the efficacy of dm-w to determine sex.
Outlook
The dynamic mechanistic capacity of dm-w reported here demonstrates that DSD can arise from altered function or interactions of a core component of a developmental system. This could stem from changed pleiotropic interactions (Pavlicev and Wagner 2012), changed stoichiometric relationships as a result of gene duplication and pseudogenization, subfunctionalization (Force et al. 1999), and neofunctionalization (Lynch et al. 2001). Related to this, empirical and theoretical observations predict that some chromosomes or genomic regions may be more likely than others to be recruited as triggers for sex determination (O’Meally et al. 2012; Brelsford et al. 2013; Blaser et al. 2014; Furman and Evans 2016). The possibility that certain genomic regions are repeatedly co-opted as sex chromosomes is consistent with the mechanism observed here that DSD involves functional or interaction changes of components of conserved genetic pathways, as opposed to convergence in different species of pathways with distinctive genetic functions (Haag and Doty 2005). Indeed, the sex chromosomes we identified here harbor several sex-related genes (supplementary table S7, Supplementary Material online).
In various ranid and bufonid frogs, orthologs of X. tropicalis chromosomes 1, 2, 3, 5, and 8 are sex chromosomes (Miura 2007; Brelsford et al. 2013; Uno et al. 2015; Jeffries et al. 2018). The diversity of sex chromosomes in pipids—including chromosomes homologous to five of the ten X. tropicalis chromosomes (2, 4, 6, 7, or 8)—expands this list, leaving homologs of only X. tropicalis chromosomes 9 and 10 as not (yet) being identified as sex chromosomes in a frog species. However, several chromosome pairs have been repeatedly co-opted as anuran sex chromosomes (e.g., homologs of X. tropicalis chromosome 5 in ranids, Jeffries et al. 2018; chromosome 1 in bufonids and ranids, Brelsford et al. 2013; chromosome 8 in ranids and pipids, Miura 2007; Uno et al. 2015; Furman and Evans 2016; and chromosome 7 in pipids and rhacophorids, Uno et al. 2015). The newly identified sex chromosome 6 of P. parva is also homologous to the sex chromosomes of some snakes (Matsubara et al. 2006). Clearly, more information on the identities of genes that trigger sexual differentiation in species that lack dm-w would allow us to evaluate the degree to which DSD generally involves modulation of core developmental genes, or alternatively is largely sculpted by inputs from mechanistically divergent systems (Cline et al. 2010; Zhang et al. 2016).
Overall, this study recovers evidence for frequent DSD of sex determination in pipid frogs, identifies new sex chromosomes and novel variation in recombination suppression on sex chromosomes, and pinpoints three of mechanisms of DSD including loss, sidelining, and empowerment of a mechanistically influential gene. We speculate that the sex chromosomes of species with inefficient sex-determination genes are less advantageous destinations for alleles with sexually antagonistic function. This points to a previously unsuspected mechanism for sex chromosome nondivergence (homomorphy) that will be fascinating to further explore in other groups with known triggers for sex determination.
Materials and Methods
Targeted Next-Generation Sequencing and Sanger Sequencing of dm-w in Xenopus
In order to survey Xenopus females for the presence of dm-w, we used oligonucleotide baits to enrich Illumina libraries for sequences similar to dmrt1, which is the gene from which exons 2 and 3 of dm-w originated (Yoshimoto et al. 2008). We captured sequences from libraries that were generated from one female from all described species in subgenus Xenopus, except X. fraseri, because we lacked a sample of high-quality genomic DNA from this species, and X. lenduensis and X. mellotropicalis where we generated libraries from male individuals. As a negative control, we captured sequences from libraries of one individual from each species in subgenus Silurana (this lineage diverged prior to the origin of dm-w). In order to evaluate female-specificity of dm-w, we used PCR and Sanger sequencing dm-w in several individuals from six species (table 1) whose phenotypic sex was determined from external morphology, or surgically after euthanasia. Additional details of samples, PCR primers, library preparation, capture, sequencing, assembly, and analysis are provided in Supplementary Material online. Capture sequences of dm-w and dmrt1 are available on GenBank (accession numbers MN030659–MN030916, supplementary table S5, Supplementary Material online) and an alignment of these data is available in supplementary data, Supplementary Material online (supplementary file S2 for exon 2 and file S3 for exon 3, Supplementary Material online).
The Sex Chromosomes of Other Pipids
The sex chromosomes of three pipid species that diverged before the origin of dm-w were also examined, including the dwarf clawed frog, H. boettgeri, the Sabana Surinam toad, P. parva, and the Gabonese clawed frog, X. mellotropicalis. The first two of these species are diploid (H. boettgeri: 2n = 24 and P. parva: 2n = 30; Cannatella and De Sa 1993) but the third (X. mellotropicalis) is allotetraploid (2n = 4× = 40; Evans et al. 2015). Half of the genome of X. mellotropicalis (i.e., one subgenome) is derived from a recent ancestor of the diploid species X. tropicalis, and its other subgenome is derived from another diploid species that is either extinct or undiscovered (Evans et al. 2015).
We generated one family from each of these three pipid species, and performed WGS on one or both parents. Assembled WGS reads were used to anchor RRGS and Sanger sequencing from both parents and their offspring to identify sex-linked genomic regions. The sex of postmetamorphic offspring was ascertained surgically after euthanasia. The H. boettgeri family included both parents, 11 daughters, and 21 sons. The P. parva family included both parents, five daughters, and seven sons. The X. mellotropicalis family included both parents, nine daughters, and nine sons. Additional details about these families and analysis are provided in Supplementary Material online.
Acknowledgments
We are especially grateful to Svante Pääbo and Max Planck Institute for Evolutionary Anthropology for hosting B.J.E. on a sabbatical in Leipzig, Germany. We thank Darcy Kelley for providing us with X. mellotropicalis and X. clivii individuals, Ben Vernot for assistance with capture probe design, Birgit Nickel for performing library capture, Brian Golding and the Shared Hierarchical Academic Research Computing Network (SHARCNET: www.sharcnet.ca) and Compute/Calcul Canada for access to computational resources, and Brian Golding, Jonathan Dushoff, Ian Dworkin, and four anonymous reviewers for enlightening discussions and comments on a draft of this article. This work was supported by the Natural Science and Engineering Research Council of Canada (RGPIN/283102-2012 and RGPIN-2017-05770 to B.J.E.), the Museum of Comparative Zoology (B.J.E.), the Max Planck Institute for Evolutionary Anthropology (M.-T.G., M.M., B.J.E.), a National Geographic Research and Exploration Grant (no. 8556-08 to E.G.), and the US National Science Foundation (DEB-1145459 to E.G.).
