EVOLUTION OF THE CLOSELY RELATED, SEX‐RELATED GENES DM‐W AND DMRT1 IN AFRICAN CLAWED FROGS (XENOPUS)

Abstract

DM‐W is a dominant, female‐specific, regulator of sex determination in the African clawed frog Xenopus laevis. This gene is derived from partial duplication of DMRT1, a male‐related autosomal gene. We set out to better understand sex determination in Xenopus by studying this pair of genes. We found that DM‐W evolved in Xenopus after divergence from the sister genus Silurana but before divergence of X. laevis and X. clivii, and that DM‐W arose from partial duplication of DMRT1β, which is one of the two DMRT1 paralogs in the tetraploid ancestor of Xenopus. Using the rate ratio of nonsynonymous to synonymous substitutions per site and multilocus polymorphism data, we show that DM‐W evolved non‐neutrally. By cloning paralogs and using a pyrosequencing assay, we also demonstrate that DMRT1 underwent phylogenetically biased pseudogenization after polyploidization, and that expression of this gene is regulated by mechanisms that vary through development. One explanation for these observations is that the expression domain of DMRT1β was marginalized, which would explain why this paralog is dispensable in Xenopus polyploids and why DM‐W has a narrow expression domain. These findings illustrate how evolution of the genetic control of stable phenotypes is facilitated by redundancy, degeneration, and compartmentalized regulation.

Expression divergence, gene duplication, positive selection, sex determination

Sex determination in frogs is genetically triggered (Hayes 1998) but is achieved via variable mechanisms, making this group a compelling but understudied model for studying evolution of this crucial phenotype. Recently, a female‐specific, W‐linked gene called DM‐W was discovered in the African clawed frog Xenopus laevis (Yoshimoto et al. 2008), supporting earlier conclusions that X. laevis females are heterogametic (Mikamo and Witschi 1963). DM‐W is the first master regulator of primary sex determination to be identified in amphibians. The 5′‐coding region of DM‐W is very similar to the 5′‐coding region of the doublesex and mab‐3 related transcription factor (DMRT1) gene, suggesting an origin of DM‐W by partial gene duplication of DMRT1 (Yoshimoto et al. 2008). This is also suggested by the exon structure of these genes, which are similar in the 5′ portion but not the 3′ portion (Fig. 1A). The homologous region of DMRT1 and DM‐W contains a DNA‐binding motif called a “DM domain” (Yoshimoto et al. 2008) and both genes bind to shared DNA sequences (Yoshimoto et al. 2010), suggesting that they are transcription factors that interact with the same regulatory elements. Analysis using full‐length sequences suggested DM‐W arose after the diversification of ancestors of frogs and fish but before divergence of the ancestors of frogs, birds, and mammals (Supplementary information of Yoshimoto et al. 2008) but the genome sequence of the Nigerian clawed frog (Silurana tropicalis) does not contain DM‐W (Yoshimoto et al. 2008), even though the individual sequenced was female (Hellsten et al. 2010).

1

Regions of homology between DM‐W and DMRT1 and proposed competitive binding during sexual differentiation of females. (A) Exons (boxes) and introns (lines) of X. laevis DM‐W and DMRT1 with homologous region used for analysis of ω in black. DM‐W exon structure follows Yoshimoto et al. (2008) and DMRT1 exon structure follows Osawa et al. (2005) and is based on S. tropicalis. Numbers below each exon reflect how many amino acids are encoded; the locations of start (arrow) and stop codons are indicated. (B) DM‐W (black) and DMRT1 (gray) bind to the same cis‐ regulatory factors indicated by black squares (Yoshimoto et al. 2010). In males, DMRT1 drives expression of male‐related genes. In females, DM‐W is expressed slightly earlier during primary gonadal differentiation and is thought to repress transcription of male‐related genes (Yoshimoto et al. 2008).

DM domain‐containing genes play a remarkably conserved role as an activator of male differentiation in metazoans including worms, flies, coral, birds, and humans (Burtis and Baker 1989; Raymond et al. 1999; Yi and Zarkower 1999; Raymond et al. 2000; Miller et al. 2003; Haag and Doty 2005). DMRT1 is broadly expressed during development of X. laevis and has been detected in unfertilized eggs, early tadpole development of both sexes, in the gonads of both sexes during primary gonadal differentiation, and in postmetamorphic testis and ovary, with expression becoming increasingly male‐biased in testis compared to ovary by 1–5 months after metamorphosis (Osawa et al. 2005; Yoshimoto et al. 2006; Yoshimoto et al. 2008). No sex difference in DMRT1 expression was detectable during gonadal differentiation (Yoshimoto et al. 2008). In contrast, DM‐W appears to be expressed only in female gonads during primary gonadal differentiation, with its peak expression level just prior to a surge of DMRT1 expression in the gonads of both sexes (Yoshimoto et al. 2008). In females, DM‐W expression may antagonize DMRT1‐activation of male‐specific genes by binding to and inhibiting regulatory regions recognized by both proteins (Fig. 1B; Yoshimoto et al. 2008, 2010).

African clawed frogs include the genera Silurana and Xenopus and a diversity of species that evolved through allopolyploidization (Evans 2008). Species of Xenopus have 36, 72, or 108 chromosomes and are tetraploid, octoploid, or dodecaploid, respectively. Xenopus laevis is tetraploid but considered “pseudotetraploid” because of disomic inheritance (Tymowska 1991). Species of Silurana have 20 or 40 chromosomes and are diploid (such as S. tropicalis), or tetraploid. Because of the unique way that allopolyloidization occurs in Xenopus (Kobel 1985, 1996a), DM‐W and linked female‐specific genes on the nonrecombining portion of the W chromosome (hereafter NRW) are not duplicated by allopolyploidization (Fig. 2). In contrast, autosomal genes such as DMRT1 (Uno et al. 2008; Yoshimoto et al. 2008) are duplicated by allopolyploidization. Thus females of tetraploid, octoploid, and dodecaploid species are expected to have only one allele of DM‐W but up to four, eight, or 12 DMRT1 alleles organized in two, four, or six paralogous loci, respectively.

2

Autosomal genes (A) such as DMRT1 are duplicated by allopolyploidization in clawed frogs but DM‐W (W) is not. Whole‐genome duplication (WGD) proceeds via F1 hybrid and backcrossed females that produce unreduced eggs. The mechanisms of sex determination in the triploid and tetraploid progeny is unclear, although both sexes do exist in laboratory‐generated polyploids. Presumably diploidization causes two Z chromosomes (Z) to become autosomes (A_Z). This figure is adapted from Kobel and Du Pasquier (1986) and Kobel (1996b).

How evolutionarily conserved is the “DM‐W/DMRT1 system” of sex determination in frogs and how did this system arise? To what extent and in what way are the components of this crucial system variable across species and ploidy levels, and how is this influenced by natural selection? To begin to address these questions, we examined the molecular evolutionary history of DM‐W and DMRT1. Using a PCR assay, we screened frog species closely related to X. laevis for evidence of DM‐W and we analyzed phylogenetic relationships between homologous portions of DM‐W and cloned paralogs of DMRT1. Using these data and polymorphism data from natural populations, we tested for evidence of nonneutral evolution of DMW and linked regions of the NRW. We also examined patterns of pseudogenization in paralogs of DMRT1 in multiple polyploid species of Xenopus. Additionally, with an aim of better understanding regulatory evolution of this system, we used a pyrosequening assay to characterize the mechanisms of expression divergence of one DMRT1 paralog at two crucial stages of development: (1) during primary sexual differentiation in tadpoles and (2) in adult testis.

