Molecular evidence for sex reversal in wild populations of green frogs (Rana clamitans)

Frog sampling

We focused on adult frogs to complement prior studies on intersex (Murphy et al., 2006; Skelly, Bolden & Dion, 2010; Smits, Skelly & Bolden, 2014), though are careful to distinguish GSD + EE here from sequential hermaphroditism which occurs in some fishes which sex reverse as adults (Bachtrog et al., 2014). Sexually-discordant adult frogs would have established their phenotypic sex opposite their genotypic sex prior to metamorphosis and not as adults. From April 25 to July 21, 2016 we collected adult phenotypic female and male green frogs from 16 ponds along a forest-suburban land use gradient (Fig. S2; Table S1). Nine ponds were studied in prior work (Smits, Skelly & Bolden, 2014; Lambert et al., 2015) and we recently acquired access to seven others. We collected frogs in buckets containing source pond water chilled on ice. Upon returning to the laboratory, we euthanized frogs with an overdose of buffered MS-222 and confirmed phenotypic sex via dissection and observation of the gonads. After euthanasia, we collected two samples of skeletal muscle from each frog in 95% ethanol and fixed each specimen in 10% buffered formalin. Pond and frog data are summarized in Table S1. In total, we collected 464 adult green frogs (μ = 29, range = 6–92) from the 16 study ponds. Of this sample, we collected a total of 129 phenotypic females (μ = 8.1, range = 2–17) and 335 phenotypic males (μ = 20.9, range = 4–75). Most phenotypic females were gravid, but a small proportion were spent, newly ovulating, or immature (Table S1). This study was approved by Yale IACUC protocols 2013-10361 and 2015-10681 and sampling was approved by CT DEEP Permit 0116019b.

Land use classification

We identified ponds surrounded either entirely by undeveloped forests or by suburban neighborhoods (Fig. 1) using modeling analyses of high-resolution remotely-sensed imagery to classify the degree of suburban land use surrounding each of our study ponds (Fig. S3). To do this, we used a 2010 statewide orthophoto of Connecticut provided by the CT Department of Energy and Environmental Protection. This orthophoto was one-m resolution, was composed of four bands (infrared and natural color), and was taken during leaf-on (Fig. S3). We performed a supervised classification using the spatial analyst tools in ArcMap 10.1 (ESRI). Our targeted land cover types included trees, lawns, buildings, and roads or other paved surfaces. We supervised the classification with representative pixels (Fig. S2) surrounding each of the 16 ponds and estimated the percent cover of each land cover type. Ultimately, we merged lawn, building, and paved surface cover types into a single “Percent Suburban” predictor variable.

Testicular histology

For histological analysis, we included all males from eleven ponds and randomly selected 22 males from the five remaining ponds with higher sample sizes (≥26 males sampled). Following Skelly, Bolden & Dion (2010), we removed the left testis from each male. At the Histology Core at the Yale School of Medicine, we retained every 20th section starting 10 sections into the gonad (Skelly, Bolden & Dion, 2010; Smits, Skelly & Bolden, 2014) and stained sections with hematoxylin and eosin. Each section was assessed for intersex (the presence of oocytes) under a Leica BF200 compound microscope. We only included samples in our analyses which had at least four high quality sections. Sections were considered high quality if the staining was consistent and the gonad was predominantly intact.

Genotyping

Tissue samples were preserved in 95% ethanol and transported to Diversity Arrays Technology Pty. Ltd. (Bruce, ACT, Australia) for genetic sequencing. We genotyped each frog at each of five recently identified sex-linked single nucleotide polymorphism (SNP) loci (RaclCT001, -003, -005, -007, and -009) (Lambert, Skelly & Ezaz, 2016) using DArTMP methods. For each locus, females are homozygous for the most common allele (i.e., X-chromosome allele) and males are heterozygous (i.e., have both a X- and Y-chromosome allele). Depending on the locus, we used either a multiplex or monoplex polymerase chain reaction (PCR) with locus-specific primers (Table S4) amplifying several base pairs around the sex-linked SNP locus. In the first round of PCR, the locus specific products were amplified in 30 cycles. In the second round of PCR, the sample-specific barcoded primers were used to identify individual frogs. The products of the second round of PCR were then sequenced on an Illumina Hiseq 2500.

Genotypic sex assignments