Methods

GENETIC SAMPLES, PHYLOGENETIC ESTIMATION, CODON ANALYSIS

We analyzed evolution and expression of DM‐W and DMRT1 using genetic samples obtained from a Xenopus colony that was in Geneva, field collections by Ben J. Evans, and animals from a commercial supplier (Xenopus Express, Brooksville, FL). Samples from a natural populations of X. laevis and X. gilli were obtained from Cape Province, South Africa in ponds near Betty's Bay as detailed elsewhere (Evans et al. 1997, 1998) and samples from the population in the Democratic Republic of the Congo (DRC) were collected by Ben J. Evans in the town of Lwiro.

Sequences from DMRT1 paralogs and DM‐W were obtained using primers detailed in Supporting information. For DMRT1, paralogs were co‐amplified and cloned using the TOPO TA cloning kit (Invitrogen). Estimation of evolutionary relationships among homologous portions of DMRT1 paralogs and DM‐W was performed using Bayesian phylogenetics as implemented by MrBayes version 3.2 (Huelsenbeck and Ronquist 2001) using a model selected by the Akaike information criterion (AIC) with MrModeltest (Nylander 2004). Dirichlet priors were used for nucleotide frequencies, two independent Markov Chain Monte Carlo runs were performed, each for 1,000,000 generations. A DMRT1 sequence from the toad Bufo marinus (Genbank accession number FJ697175) was used as an outgroup. Convergence was assessed using plots of likelihood values and parameter estimates from the posterior distribution using R (R Development Core Team 2005) and based on this, a burn‐in of 100,000 generations was discarded.

Analysis of the dN/dS ratio (ω) of the homologous region of DM‐W and DMRT1 was performed with the codeml module of PAML version 4.2 (Yang 1997) using models described in further detail below. The coding sequence analyzed included 98 codons. It did not include the first 17 codons of exon 2 because our forward PCR primer annealed in this region and it did not include the last six codons of exon 3 because this region was deemed nonhomologous. We included DMRT1 sequences amplified from cDNA so there was not missing data surrounding the junction of exons 2 and 3 in this analysis.

To explore whether changes within the DM domain of DM‐W are unique to this gene (and therefore potentially functionally significant), we compared this 68 amino acid long region to the DM domain of the DMRT1 gene in other lineages including DMRT1 paralogs of six fish (Gasterosteus aculeatus AY867870, Paralichthys olivaceus EU490514, Takifugu rubripes NM_001037949, Xiphophorus maculates AF29187, Epinephelus merra EU555179, Oryzias latipes AY442916), chicken (Gallus gallus AF123456), a toad (B. marinus FJ697175), a crocodile (Crocodylus palustris EU531727), and various mammals (Homo sapiens NC_000009, Pan troglodytes XM_528528, Mus musculus NM_015826, Bos taurus XM_002689628, Oryctolagus cuniculus XM_002708188, Canis familiaris XM_846402, and Sus scrofa NM_214111). We also compared the DM domain of DM‐W to the DM domain of the DMY gene of O. latipes (NM_001104680), which is a male‐specific sex determining gene, in the context of each genes’ closest paralog.

MULTILOCUS POLYMORPHISM DATA

By analyzing expressed sequence tag databases with the assistance of Frédéric Chain, we identified genes in X. laevis that we suspected to be singletons due to loss of one paralog after whole‐genome duplication. To test for evidence of recent nonneutral evolution, we sequenced these putative single‐copy genes and also portions of DM‐W and linked flanking regions of the NRW. After accounting for differences in mutation rate and the number of silent sites sequenced, our null expectation is that the level of polymorphism at the NRW should be one‐fourth that of the other loci as a consequence of differences in mode of inheritance and copy number. The likelihood of neutral evolution, where the effective population size of the NRW is one‐fourth that of autosomal loci (i.e., the inheritance scalar = 0.25), and nonneutral evolution (the inheritance scalar < 0.25), was estimated using a maximum likelihood version of Hudson, Kreitman, Aguade test (Hudson et al. 1987) with mlHKA version 2 (Wright and Charlesworth 2004). Divergence from an outgroup was used to accommodate differences in rates of evolution among loci. Xenopus gilli and X. laevis are sister species and sequences from each were used as an outgroup for the other. PCR primers are provided in Supporting information.

PHYLOGENETICALLY BIASED PSEUDOGENIZATION

To better understand pseudogenization after genome duplication in DMRT1, we sequenced 367 clones amplified from DMRT1 exon 2 and 155 clones amplified from DMRT1 exon 3 from genomic DNA of a total of 18 species including tetraploids: X. laevis, X. borealis, X. muelleri, X. pygmaeus, X. fraseri, X. largeni, X. clivii, and an undescribed species X. new tetraploid 1 (Evans et al. 2004), octoploids: X. itombwensis, X. boumbaensis, X. vestitus, X. andrei, X. wittei, X. ameiti, and an undescribed octoploid X. new octoploid 2 (Evans et al. 2010), and dodecaploids: X. longipes, X. ruwenzoriensis, and an undescribed dodecaploid X. cf. boumbaensis (Evans 2007). We used a model‐based approach implemented by the program BayesTraits (Pagel and Meade 2006) to test for phylogenetic bias in the pattern of pseudogenization in this gene. We compared the likelihood of a model with one rate of pseudogenization in both major paralogous lineages of DMRT1 (α and β) to the likelihood of a model with a separate rate of pseudogenization in each lineage. The rate of change of pseudogenes to functional paralogs was fixed to zero under the assumption that once a paralog became a pseudogene, functional resuscitation is not possible. A phylogeny for this analysis was estimated using up to 3461 base pairs of sequence from cloned paralogs of the recombination activating genes 1 and 2, as detailed in Evans (2007). Following model selection and analytical methods of Evans (2007), phylogenetic relationships were estimated with the data partitioned by codon position using MrBayes (Huelsenbeck and Ronquist 2001). Two independent runs were performed, each for 2,000,000 generations, and 100,000 generations discarded as burn‐in. For the BayesTraits analysis, a chronogram was constructed from the consensus topology recovered from the Bayesian analysis using penalized likelihood as implemented by r8s version 1.71, with the smoothing parameter obtained from a cross‐validation procedure (Sanderson 1997, 2002).

EXPRESSION DIVERGENCE

Gene regulation is orchestrated by cis‐acting factors that independently affect transcription of each allele and by trans‐acting factors that affect expression of both alleles. To quantify the degree to which cis‐ and trans‐ acting factors drive DMRT1expression divergence between species, we used an experimental approach that compares the expression ratios of species‐specific alleles in each species and in their F1 hybrid (Wittkopp et al. 2004, 2008; Landry et al. 2005). In hybrids, trans‐acting factors from both parental species interact with the cis‐acting factors of species‐specific alleles, so divergent expression of species‐specific alleles can be attributed exclusively to cis‐ divergence. In the parental species, expression divergence is the culmination of cis‐ and trans‐ factors. Comparison of the expression ratio of species‐specific alleles in hybrids to the parental species therefore allows one to dissect apart and quantify cis‐ and trans‐ contributions to expression divergence between species. For example, trans‐ only divergence is inferred if species‐specific alleles are expressed at similar levels in hybrid individuals but at different levels in the parental species. Note that this approach compares expression ratios of species‐specific alleles (which are identified by species‐specific nucleotide polymorphisms), and does not compare the magnitude of expression. This is because the magnitude of expression between parental individuals and hybrids could vary as a consequence of upstream factors that impact entire genetic networks—even in the absence of expression divergence between species—whereas the expression ratio in parental individuals and hybrids should be affected by mechanisms that operate at the level of an individual locus (Wittkopp et al. 2004). We performed this assay using X. laevis and X. borealis on tissue from adult testis and tissue from tadpole stage 53, which is after the gonads develop but before they differentiate in X. laevis.

The ratio of expression of species‐specific alleles of DMRT1α in the parental species and their F1 hybrid was quantified using a Biotage PSQ96 pyrosequencer based on species‐specific and paralog‐specific single nucleotide polymorphisms in exon 3 (Supporting Information). Expression ratios were estimated in F1 hybrids generated from a cross between an X. laevis female and an X. borealis male (six tadpoles and two adults) or the reciprocal cross (four tadpoles and two adults), and compared to ratios the corresponding expression ratios in mixtures of tissues from each parental species (seven tadpole and four adult parental mixes). According to the manufacturer's protocol, genomic DNA is efficiently removed from the RNA extraction but we nonetheless implemented an optional DNAse digestion step in this procedure. Expression ratios in parental mixes were normalized by the amount of genomic DNA in each mix (which was extracted from an aliquot taken before the DNAse digestion step) following Landry et al. (2005). For genomic DNA normalization we used data from DMRT1α, one paralog of the recombination activating gene 2 (RAG2), and both paralogs of the recombination activating gene 1 (RAG1). Pyrosequencing primers and normalized expression ratios are provided in Supplemental Information and elsewhere (Anderson and Evans 2009).

Statistical analysis of mechanisms of expression divergence followed Landry et al. (2005) and used the “proc mixed” implementation of restricted maximum likelihood of SAS version 9.1.3 (SAS Institute) with modified scripts provided by Patricia Wittkopp. Student's t‐tests were computed within the mixed procedure and locus‐level significance was interpreted after sequential Bonferroni correction for two tests (Rice 1989). Application of the Bonferroni correction makes the detection of antagonistic cis‐ and trans‐ divergence more conservative (Anderson and Evans 2009). Conclusions were identical to those drawn from analysis using a standard Student's t‐test.

Results

DM‐W AROSE AFTER DIVERGENCE OF XENOPUS AND SILURANA BUT BEFORE DIVERGENCE OF X. LAEVIS AND X. CLIVII

The DM‐W gene was first identified in X. laevis but a search of the genome of S. tropicalis failed to recover this gene (Yoshimoto et al. 2008). To investigate what other species have DM‐W, we attempted to amplify and sequence portions of DM‐W orthologs in species that are more closely related to X. laevis than to S. tropicalis. We found DM‐W in X. andrei, X. clivii, X. gilli, X. itombwensis, X. laevis, X. largeni, X. pygmaeus, and X. vestitus (Genbank accession numbers HQ220848–54; HQ220877–8). These species are tetraploid except X. itombwensis and X. vestitus, which are octoploid.

Based on a previously published phylogeny of clawed frogs (Fig. 3), the most distantly related species from X. laevis in which we detected a DM‐W ortholog was X. clivii. The phylogenetic position of X. clivii is unclear. Mitochondrial DNA strongly supports a sister relationship between X. clivii and Clade 1 in Figure 3 (Evans et al. 2004) but two tightly linked nuclear loci (RAG1 and RAG2) strongly support a sister relationship between X. clivii and Clade 2 in Figure 3 (Evans 2007).

3

DM‐W was detected in a subset of Xenopus species whose most recent common ancestor and descendants are indicated in gray. A “D” next to species names indicates species from which portions of DM‐W was amplified. Nodes with 1 or 2 inside refer to Clades 1 and 2, respectively. An asterisk follows species names for which amplification of DM‐W was not attempted. This phylogeny represents reticulating relationships among species that stem from allopolyploidization; allopolyploid speciation is indicated by asterisks on internal nodes and the ploidy level of each species is indicated after the species name. Maternal or biparental relationships are indicated with solid lineages and paternal contributions to allopolyploid speciation events are indicated with dashed lineages. Resolved nodes have >95% posterior probability and are based on the linked autosomal genes RAG1 and RAG2 and mitochondrial DNA; details of phylogenetic estimation are given in Evans et al. (2004, 2005) and Evans (2007). Strongly supported but conflicting relationships were recovered for X. clivii with respect to numbered Clades 1 and 2 (see text) so this relationship is represented as a polytomy.

We were not able to amplify DM‐W from the dodecaploid species X. longipes and X. ruwenzoriensis or the tetraploid species X. muelleri and X. borealis. The dodecaploids are closely related to species that have DM‐W (Fig. 3) so our failure to amplify this gene in these species could be either due to sequence divergence or to loss of DM‐W. We put extensive effort into attempts to amplify DM‐W in X. borealis using a total of 20 primer combinations that independently targeted each coding exon of DM‐W over a gradient of different annealing temperatures. As a positive control, amplification of both paralogs of DMRT1 from X. borealis was successful. We therefore suspect that DM‐W may not be present in this species or potentially in other closely related species in Clade 2 in Figure 3 (i.e., X. muelleri and an undescribed tetraploid species).

Because DM‐W is a chimeric protein formed only partially from DMRT1 (Yoshimoto et al. 2008), evolutionary relationships among homologous portions of these genes (Fig. 1A) can provide insights into the origin of DM‐W. To this end, we sequenced or downloaded homologous portions of DMRT1 from paralogs of the tetraploid species X. laevis, X. borealis, and X. muelleri (accession numbers NM_001096500, NM_001085483, HQ225638–41), and compared these sequences to the homologous region of X. laevis DM‐W (accession number AB259777). For phylogenetic analysis, the AIC selected the general time reversible model with a gamma distributed among site rate heterogeneity parameter, and base frequencies fixed at equal frequencies. Phylogenetic analysis points to a much more recent origin than was previously suggested (Fig. 4), and indicates with strong statistical support that the 5′ portion of DM‐W is closely related to DMRT1β of Xenopus tetraploids. (The designation of paralogs as α or β is arbitrary and we followed Genbank annotations here). This phylogeny establishes a recent origin of DM‐W after divergence of Silurana and Xenopus. Because the relationships within the clade containing DM‐W and DMRT1β are unresolved, we cannot determine from this analysis whether DM‐W appeared before or after divergence of the ancestor of X. laevis from the ancestor of (X. borealis+X. muelleri). In an attempt to address this, we sequenced portions of the intron between exon 2 and 3 in DMRT1 paralogs and DM‐W. However, extensive divergence and insertion/deletion polymorphisms prevented us from confidently assessing homology (data not shown), so we restricted our analysis to the coding region as presented above.

4