Many anamniote taxa have homomorphic sex chromosomes which have traditionally made it challenging identify the sex chromosomes and the master SD gene for any given taxon (Alho, Matsuba & Merila, 2010; Rodrigues et al., 2014, 2018). For amphibians, an entirely non-recombining master SD gene has been discovered in only one amphibian species (Xenopus laevis) to date (Yoshimoto et al., 2008, 2010; Kloc & Kubiak, 2014). Through a variety of approaches (candidate genes, microsatellites, GBS), several studies have identified sex-linked genetic loci in a small number of amphibian taxa (Berset-Brandli et al., 2006; Stock et al., 2011; Lambert, Skelly & Ezaz, 2016; Brelsford, Dufresnes & Perrin, 2016; Brelsford et al., 2017). Individual sex-linked loci are likely to experience recombination and will therefore inaccurately identify genotypic sex at varying frequencies when analyzed individually. However, probabilistic Bayesian statistical methods are valuable techniques which allow us to statistically assign genotypic sex from multiple sex-linked genetic markers even in the absence of direct recombination rates at each locus from parent and offspring linkage analyses (Alho, Matsuba & Merila, 2010). By extension, such models provide inference into sex reversal by identifying sexual genotype-phenotype discordance, accounting for variation among loci in genotypic sexing reliability due to recombination. We take this approach in this study.

We used Bayes’ theorem (below) to estimate the probability that each frog was a given genotypic sex and use Bayesian generalized linear models to estimate the probability that sex reversal was related to land use.

We used prior combinations of genotypes across loci (Table S5) in the original set of 23 adult females and 54 adult males used to identify these sex-linked markers (Lambert, Skelly & Ezaz, 2016). We assumed the probability of being either sex was 0.5 but found that the value of P(Female) and P(Male) do not heavily influence the outcome of this analysis.

For the frogs here, we observed 31 different combinations of genotypes across loci for the sequenced markers. A total of 18 of these combinations were present in the original set of frogs used to identify the markers. These combinations were present in most sampled frogs in the study reported here. However, thirteen combinations, representing 88 individual frogs, were not previously observed. For these 13 combinations, we were not able to use prior combinations of genotypes. Instead, we used the prior frequencies of individuals with a given number of male and female genotypes for across all markers (e.g., an individual displaying a female genotype at three loci and a male genotype at two loci) regardless of the order of each locus. In most cases, the phenotypic sex and probable genotypic sex of these individuals were concordant.

Because of the possible effects of recombination influencing interpretation of sex reversal, we performed a second, but otherwise identical, analysis using just RaclCT001 and RaclCT003 (Table S6) because through our previous study we found that these two loci confidently assigned genotypic sex on their own (Lambert, Skelly & Ezaz, 2016). We constrained our analysis to individuals which had sequence data for both loci and individuals which had genotypes observed in our previous dataset. These criteria ultimately only excluded 13 phenotypic females and eight phenotypic males.

Statistical analyses of sex

We analyzed putative associations between frequencies of sex reversal and intersex with land use within a Bayesian probabilistic framework following McElreath (2015) and using the maximum a posteriori (map) function in the R (v 3.4.0) package “rethinking.” To do so, we modeled the probability of female-to-male sex reversal, male-to-female sex reversal, or intersex as binomial responses and used flat prior distributions for intercepts and land use slopes. For each the probability of female-to-male sex reversal, male-to-female sex reversal, and intersex, we developed two models. One model only included only an intercept term and the other include an intercept term and a term for suburban land use. We similarly modeled whether intersex frequencies were associated with female-to-male, male-to-female, or total sex reversal frequencies.

We evaluated performance among models using the widely applicable information criterion (WAIC) which takes the averages of log-likelihoods over the posterior distribution of a given model. WAIC values are particularly helpful inference metrics for evaluating Bayesian probabilistic generalized linear models because they do not rely on a Gaussian posterior distribution, therefore making them useful for binomial models like those used here. Furthermore, we assessed the relative predictive power of different models by comparing changes in WAIC (dWAIC in Table S3) and weighted model probabilities using Akaike weighted ranking (McElreath, 2015).

We note that, because sex reversal likely occurs during larval development and because the adult population represents multiple aggregated cohorts, it would be challenging to interpret any associations or lack of associations between sex reversal frequencies and environmental variables (e.g., temperature or pH) which may vary between years. Because of this, we limited our analyses to suburban land use which has remained a relatively constant measurement of the relative anthropogenic impact among populations over the past several years.