DM‐W originated after divergence of Silurana and Xenopus. Evolutionary relationships of X. laevis DM‐W with respect to homologous of DMRT1 paralogs (α and β) of the tetraploid species X. laevis, X. borealis, and X. muelleri, and the diploid species S. tropicalis show that this gene is most closely related to Xenopus DMRT1β. DMRT1 from Bufo marinus is used as an outgroup. Nodes with posterior probabilities ≥0.99 are indicated with black dots. This analysis fails to resolve whether DM‐W originated before or after the divergence of ancestors of X. laevis and (X. borealis+X. muelleri). To evaluate alternative scenarios for the origin of DMW, we filtered the post‐burn‐in posterior distribution of tree topologies with constraint trees A–C. The posterior probabilities of constraints A, B, and C were 0.38, 0.245, and 0.375, respectively.

We explored alternative scenarios for the origin of DM‐W by filtering the post‐burn‐in posterior distribution of tree topologies recovered from the Bayesian analysis with constraint trees depicted in Figure 4. Constraint A retains topologies in which DMRT1β of X. laevis, X. borealis, and X. muelleri are monophyletic; these topologies are consistent with an origin of DM‐W before the divergence of the most recent common ancestor (MRCA) of X. laevis and X. borealis. Constraint B retains topologies in which X. laevis DMRT1β and all of the DM‐W sequences are monophyletic; these topologies are consistent with an origin of DM‐W after the divergence of the MRCA of X. borealis and X. laevis. Constraint C retains topologies in which DMRT1β, of X. borealis and X. muelleri and all of the DM‐W sequences are monophyletic. Topologies consistent with Constraint C require an explanation involving ancestral polymorphism or gene transfer via hybridization because they would suggest an origin of DM‐W after the divergence of the MRCA of X. laevis and X. borealis, but in the X. borealis lineage, even though the X. laevis lineage carries DM‐W. The posterior probabilities of constraints A, B, and C were 0.38, 0.245, and 0.375, respectively. This indicates that there is insufficient statistical power to distinguish between these possibilities with the available data. We repeated this exercise using additional data from partial sequences from DMRT1 paralogs from other species (discussed below) and recovered essentially identical results (Fig. S1).

Overall, these analyses pin down the origin of DM‐W to a window of time after divergence of Silurana+Xenopus but before divergence of X. clivii+ Clade 1 (Fig. 3). Rough divergence time estimates based on mitochondrial DNA (Evans et al. 2004) suggest that DM‐W therefore arose about 32–64 million years ago. Another divergence time estimate based on nuclear DNA suggests DM‐W could be even younger, about 13–38 million years ago (Chain and Evans 2006).

DM‐W EVOLVED UNDER NATURAL SELECTION

How did DM‐W seize control of the sex determination system? One possibility is that this was catalyzed by natural selection on DM‐W or a linked region (e.g., Rice 1986; van Doorn and Kirkpatrick 2007; Vuilleumier et al. 2007). To explore this possibility, we used two approaches to test for evidence of natural selection on DM‐W and/or the NRW.

We first analyzed the rate ratio of nonsynonymous to synonymous substitutions per site (ω) and recovered significant support for positive selection during the ascendancy and evolution of DM‐W. Using the codeml module of PAML (Yang 1998), we performed a branch‐site test for positive selection (Yang et al. 2005) on a subset of the data analyzed in Figure 4 including homologous regions of both DMRT1 paralogs from X. laevis and X. muelleri, DMRT1 paralog α of X. borealis, DMRT1 orthologs of B. marinus and S. tropicalis, and the DM‐W gene of X. laevis. DM‐W sequences from the other species (X. gilli, X. vestitus, X. largeni, X. pygmaeus, X. itombwensis, X. andrei, and X. clivii) were excluded because of missing data surrounding the junction of exons 2 and 3. Xenopus borealis DMRT1 paralog β was excluded because it had premature stop codons at amino acid positions 36 and 59.

This test compares models in which four classes of sites are considered. Sites in one class evolve under purifying selection, sites in another class evolve neutrally, and sites in two other classes evolve under positive selection. The proportion of sites in each class is estimated separately for different portions of the tree (in this case for the DM‐W branch and for the rest of the topology) and the likelihood of this model is compared to a null model in which no positive selection occurs and sites either evolve neutrally or under purifying selection (Yang et al. 2005). Comparison of the alternative model to the null model recovers significant evidence of positive selection (P < 0.01, Table 1). In the alternative model, the maximum likelihood estimate of the proportion of sites under positive selection is 0.134 with ω of these sites equal to 10.058.

1

Codon analysis of DMW and DMRT1 supports positive selection on sites in the 5′ region of DMW. Results are shown for a site‐branch test, a branch test, and a branch test for neutrality, including the likelihood of the null model (−lnL Ho), likelihood of the alternative model (‐lnL Ha), degrees of freedom (df), and P‐value (P). The P‐value of the site‐branch test is based on a 50:50 mixture of point mass 0 and χ₁². P‐values for the other tests are based on a standard) χ₁² approximation.

Model	−lnL Ho	−lnL Ha	df	P
Site‐branch	877.363	873.68	1	<0.01
Branch	889.151	879.181	1	<0.00001
Test for neutrality	879.256	879.181	1	0.698

Analyzing the same data using a model with a different ω across all sites on the DM‐W branch is significantly preferred over a model in which ω is the same in DM‐W and DMRT1 (P < 0.00001, Table 1). The maximum likelihood estimate of the ω over all sites on the branch leading to X. laevis DM‐W is 0.8152 as compared to a ratio of 0.0715 across all sites over the rest of the phylogeny. This is also more likely than the neutral expectation (where ω= 1 for the DM‐W branch), but not significantly so (P= 0.698, Table 1). In this analysis, a total of 98 codons were analyzed, 70 of which occur in the DM domain. Of 14.67 inferred nonsynonymous substitutions, slightly fewer occurred in the DM domain than expected (nine substitutions), but not significantly so (P= 0.4279, G test). By contrast, over a similar period of time only four or five amino acid substitutions occurred in the 98 codon region of X. muelleri DMRT1β and X. laevis DMRT1β, respectively, and only two occurred in X. muelleri DMRT1α and X. laevis DMRT1α.

Interestingly, there is considerable variation among species in the amino acid sequence of DM‐W. Seven of 111 codons sequenced from X. laevis and X. gilli encoded different amino acids in DM‐W and these divergent codons were evenly distributed among the exons, including in the DM domain. X. laevis and X. gilli diverged ∼17 million years ago (Evans et al. 2004). In X. clivii DM‐W from exon 2, 8 of 55 amino acids are diverged from the homologous portion of X. laevis DM‐W. These species diverged ∼32 million years ago (Evans et al. 2004).

We then used the Hudson, Kreitman, Aguade (HKA) test (Hudson et al. 1987) to test for more recent evidence of natural selection based on a multilocus polymorphism dataset, and again recovered a significant departure from neutrality. We collected polymorphism data from 23 genes from a natural population of X. laevis and X. gilli from Cape Province, South Africa. Although both species are tetraploid, we sequenced loci that are either single copy due to post‐polyploidization deletion, or with sufficiently diverged primer sites that made possible amplification of only one pair of alleles. We compared silent site polymorphism at these loci to silent polymorphism at DM‐W and linked female‐specific flanking regions of the NRW (amplified using primers designed from accession AB365520) in both species (Table 2).

2

Polymorphism data from a natural population of X. laevis from South Africa including names of unlinked loci (locus), the number of base pairs sequenced (bp), number of silent sites (silent), number of segregating sites for X. laevis and X. gilli (SXL and SXG, respectively), number of chromosomes sequenced for X. laevis and X. gilli (nXL and nXG respectively), divergence from the outgroup (diverge), and Genbank accession numbers (Accession). For DMW, different amounts of polymorphism data were obtained for X. laevis and X. gilli and separated by a forward slash, respectively. Acronyms for locus names refer to names in Supporting information.

locus	bp	silent	SXL	nXL	SXG	nXG	diverge	Accession
AR	402	99	0	30	0	32	5	HQ220654‐84
BTB domain protein 6	525	121	0	26	2	24	6	HQ220685‐709
LOC100487412	481	458	–	–	3	30	34	HQ225642‐57
efcab5	483	473	30	16	4	14	17	HQ220931‐42
femlc	497	112	1	20	4	26	7	HQ221309‐31
LOC100488701	581	581	19	30	0	18	14	HQ220994‐1016
gja3	717	717	19	14	–	–	21	HQ220943‐50
kiaa1919	660	144	0	10	2	26	7	HQ220951‐64
LOC100145442	267	267	9	18	–	–	15	HQ220965‐93
LOC100158283	690	140	15	14	0	26	15	HQ220974‐93
mastl	539	119	2	24	5	26	8	HQ221046‐70
mogs	705	170	4	20	0	24	13	HQ221071‐92
nfil3	337	77	1	24	0	24	3	HQ221093‐120
pcdh1	507	122	1	20	2	26	5	HQ221121‐43
pigo	494	126	3	22	1	26	4	HQ221144‐67
prmt6	666	157	5	22	6	22	13	HQ221168‐90
RAG‐2	1062	237	11	22	5	18	18	HQ221191‐210
rassf10	489	98	3	24	3	26	7	HQ221332‐56
zbed4	472	110	2	24	3	22	4	HQ221262‐84
LOC100496658	435	435	7	32	0	26	20	HQ221017‐45
C7orf25	532	119	2	24	0	24	4	HQ 220710‐32
sugp2	463	113	3	24	0	26	7	HQ221211‐35
znf238.2	526	118	0	24	2	24	11	HQ221285‐308
DMW	2781/1835	2549/1677	0*	12	5	14	58/49	HQ220822‐47;
								HQ220855‐76;
								HQ220879‐930

*

Two segregating sites were observed in an insertion that occurred in an ancestor of the South African population of X. laevis after divergence from X. gilli.

In the South Africa population of X. laevis, very little polymorphism was detected at the NRW; over 2500 silent sites from the NRW contained only two polymorphisms (Table 2). Both polymorphisms were in an insertion/deletion that was present in the South Africa population of X. laevis but absent in X. gilli. Because gaps were deleted to estimate divergence with respect to an outgroup, this meant that there were no polymorphisms in the region of the X. laevis NRW that was homologous to the outgroup (X. gilli). These data departed significantly from the neutral expectation (P= 0.0034, degrees of freedom (df) = 1, mlHKA). The maximum likelihood estimate of the DM‐W inheritance scalar was 0.003 as compared to the neutral expectation of 0.250. When the two polymorphisms with no outgroup homology are included in the analysis, departure from the neutral expectation is not significant (P= 0.319, df = 1) although the maximum likelihood estimate of the inheritance scalar (0.079) is still below the neutral expectation. Both X. laevis polymorphisms were completely linked and had intermediate frequency (five of 12 individuals).

To further understand the evolution of polymorphism on the X. laevis NRW, we sequenced portions of the NRW from a divergent population of X. laevis from the Democratic Republic of the Congo (DRC). By sequencing the NRW from this more closely related population, we were able to determine that the insertion containing the polymorphisms in the South Africa X. laevis population occurred after divergence from the DRC population, as opposed to the alternative that a deletion occurred in the X. gilli lineage. This is because, like X. gilli, the X. laevis population from the DRC lacks the region carrying the two NRW polymorphisms in the X. laevis population from South Africa. Twenty‐nine of 1549 sites were diverged between the two X. laevis populations in the 5′ upstream region of DM‐W. No polymorphic sites were detected in the NRW of 11 individuals from the DRC population (data not shown). This could be due to small population size, natural selection, or both, but we did not investigate these possibilities here.

In contrast to the South Africa population of X. laevis, the NRW of the X. gilli population was highly polymorphic. The NRW of the X. gilli contained five polymorphic sites out of about 1700 silent sites, even though polymorphism at autosomal loci was similar in magnitude to the South Africa population of X. laevis (Table 2). In all of the X. gilli NRW polymorphisms, the rare allele was derived relative to X. laevis, and four of these rare derived substitutions were present in only one individual. These data from X. gilli do not depart significantly from the neutral expectation (P= 0.260, df = 1) and the maximum likelihood estimate of the inheritance scalar (0.691) was greater than the neutral expectation of 0.250.

Out of eight substitutions in the DM domain of DM‐W, four were at positions that were invariant (conserved) in the DM domains of DMRT1 in the other species we examined, and four were in positions that were variable (non‐conserved) in the DM domains of DMRT1 in these species. By comparison, the DM domain of the DMY gene of the fish O. latipes, a male‐specific sex determining gene, has 13 amino acid substitutions in this region relative to its O. latipes DMRT1 paralog. Two are in positions that are conserved in other DMRT1 loci, and the others are at non‐conserved positions in DMRT1. Thus although DMY is more diverged from its closest DMRT1 paralog than is DM‐W in terms of the number of amino acid substitutions, DM‐W has more divergence at highly conserved sites in DMRT1 than does DMY.

Taken together, our tests for natural selection uncover a dynamic evolutionary history of DM‐W. It appears that the 5′ portion of this gene was subject to positive selection during early evolution that increased ω of some sites beyond the neutral expectation. More recently, DM‐W and linked portions of the NRW continued to evolve non‐neutrally in X. laevis leading to a dearth of polymorphism in this region of the genome. In X. gilli, however, recent evolution of the NRW did not depart from the neutral expectation.

BIASED PSEUDOGENIZATION OF DMRT1 PARALOGS AFTER ALLOPOLYPLOIDIZATION

How did the male‐related gene DMRT1 respond to the ascendance of DM‐W and genome duplication? One step to further understanding the evolution of this system is quantify how many paralogs of DMRT1 are functional in polyploid species of Xenopus. To accomplish this, we amplified, cloned, and sequenced portions of DMRT1 paralogs in 18 polyploid species (Table 3). We focused our efforts on exons 2 and 3 because these 5′ regions contain the DM domain, which is key to the DNA‐binding function of the protein (Yoshimoto et al. 2010).

3

DMRT1 paralogs cloned from representative tetraploid, octoploid, and dodecaploid species oKenopus. Cloned exons that appear functional are indicated with a “Y,” missing data indicated with a dash, and degenerate exons (with premature stop codons, frameshift mutations, or large deletions) are indicated with a “D” with a forward slash dividing results from exon 2 and 3 because each exon was independently amplified for cloning. Tw paralogs (α, β) are expected in tetraploids, four (α1, α2, β1, β2) in octoploids, and six (α1, α2, α3, β1, β2, β3) in dodecaploids. In some instances paralogous divergence could not be positively distinguished from allelic variation and is indicated with a “?.” Also indicated are the proportion of paralogs with data that appear functional at the coding level (proportion) and Genbank accession numbers (Accession).

Tetraploid species	α	β	Proportion	Accession
				NM_001096500,