Water chemistry

Our prior work demonstrated that the presence and diversity of endocrine disrupting chemicals capable of impacting vertebrate sexual differentiation are positively correlated with a suburban land use gradient (Lambert et al., 2015). Similarly, basic water quality parameters are also well known to be correlated with urban land use gradients (Dow & Zampella, 2000; Conway, 2007; Zampella et al., 2007; Brans et al., 2017), including in our own work (Lambert et al., 2015; Holgerson et al., 2018). As an additional layer to this study, we measured multiple water chemistry parameters to evaluate its association with the land use gradient as a relative measure of human influence among ponds. At most visits to each pond, we measured water conductivity (specific conductance), pH, and dissolved oxygen. For the first two measurements, we used an Oakton PCSTestr 35 Multiparameter probe. For dissolved oxygen, we used an YSI ProODO Handheld Optical meter. Conductivity in particular is a useful proxy for myriad contaminant sources entering urban ponds (Dow & Zampella, 2000; Conway, 2007; Zampella et al., 2007).

To characterize how land use might be associated with different water variables (conductivity, pH, dissolved oxygen, and temperature) we used linear models which included the percent of suburban land cover surrounding each pond as well as sampling date as predictor variables. If sampling date was not significant (at p < 0.05) in the model, we instead used linear mixed effects models with the “lmer” function in the lme4 package in R, treating sampling date as a random effect to account for repeated measures, and treated suburban land use as the only fixed effect. For mixed effects models, we calculated an R² with the “r.squaredGLMM” function in the R package MuMIn.

Sex phenotype-genotype discordance

We inferred sex reversal from discordances between frog phenotypic sex and genotypic sex across five sex-linked markers, finding evidence that sex reversal occurs across most populations and in both directions. From 464 frogs (Table S1) collected across 16 populations, we observed sexual phenotype-genotype discordance in 4.5% (n = 21) of frogs and in 75% of populations (Fig. 2). Nine populations exhibited a single direction of sexual discordance while three populations exhibited bidirectional sexual discordance (Fig. 2). Across populations, 8.5% (n = 11, 89% credible interval 5.4–13.4%) of genotypic females had a male phenotype, whereas 3% (n = 10, 89% credible interval 1.8–4.9%) of genotypic males had a female phenotype.

Contrary to expectations that human land use drives sex reversal, we observed no such relationships between sex reversal frequencies, in either direction, and land use (Fig. S1). We detected sex phenotype-genotype discordant frogs in all uncontaminated forested ponds. Bayesian analyses indicated that suburban land use was not associated with frequencies of sexually-discordant genotypic females (XX♂; posterior probability = 0.500, 89% credible interval 0.495–0.510) nor sexually-discordant genotypic males (XY♀; posterior probability = 0.503, 89% credible interval 0.495–0.507) (Table S5; Fig. S1). In all cases, models only containing intercepts carried the majority of Akaike model weighting compared to models including suburban land use (Table S2). Similarly, WAIC values across models for either sex reversal type were relatively similar and standard errors of each models’ WAICs were relatively high (Table S2). Importantly, changes in WAIC values (dWAIC) between the suburban model and the better model (lowest WAIC) with only an intercept were low. Additionally, the standard errors of estimated dWAIC were relatively high compared to actual dWAIC values. All metrics suggest no influence of land use on sex reversal. While it is possible that forested ponds studied here were exposed to some unmeasured contaminants, our previous work demonstrated that suburban land cover is a reliable proxy for the presence and extent of contamination in these ponds (Lambert et al., 2015).

To address the possibility that recombination could obscure the genotype-phenotype relationship leading to false positives for sex reversal, we performed a secondary analysis using a restricted dataset including only the most tightly sex-linked markers (those with alleles showing the strongest association with sex; see Methods). This additional analysis provided qualitatively identical results (Table S3). Every individual identified as sex-reversed in this restricted analysis was previously assigned as sex-reversed in the original analysis. Of the 21 individuals excluded in this secondary analysis (see Methods for exclusion criteria), three had previously been assigned as sex-reversed and the remaining 18 were sex-concordant. No individuals present in both analyses were identified as sex-reversed in the first analysis but not in the second analysis. Across populations, 7% (n = 8, 89% credible interval 4.0–12.0%) of genotypic females displayed a male phenotype, whereas 3% (n = 10, 89% credible interval 1.8–5.0%) of genotypic males were phenotypically female. As with the prior analysis, we detected sexually-discordant frogs in all uncontaminated forested ponds. Model comparisons indicated a slight advantage of the model including suburban land cover (negative effect; posterior probability = 0.487, 89% credible interval 0.475–0.50) over an intercept-only model, suggesting a slight correlation between suburban land use and female-to-male sex reversal. However, the credible interval overlapped with 0.50, all model comparison metrics (Table S3) were only modestly different between the two models, and the proportional change in log odds was only 0.95 indicating the relationship between female-to-male sex reversal and suburban land use was minimal. Bayesian analyses for sexually-discordant genotypic males (XY♀) was not influenced by suburban land use (posterior probability = 0.503, 89% credible interval 0.495–0.507) (Table S3).

Intersex

Of 246 phenotypic males histologically examined, 23 had intersex testes (9.3%, 89% credible interval 6.8–12.9%). Of the intersex frogs, sex-linked markers indicated that 91% (n = 21) were genotypically male while 9% (n = 2) were genotypically female, although a higher proportion of sexually-discordant phenotypic males (XX♂) were intersex compared to sexually-concordant phenotypic males (XY♂) (Fig. 3). While sexually-discordant phenotypic male (XX♂) frogs were limited in number, the posterior probability for genotypic sex influencing intersex frequencies was 0.28 (89% credible interval 0.10–0.60) and the proportional change in log odds of being intersex as a genotypic female (XX♂) compared to being a genotypic male (XY♂) was 0.40. This indicates that the odds of a sexually-concordant male (XY♂) being intersex is 60% lower than the odds for a sexually-discordant genotypic female (XX♂) such that 20% of sexually-discordant genotypic females (XX♂) were intersex whereas only 9% of sexually-concordant males (XY♂) were intersex. We note, though, that although XX♂ have a higher probability of being intersex, most intersex frogs were XY♂. Because most sexually-concordant males (XY♂) develop typical testicular morphology whereas many sex-reversed genotypic females (XX♂) likely retain oocytes in their testes, this result reflects differences in the frequency of XX♂ and XY♂ in wild populations.

We identified intersex phenotypic males in 10 of 16 populations (Fig. 2; Table S1). Half of ponds surrounded by suburban land cover had intersex males; for suburban ponds with intersex males, intersex frequencies ranged from 5.3% to 33.3% (Table S1). All four forested ponds harbored intersex males with intersex frequencies ranging from 4.3% to 19% of males in a given pond (Table S1). Bayesian models suggested intersex frequencies were not associated with human land use (posterior probability = 0.498, 89% credible interval 0.493–0.503; Table S2; Fig. S1). WAIC and dWAIC standards errors were again high and Akaike weighting again did not favor the model including suburban land use (Table S2; Fig. S1). This result is mirrored by data (H.E. Bement & D.K. Skelly, 2011, unpublished data) showing equal intersex frequencies in green frog populations from each of three forested (range 6.7–25% intersex) and three suburban ponds (range 4.3–17% intersex). Models also indicated that intersex frequencies were not associated with rates of female-to-male (posterior probability 0.49, 89% credible interval 0.46–0.52), male-to-female (posterior probability 0.53, 89% credible interval 0.49–0.56), or total sex reversal frequencies (posterior probability 0.51, 89% credible interval 0.47–0.54).

Water chemistry

Dissolved oxygen, pH, conductivity, and temperature all increased with increasing suburban land cover surrounding ponds. The model (full model R² = 0.29, p < 0.001) for dissolved oxygen included both Percent Suburban (p = 0.002, estimate = 0.06 ± 0.02) and Sampling Date (p < 0.001, estimate = −0.07 ± 0.02). Dissolved oxygen generally declined throughout the season by 0.7 mg/L every 10 days and, controlling for Sampling Date, dissolved oxygen increased by 0.6 mg/L for every 10% increase in landscape composition that was suburban land cover.

For pH, the model (full model R² = 0.22, p < 0.001) included both Percent Suburban (p < 0.001, estimate = 0.02 ± 0.003) and Sampling Date (p = 0.04, estimate = 0.006 ± 0.003). Generally, pH increased throughout the season by 0.06 every 10 days and, controlling for Sampling Date, pH increased by 0.2 for every 10% increase in suburban land cover composition.