X. laevis	Y/Y	Y/Y	2 out of 2	NM 001085483
X. largeni	Y/−	Y/− and D*/−	2 out of 2	HQ220775‐7
X. muelleri	Y/Y	Y/Y	2 out of 2	HQ225640‐1
X. borealis	Y/Y	D/Y	1 out of 2	HQ225638‐9
				HQ220733‐4,
X. fraseri	Y/−	D/Y	1 out of 2	HQ220792
				HQ220737;
X. pygmaeus	−/Y	D/Y	1 out of 2	HQ220790‐1
X. new tetraploid 1	Y**/−	Y/−	2 out of 2	HQ220778‐83
X. clivii	Y**/Y	−/Y	2 out of 2	HQ220735‐6;
				HQ220788‐9
Octoploid species	α1	α2	β1	β2
						HQ220738‐42;
X. vestitus	Y/Y	Y/Y	Y/Y	D***/−	3 out of 4	HQ220793‐5
						HQ220745‐9;
X. itombwensis	Y/Y	Y/Y	Y/Y	Y/Y	4 out of 4	HQ220796‐800
						HQ220750‐2;
X. boumbaensis	Y/Y	?/Y	Y/Y	−/−	3 out of 3	HQ220801‐4
X. andrei	Y/−		Y/−	Y/−	3 out of 4	HQ220772‐4;
						HQ220753‐4;
X. amieti	−/Y	Y/Y	Y/Y		3 out of 3	HQ220805‐7
X. wittei	Y/−	Y/−	D/−	Y/−	4 out of 4	HQ220784‐7
						HQ220743‐4;
X. new octoploid 2	−/Y	−/Y	Y/Y	Y/?	3 or 4 out of 3 or 4	HQ220818‐21
Dodecaploid species	α1	α2	α3	β1	β2	β3
								HQ220755‐61;
X. longipes	Y/Y	Y/Y	Y/Y	Y/Y	Y/Y	D/Y	5 out of 6	HQ220808‐13
								HQ220762‐6;
X. ruwenzoriensis	Y/−	Y/−	Y/−	Y/Y	−/−	D/−	4 out of 5	HQ220814
								HQ220767‐71;
X. cf. boumbaensis	Y/Y	Y/−	−/−	D/Y	−/?	Y/−	3 out of 4 or 5	HQ220815‐7

*

Segmental duplicate of exon 2 paralog β is degenerate.

**

Segmental duplicate of exon 2 paralog α not degenerate.

***

X. vestitusβ2 is heterozygous for stop codon in exon 2.

Copy number and evolutionary relationships of DMRT1 paralogs are generally consistent with duplication by whole‐genome duplication (Table 3). As expected, two DMRT1 paralogs were usually detected in tetraploid species, up to four were detected in octoploid species, and up to six were detected in dodecaploid species. In some species, fewer paralogs were detected than expected; these could have been lost by deletion or not cloned in our assay. In three tetraploid species, an extra divergent sequence was identified (Table 3); these could stem from segmental duplication of one paralog or possibly be examples of very high allelic variation. Two of the extra divergent sequences were DMRT1α paralogs and one was a DMRT1β paralog.

Translation of the cloned DMRT1 sequences indicates multiple paralogs of multiple species are degenerate due to stop codons, frameshift mutations, or large deletions (Table 3); pseudogenization therefore has a marked impact on the stoichiometry of DMRT1 and DM‐W in Xenopus females. Eight instances of degeneration were observed in DMRT1β and a ninth in a putative singleton duplicate copy of DMRT1β in X. largeni. To our surprise, however, degeneration was never observed in DMRT1α or in either of the putative singleton duplicates of DMRT1α (Table 3). Under the assumption of allopolyploidization, the “β” paralogs are inherited from one of the two diploid ancestors whose genomes were fused to form an ancestral allotetraploid Xenopus. The eight to nine instances of degeneration were observed in exon 2 of DMRT1β (and never in exon 3), and each instance was species specific (e.g., a premature stop codon was never observed in the same position in multiple species).

Using BayesTraits (Pagel and Meade 2006), we recovered support for significant phylogenetic bias in the pattern of pseudogenization in DMRT1. When missing data were coded as missing, a model with a different rate of pseudogenization in DMRT1 paralog α and β (−lnL =−30.8618) was significantly preferred over a model with one rate of pseudogenization in both paralogs (−lnL =−34.3794, P= 0.00776, df = 1). In the more parameterized model the rate of pseudogenization was 10 times faster in DMRT1 paralog β than in DMRT1 paralog α. As a conservative measure, we also performed this test with missing data coded as functional (there were more missing data from DMRT1β than DMRT1α), and the result was the same (P= 0.00781). This result is consistent with the notion that natural selection preferentially targeted DMRT1β for pseudogenization in multiple polyploid species or that natural selection preferentially favored DMRT1α to remain functional. A caveat is that we may have underestimated the number of DMRT1 pseudogenes if degenerative mutations are present regions that we did not sequence, if seemingly functional paralogs are not expressed due to mutations in the regulatory region, or if we failed to clone some pseudogenized or deleted paralogs.

DMRT1α EXPRESSION IS REGULATED BY DISTINCT MECHANISMS OVER DEVELOPMENT

Is it possible that the pattern of biased pseudogenization in DMRT1β is somehow linked to the origin of DMW from partial duplication of DMRT1β? To explore this question, we examined mechanisms of expression divergence (cis‐ and/or trans‐) of a DMRT1α, in X. laevis and X. borealis during two developmental stages (tadpole and adult). As discussed earlier and in Table 3, DMRT1α paralogs appear to be fully functional at the sequence level in both of these species.

A first concern in this analysis is whether expression levels in hybrid individuals differ depending on the direction of cross (i.e., whether the mother of the hybrid individuals is X. laevis or X. borealis). However, tests for parent of origin effects were not significant (Supporting information) so data from both directions of hybrid cross were considered jointly. Another concern is whether expression levels differ in males and females. In hybrid tadpoles with a X. borealis mother, we were not able to test for sex effects because we do not know the sex‐determining gene in this species. However, we did test for sex effects in hybrid tadpoles when the mother was X. laevis (Supporting information); these tests also were not significant after Bonferroni correction, which is consistent with Yoshimoto et al. (2008), allowing us to combine expression data from male and female tadpoles. Because the sex effect was individually significant before Bonferroni correction, we repeated analysis presented below using only males and conclusions were essentially identical (Supporting information).

Comparison of expression levels of alleles from each species indicates that DMRT1α expression in X. laevis and X. borealis is significantly different during primary gonadal differentiation in tadpoles, but not in adult testis (indicated by significant “P= 0” tests in tadpoles but not adults; Table 4). Significant expression divergence of DMRT1α in tadpoles during primary gonadal differentiation is interesting from the standpoint that DM‐W is expressed in X. laevis females but possibly not in X. borealis if this species lacks this gene. Expression of DMRT1〈 is 1.74 fold higher in X. laevis than in X. borealis; this is indicated by the log₂‐transformed X. laevis/X. borealis parental ratio of 0.8027 in Table 4. Because DMRT1β is a pseudogene in X. borealis (Table 3), the disparity in DMRT1 expression between X. laevis and X. borealis during primary gonadal differentiation could be even higher if in X. laevis the DMRT1α and DMRT1β loci are co‐expressed, translated, and functionally interchangeable.

4

Log₂‐transformed parental and hybrid expression ratios (X. laevis/X. borealis) of DMRT1α and probabilities of the null hypothesis that this ratio is equal to zero (P= 0), that ratios in parental species and hybrids are equal (P= H), and that the hybrid ratio is equal to zero (H = 0). Interpretations include cis‐ divergence (C) and trans‐ divergence (T) with each mechanism causing upregulation of the X. laevis (XL) or X. borealis (XB) allele. Inferences are based on two tests, (P= H) and (H = 0), and significant departures of these test after Bonferroni correction for two tests are indicated with asterisks. Significant departure of (P= 0), is also indicated with asterisks.

	Parental ratio	P value (P=0)	Hybrid ratio	P value (P=H)	P value (H=0)	Interpretation
Tadpole	0.8027	0.0068*	−0.3982	0.0070*	0.0053*	C: XB up, T: XL up
Adult	1.1491	0.0846	3.0141	0.0250*	0.0001*	C: XL up, T: XB up

How can we account for these expression differences? In the simplest scenario, a single mutation could cause different levels of expression in two species; such a mutation could occur in cis‐ (e.g., affecting an upstream promoter) or in trans‐ (e.g., affecting expression level or binding capacity of a transcription factor). If DMRT1α expression were governed by the same regulatory elements throughout development and across multiple tissue types, one would expect a similar direction of divergence (e.g., upregulation of DMRT1α in X. laevis compared to X. borealis) and a similar relative magnitude of divergence (e.g., twofold) in all developmental stages and all tissue types where the gene is expressed. Contrary to these expectations, we found that DMRT1α regulation occurs in a developmental‐stage specific manner. In tadpoles, DMRT1α upregulation of X. laevis compared to X. borealis alleles is attributable to trans‐ acting factors (this is indicated by a significantly higher parental ratio than the hybrid ratio in tadpoles according to the “P= H” test; Table 4). But this is counteracted to some degree by antagonistic cis‐ upregulation of the X. borealis DMRT1α allele (because the log₂‐transformed DMRT1〈 hybrid tadpole ratio is significantly higher than 0 according to the “H = 0” test; Table 4). In adult testis, however, the opposite mechanisms operate: trans‐ divergence upregulates the X. borealis allele (the adult hybrid ratio is significantly higher than the adult parental ratio according to the “P= H” test; Table 4) and cis‐ divergence upregulates the X. laevis allele (the adult hybrid ratios are both significantly higher than zero according to the “H = 0” test; Table 4).

In X. laevis, DMRT1 is autosomal (Uno et al. 2008; Yoshimoto et al. 2008). But if DMRT1 inheritance is not biparental in X. borealis (e.g., if a sex‐linked allele of one paralog was deleted), this could affect conclusions relating to mechanisms of expression divergence in tadpoles because some X. borealis tadpoles would be hemizygous for one DMRT1 paralog, and because some hybrids would lack one DMRT1 paralog from X. borealis. To investigate this possibility, we co‐amplified and sequenced exon 2 of DMRT1α and DMRT1β from genomic DNA of F1 hybrids generated by crossing an X. borealis female to an X. laevis male, and also from the reciprocal cross with an X. laevis mother. If both paralogs of DMRT1 are biparentally inherited in both species, then hybrids from both crosses should both inherit an allele of both paralogs from both species. Alternatively, if both alleles of one DMRT1 paralog were not biparentally inherited in X. borealis, 50% of the F1 hybrids from one of these crosses would inherit an allele from only one paralog from X. borealis. This is because they would either receive different sex chromosomes from their X. borealis mother (if females are heterogametic) or from their X. borealis father (if males are heterogametic). We used adult female and male hybrids (12 females, six males from the cross with the X. borealis mother and seven of each sex from the cross with the X. borealis father) to test this. Species‐ and paralog‐specific single nucleotide polymorphisms allowed us to confirm that F1 hybrids inherit both alleles from DMRT1 paralogs from both species. We note that biparental inheritance of both alleles of both DMRT1 paralogs in X. borealis does not necessarily mean that DMRT1 is not sex‐linked in this species.

Overall then, analysis of DMRT1α expression divergence illustrates two main points with respect to the evolution of variation in mechanisms of sex determination in clawed frogs. First, DMRT1α is expressed significantly higher during gonadal differentiation in a DM‐W‐containing species (X. laevis) compared to one that may lack DM‐W (X. borealis) but not significantly different in adult testis. Second, genetic mechanisms that govern expression of DMRT1α during primary gonadal differentiation in tadpoles are at least partially independent from those that govern expression in adult testis. We relate these expression findings to our results concerning biased pseudogenization of DMRT1 below.

Discussion and Conclusions

Almost every vertebrate species has separate sexes, but mechanisms by which sexual differentiation is orchestrated vary considerably. In frogs, variation in male versus female heterogamy is observed among frog families and genera (Hillis and Green 1990) and within the species Rana rugosa (Miura 2008). Presence/absence polymorphism of a W sex chromosome occurs in Leiopelma hochstetteri (Green et al. 1993), and mapping of sex‐linked genes demonstrates that ranids and pipids have different sex chromosomes (Uno et al. 2008). Here we illustrate that mechanisms of sex determination also vary within the African clawed frogs (Xenopus+Silurana)—and possibly within the genus Xenopus—as a consequence of the appearance of a novel female‐specific sex‐determining gene called DM‐W. In X. laevis, DM‐W triggers primary (gonadal) female development, and this may occur via competitive inhibition in females of genes that are activated by DMRT1 in males (Yoshimoto et al. 2008, 2010). In contrast to analysis of the full DM‐W and DMRT1 genes (Yoshimoto et al. 2008), analyses presented here demonstrate that this gene evolved in Xenopus after divergence from Silurana but before divergence of X. laevis and X. clivii. That we did not succeed in amplifying DM‐W in X. borealis despite considerable effort opens the possibility that this and other closely related species in Clade 2 (Fig. 3) do not have DM‐W. Similar lines of evidence suggested that the mammalian sex‐determination locus SRY was not present in platypuses (Grützner et al. 2004) before confirmation with the complete genome sequence (Warren and others 2008), or in various vole species (Just et al. 1995), and that DMY (a paralog of DMRT1) is not present in species closely related to the medaka (Kondo et al. 2003, 2004).

At this point, it is not clear whether DM‐W evolved in a diploid genome with 18 chromosomes or a tetraploid genome with 36 chromosomes. The probability of fixation of DM‐W in a diploid versus a tetraploid is presumably influenced by multiple factors such as dominance interactions between DM‐W and the ancestral trigger for sex determination, whether a tetraploid genome is polysomic or disomic (Otto and Whitton 2000), whether DM‐W arose on the same chromosome or a different chromosome as the ancestral sex‐determining locus, whether DM‐W conferred advantages to one or both sexes, and whether the ancestral heterogametic sex was female or male. If DM‐W did evolve after divergence of the most recent common ancestor of X. laevis and X. borealis, various scenarios depicted in Figure S2 involving multiple episodes of allo‐tetraploidization in Xenopus are possible. Identification of the sex‐determining locus in S. tropicalis, and a complete genome sequence for X. laevis will assist in exploring these possibilities. Also relevant to the question of exactly when DM‐W arose is the relationship among X. clivii and Clades 1 and 2 in Figure 3. Most relevant is the relationship among genomic regions linked to DM‐W and DMRT1β in these lineages. We attempted to address this question by sequencing portions of the intron between exons 2 and 3 of DM‐W and DMRT1β in X. laevis, X. borealis, X. muelleri, and X. clivii but divergence prevented us from making reliable homology statements for these sequences. Analysis of the available data with unambiguous alignment provided insufficient statistical support to distinguish between an origin before or after divergence of the MRCA of X. laevis and X. borealis.