In the linear regression (full model R² = 0.44, p < 0.001) for conductivity, Percent Suburban was significant (p < 0.001) but Sampling Date was not (p = 0.10). A linear mixed effects model treating Sampling Date as a random effect found a significant positive correlation between conductivity and the fixed effect Percent Suburban (R² = 0.43, p < 0.001, estimate = 6.5 ± 0.9). Conductivity increases on average by 65 microsiemens for every 10% increase in landscape comprised of suburban land cover.

For temperature, the model (full model R² = 0.52, p < 0.001) included both Percent Suburban (p < 0.001, estimate = 0.06342 ± 0.01528) and Sampling Date (p < 0.001, estimate = 0.12419 ± 0.0141). Generally, temperature increased throughout the summer by 1.2 °C every 10 days throughout the summer and, controlling for Sampling Date, temperature increased by 0.6 °C for every 10% increase in landscape composition that was suburban land cover.

These data help reconfirm the use of our suburban land use gradient as a proxy for contamination. Chemical data from forested and suburban pond studied here and in prior work can be found in Lambert et al. (2015).

Potential causes of sex reversal

The bidirectional sexual discordance observed here perhaps indicates that green frog sexual development may be influenced by multiple factors. Laboratory experiments show the direction of amphibian sexual differentiation away from the genotypically predisposed sex (e.g., genotypic female to phenotypic male or genotypic male to phenotypic female) can vary by species, demes, and environmental conditions (Miura et al., 2016; Tamschick et al., 2016a). The bidirectional sex reversal observed here in adult green frogs could possibly indicate that amphibian sexual differentiation is guided by multiple environmental (temperature and chemical) factors or yearly variation in the environment. If this were true, it may be challenging to identify specific mechanisms of sex reversal. Unfortunately, green frogs require 2–3 years post-metamorphosis to mature and survive for up to 7 years (Shirose & Brooks, 1995) and so we are unable to assess sex reversal frequencies as a function of a single year’s water conditions in the adults sampled here due to multiple overlapping cohorts. Additionally, we cannot discount that frogs may have moved from the environment which led to their intersex or sex-reversed condition before they were collected. Additional studies assessing patterns in environmental conditions with sex reversal frequencies in wild cohorts of larval or metamorphosing amphibians for which the collection site and the larval environment are known to be the same will be critical for identifying putative drivers (e.g., temperature, pH, dissolved oxygen, organic chemicals, etc.) of sex reversal as well as contributions to cohort sex ratio variation.

Intersex

While we identified both genotypically female and male intersex frogs, intersex and sexual genotype-phenotype discordance frequencies were not correlated, perhaps suggesting the mechanisms underlying intersex and sex reversal may be different. By extension, this suggests that we cannot conflate intersex and sex reversal as has been suggested previously (Murphy et al., 2006). The few experiments which have assessed both intersex and genotypic sex have not reported relationships between sex reversal and intersex (Hayes et al., 2010; Tamschick et al., 2016a, 2016b). Future experimental and field studies should assess any relationship between intersex and sex reversal to better understand the conditions and mechanisms underlying these processes. We note that intersex is likely a terminal phenotype in most amphibians and, at least in green frogs, is not evidence of sequential hermaphroditism which currently is not known to occur in amphibians. Phenotypic males with intersex testes may represent two developmental scenarios. In one pathway, larval genotypic females either partially or completely develop their ovaries but are later environmentally-induced to sex reverse, degenerating their ovaries and developing testes prior to metamorphosis. This pathway has been experimentally demonstrated in green frog tadpoles previously (Mintz, Foote & Witschi, 1945). Alternatively, genotypic males may experience abnormal gonadal differentiation but do not flip their sexual trajectories toward being phenotypic females (Tamschick et al., 2016a, 2016b). Both field-based and experimental research could provide exciting insight into the relative contribution of these two developmental pathways to intersex frequencies in the wild.

In recent years, discussions about amphibian sexual development have often been guided by striking results from laboratory-based ecotoxicological studies which have led to the impression that sex reversal and intersex are predominantly the result of human factors (Hayes et al., 2002, 2010; Pettersson & Berg, 2007; Tamschick et al., 2016a, 2016b). However, previous work in other frog species throughout the past century (Witschi, 1914, 1929, 1930; Hsu, Yu & Liang, 1971) including a field survey in Africa (Du Preez et al., 2009), also indicates that intersex can be widespread among wild populations and may not be necessarily related to contaminants but rather natural environmental factors like temperature. It currently remains unclear what environmental conditions cause intersex to develop in green frogs and other taxa. Our findings that sexual genotype-phenotype discordance and intersex frequencies were not associated with land use is therefore intriguing and perhaps surprising. These results should prompt a reconsideration of the perception that these sexual characteristics are predominantly of anthropogenic origin.