Multiple mechanisms have been proposed to account for variation among species and within populations in the mechanisms of sex determination. For example, sexual antagonism and other types of genetic conflict could potentially drive sex chromosome turnover (Rice 1986; Werren and Beukeboom 1998; van Doorn and Kirkpatrick 2007). Sex chromosomes have a smaller effective population size than the autosomes (barring extreme sex‐specific demography) and genetic drift thus has a more profound impact on their evolution. Genetic drift, possibly combined with sex‐ratio selection, could also contribute to sex chromosome turnover (Bull 1983; Kozielska et al. 2006; Vuilleumier et al. 2007). Another possibility is that sexual selection and addition of upstream elements drives evolution of sex determination systems (Pomiankowski et al. 2004; Wilkins 1995). Consistent with a role for natural selection in sex chromosome turnover and evolution, we found that some sites in the 5′ portion of DM‐W evolved under positive selection since divergence from DMRT1β and that DM‐W and linked regions have a significantly lower level of polymorphism in the South Africa population of X. laevis. Other examples of rapidly evolving sex determining loci have been reported (Tucker and Lundrigan 1993; Whitfield et al. 1993; Zhang 2004; Graves 2008), but it is still unclear whether there is a general role for natural selection during sex chromosome turnover and evolution, or whether purifying selection on these loci tends to be relaxed.

Arguing against relaxed purifying selection on DM‐W, we found a dearth of molecular polymorphism in the NRW of a natural population of X. laevis from South Africa. The NRW of a DRC population of X. laevis is diverged from the South African population, so it appears that non‐neutral evolution of DM‐W or linked portions of the NRW occurred in the South African population after divergence from the DRC population, as opposed to being a species‐wide phenomenon. It is also worth pointing out that the DRC population may warrant separate species status (Measey and Channing 2003), in which case one would not expect a selective sweep in the X. laevis from South Africa to involve the DRC population. An alternative demographic explanation for these data is a higher variance in reproductive success in females than males (Lande and Barrowclough 1987; Engen et al. 2007). This could happen, for instance, if females skipped some mating cycles to develop a full clutch of eggs. However, this was not the case in a population of X. gilli, which had a higher level of polymorphism than the neutral expectation (although not significantly so).

MARGINALIZED EXPRESSION OF DMRT1β?

In general, the probability that one paralog of a duplicate pair is a pseudogene increases with the age of the duplicates (Lynch and Conery 2000). Given that ∼60% of duplicate genes in the pseudotetraploid species X. laevis have degenerated to singletons (Sémon and Wolfe 2008), it is no surprise that some DMRT1 paralogs of the Xenopus species we surveyed are pseudogenes. What is surprising, however, is that pseudogenization of DMRT1 paralogs occurred in a phylogenetically biased pattern in multiple polyploid species, affecting only paralogs closely related to DM‐W (that is, the DMRT1β paralogs). A similar pattern of biased pseudogenization was observed at the RAG1 locus in African clawed frogs (Evans 2007). How can these findings be explained?

Analysis of DMRT1 in S. tropicalis indicates that this gene is expressed during much of tadpole development, including well before gonadal differentiation and adult testis (Yoshimoto et al. 2006). A broad expression pattern of DMRT1 is also observed in X. laevis based on real‐time PCR assays (Osawa et al. 2005; Yoshimoto et al. 2006, 2008), suggesting that a broad expression domain is the ancestral condition of DMRT1 in African clawed frogs. However, these X. laevis assays probably jointly quantified expression of DMRT1α and DMRT1β, and it is therefore possible that the expression domain of each of these paralogs is diverged.

One explanation for these observations is that the expression domain of DMRT1β was marginalized soon after it was formed by whole‐genome duplication, making this paralog dispensable in Xenopus polyploids. Consistent with this speculation, our results demonstrate that expression of DMRT1α is orchestrated by a nonidentical suite of regulatory mechanisms during development that are therefore distinct and independently mutable (Table 4). If this developmentally compartmentalized regulation is the ancestral condition for both DMRT1 paralogs in Xenopus, this opens the possibility that mutation could have degenerated regulatory machinery of DMRT1β that drives expression in adult testis and other tissues, while leaving intact expression during primary gonadal differentiation. Also consistent with the hypothesis that the expression domain of DMRT1β was marginalized soon after genome duplication is the observation that the expression domain of DM‐W—a gene formed from partial gene duplication of the 5′ region of DMRT1β and potentially also portions of its regulatory region—is restricted to a brief window during primary differentiation in the female gonad (Yoshimoto et al. 2008). Thus, we speculate that the narrow expression domain of DMW could offer clues into why DMRT1β paralogs appear to be dispensable in Xenopus polyploids.

Future work aimed at independently characterizing the expression domain of DMRT1α and DMRT1β in multiple Xenopus species could further test the hypothesis that the ancestral expression domain of the progenitor locus of DMRT1β and DMW was marginalized prior to the origin, functional diversification, and biased pseudogenization of DMRT1β. It would also be useful to evaluate whether the upstream portion of DMRT1β is homologous to that of DM‐W when X. laevis genome sequences are available. We were not able to detect homology between available sequences from the NRW of X. laevis and sequences upstream of S. tropicalis DMRT1 obtained from the genome sequence (Hellsten et al. 2010).

Why is DMRT1α expression in tadpoles upregulated in X. laevis compared to X. borealis? At this point, we do not know, but one possibility is that this is a consequence of sexual antagonism in X. laevis where upregulation of DMRT1 is favored in males. If this is the case, developmental compartmentalization of DMRT1α regulation could have facilitated regulatory response to the invasion of DM‐W by decoupling regulatory control in tadpoles and adults (we did not detect significant difference in expression level of DMRT1α in adult X. laevis and X. borealis). Alternatively this expression divergence could be unrelated to DM‐W and instead be due to neutral drift, or related to other physiological differences between X. laevis and X. borealis.

We did not consider mechanisms of expression divergence between X. laevis and X. borealis DMRT1β because of premature stop codons in the X. borealis paralog. However, in X. muelleri, a close relative of X. borealis, the reading frames of both DMRT1 paralogs appear intact. If Clade 2 (Fig. 3) really does lack DM‐W, an interesting direction for further work would compare expression levels and mechanisms of expression divergence of both paralogs of this species to a DM‐W‐containing species such as X. laevis during gonadal differentiation. Another question that remains unanswered is how expression stoichiometry of DMRT1 and DM‐W expression varies in species with higher ploidy levels at the RNA and protein level.


Associate Editor: C. Peichel

ACKNOWLEDGMENTS

We thank F. Chain for identifying singleton genes in X. laevis, P. Wittkopp for providing SAS scripts for analysis of expression divergence, A. Smith for assistance with pyrosequencing, and Catherine Peichel and anonymous reviewers for constructive comments.

Author notes

Current address: Center for Ecology and Evolution, 335 Pacific Hall, 5289 University of Oregon, Eugene, Oregon 97403‐5289.