Evolutionary and ecological implications

Despite having experienced at least 32 transitions between GSD (e.g., XY vs ZW sex chromosome) systems (Evans, Alexander Pyron & Wiens, 2012), amphibians are not generally believed to naturally undergo sex reversal (but see Perrin, 2009; Alho, Matsuba & Merila, 2010; Rodrigues et al., 2018). Yet, most mechanistic models of sex determination transitions invoke a role for environmentally-sensitive sex determination (Bull, 1981; Grossen, Neuenschwander & Perrin, 2011; Quinn et al., 2011; Schwanz et al., 2013; Muralidhar & Veller, 2018). Our findings in conjunction with research on European common frogs should encourage future investigation in to whether sex reversal influences transitions among amphibian SD modes.

Numerical modeling suggests that the frequency of sexually-discordant genotypic females (XX♂) reported here should have minimal effects on population size or viability unless sexually-discordant genotypic females (XX♂) have reduced fitness (Cotton & Wedekind, 2009; Wedekind, 2010). The low frequency of sexually-discordant genotypic males (XY♀) in green frogs, however, may permit the type of sex-chromosome rejuvenation suggested by the “fountain of youth” hypothesis (Perrin, 2009; Grossen, Neuenschwander & Perrin, 2012).

Sex chromosome recombination

The main complication with identifying sex reversal in taxa like amphibians, fishes, and non-avian reptiles is having to identify and rely on sex-linked genetic loci which may or may not be the SD locus. While “the SD locus” should experience no recombination, other loci on the sex chromosomes will experience varying degrees of recombination depending on their proximity to the SD locus as well distance between the loci and the centromere. Therefore, while a number of loci may all be sex-linked based on their corresponding location to SD locus, certain individual sex-linked loci will show different degrees of recombination and therefore sex-linkage. Because of this heterogeneity in linkage, alleles associated with the heterogametic sex chromosome (e.g., Y or W) will possibly be observed in the opposite genotypic sex and the heterogametic sex will occasionally exhibit the homogametic sex’s genotype at certain sex-linked loci (Hussain et al., 1994). Recombination rates for sex-linked loci can be directly estimated with family inheritance data from both parents and their offspring (Matsuba, Miura & Merila, 2008). While having directly estimated recombination rates would be useful to our study, we are not prevented from identifying informative sex-linked loci in their absence (Lambert, Skelly & Ezaz, 2016; Brelsford et al., 2017).

Bayesian probabilistic modeling approaches like the one used here and first used by Alho, Matsuba & Merila (2010) to infer sex reversal in European common frogs are useful tools for inferring the genotypic sexes of individuals across multiple loci by accounting for variation in sex-linkage and, by extension, recombination. As we discussed previously (Lambert, Skelly & Ezaz, 2016), the sex-linked loci identified for green frogs show varying degrees of sex-linkage, likely due to recombination. RaclCT001 is a particularly interesting locus because its alleles were perfectly sex-linked (all females were homozygous for the common X-allele and males were heterozygous) across 23 adult females and 54 adult males in the original dataset (Lambert, Skelly & Ezaz, 2016). And in this current study of over 400 frogs, RaclCT001 was again perfectly sex-linked outside of inferred sex-reversed individuals suggesting that it is either at or very near the SD locus in these populations and experiences minimal if any recombination. By analyzing across multiple loci, and in particular locus RaclCT001, we can accurately assign the genotypic sex to green frogs in our study population even if multiple loci for an individual are apparently discordant. Although we cannot fully discount the possibility that sex chromosome recombination influenced our results, the methods we employed account for variation in sex-linkage and allow for useful inference into individuals’ genotypic sex, and therefore sex reversal, in the absence of a well-established SD locus. Recombination likely accounted for variation in sex-linkage at several of the loci (RaclCT005, -007, and -009) used here but was unlikely to be so extreme to as result in our analyses mis-assigning genotypic sexes. Our analyses using either the complete sex-linked markers set or only the most tightly sex-linked loci produced similar results, validating our interpretations of sex reversal. However, future work in this system using family data will be useful for unequivocally disentangling the effects of recombination in these markers, if there is any.