Sex-linked markers in the North American green frog (Rana clamitans) developed using DArTseq provide early insight into sex chromosome evolution

Specimen selection

The range of R. clamitans occupies much of eastern North America, extending from north Florida into southeastern Canada (Fig. 1a); our sampling was located in Connecticut, a state in the northeastern United States (Fig. 1a, b). We used two developmental stages: tadpoles reared to sexual differentiation in controlled-environment mesocosms and wild-caught adults. We predict that tadpoles reared in carefully constructed mesocosms should experience minimal sex reversal, maximizing the probability of genotype-phenotype concordance [30]. We collected portions of freshly-laid clutches from three ponds (1 from Septic7, 3 from Forest5, and 2 from Forest6; Fig. 1b), hatched the embryos in the lab at 19 °C on a 12L: 12D cycle, and reared tadpoles in outdoor mesocosms. We did not use any leaf litter to our mesocosms as leaf leachate has been shown to influence sex steroid pathways [31]. Tadpoles of this species can often be sexed within a few months of hatching despite taking over a year to metamorphose [32]. At ~4 months old, we euthanized tadpoles, sexed them under a dissecting microscope, and collected tail muscle tissue in 90 % ethanol for genetic analysis. Not all tadpoles were fully sexually differentiated at this point and so we focused only on tadpoles with clearly differentiated testes and ovaries (Fig. 2). While not random, we selected tadpoles that had clearly sexually differentiated under controlled conditions. We used 1–3 tadpoles of each phenotypic sex from each clutch for analyses. For adults, we capitalized on Yale Peabody Museum tissues (Additional file 1: Table S1) preserved from wild-caught adults in south-central Connecticut in 2014 as well as an additional five adults from Sewer3 pond collected in 2015. These individuals made up the majority of our specimens and represented variation in local geography and human land use (Fig. 1b). While wild-caught adults may have experienced sex reversal, our goal was to test whether DArTseq could identify sex-linked markers despite the presence of possible sex-reversed individuals and whether we could identify sex-reversed individuals as outliers in our data. We used 1–5 females and 3–11 males from each pond. We used 17 samples in duplicate for sequencing quality control. The source ponds for all specimens ranged from less than 1 km to over 52 km apart (Table 1).

	Septic7	Forest5	Septic4	Forest3	Sewer3	WNHL	YP10	Forest6
Septic7	0	36,081	44,753	31,080	52,290	3519	30,836	31,184
Forest5	36,081	0	43,348	22,776	48,315	32,851	30,782	31,650
Septic4	44,753	43,348	0	21,181	7555	44,647	15,273	14,592
Forest3	31,080	22,776	21,181	0	27,224	29,533	8266	9102
Sewer3	52,290	48,315	7555	27,224	0	52,126	22,381	21,779
WNHL	3519	32,851	44,647	29,533	52,126	0	30,197	30,633
YP10	30,836	30,782	15,273	8266	22,381	30,197	0	871
Forest6	31,184	31,650	14,592	9102	21,779	30,633	871	0

Genotyping-by-sequencing (GBS)

DNA extractions and sequencing were performed using DArTseq™ (Diversity Arrays Technology Pty Ltd). DArTseq™ represents a combination of a DArT complexity reduction methods and next generation sequencing platforms [27]. Therefore, DArTseq™ represents a new implementation of sequencing of complexity reduced representations [33] and more recent applications of this concept on the next generation sequencing platforms [28, 29], conceptually similar to double digest RADseq. A detailed description of the DArTseq™ methodology can be found in Kilian et al. [27]. While RADseq typically yields markers 85 bp or longer (e.g., [21, 23, 24]), DArTseq™’s criteria (see below) often produce sequences 69 bp or shorter in length. Compared to other similar approaches, DArTseq™ yields a lower density of markers (10’s of thousands and up to 35,000 loci versus >800,000 loci with a GBS approach) but has substantially higher coverage and results in less missing data [34]. In addition to these advantages of DArTseq™, we chose this platform because it can directly score samples as heterozygous/homozygous at each locus with the lower density approach [34] and due to the ease of optimizing the platform for each species as well as its capacity to produce thousands of short, high quality polymorphic loci using a custom analytical pipeline [35–37]. Four methods of complexity reduction were tested in R. clamitans (data not presented) and double digestions with PstI-SphI method was selected.

DNA samples were processed in digestion/ligation reactions principally as per Kilian et al. [27] but replacing a single PstI-compatible adaptor with two different adaptors corresponding to two different Restriction Enzyme (RE) overhangs. The PstI-compatible adapter was designed to include Illumina flowcell attachment sequence, sequencing primer sequence and “staggered”, varying length barcode region, similar to the sequence reported by Elshire et al. [29]. The reverse adapter contained flowcell attachment region and SphI-compatible overhang sequence. Only “mixed fragments” (PstI-SphI) were effectively amplified in 30 rounds of PCR using the following reaction conditions: 94̊ C for 1 min then 30 cycles of 94 °C for 20 s, 58 °C for 30 s, 72 °C for 45 s, and 72 °C for 7 min. After PCR, equimolar amounts of amplification products from each sample of the 96-well microtiter plate were bulked and applied to c-Bot (Illumina) bridge PCR followed by sequencing on Illumina Hiseq2500. The sequencing (single read) was run for 77 cycles.

Sequences generated from each lane were processed using proprietary DArT analytical pipelines (PL). In the primary pipeline the Fastq files were first processed to filter away poor quality sequences, such as those with reproducibility below 90 %, read depth lower than 3.5 for SNPs or 5 for presence-absence markers, and applying more stringent selection criteria to the barcode region compared to the rest of the sequence. In that way the assignments of the sequences to specific samples carried in the “barcode split” step were very reliable. No samples were dropped due to low coverage across loci but individual sequences were removed if they did not meet the above criteria. Approximately 2,500,000 sequences per barcode/sample were identified and used in marker calling. The average read depth across loci was 9.2 reads per individual per locus for reference alleles and 6.5 for SNP alleles. Finally, identical sequences were collapsed into “fastqcoll files”. The fastqcoll files were “groomed” using DArT PL’s proprietary algorithm which corrects low quality base from singleton tag into a correct base using collapsed tags with multiple members as a template. The “groomed” fastqcoll files were used in the secondary pipeline for DArT PL’s proprietary single nucleotide polymorphism (SNP) and SilicoDArT (presence/absence of restriction fragments in representation; PA markers) calling algorithms (DArTsoft14). For SNP calling all tags from all libraries included in the DArTsoft14 analysis are clustered using DArT PL’s C++ algorithm at the threshold distance of 3, followed by parsing of the clusters into separate SNP loci using a range of technical parameters, especially the balance of read counts for the allelic pairs. Additional selection criteria were added to the algorithm based on analysis of approximately 1000 controlled cross populations. Testing for deviations from Hardy-Weinberg equilibrium of alleles in these populations facilitated selection of technical parameters to effectively discriminate true allelic variants from paralogous sequences. In addition, multiple samples were processed from DNA to allelic calls as technical replicates and scoring consistency was used as the main selection criteria for high quality/low error rate markers. Calling quality was assured by high average read depth per locus (average across all markers was over 30 reads/locus).

Marker selection

Our goal was to isolate single nucleotide polymorphisms (SNP) and presence-absence (PA) markers that are sex-linked in R. clamitans. For SNP loci, we define the “reference allele” as the allele that was sequenced most often in our data. For sex-linked markers in an XX-XY sex determining system, reference alleles will be those found on the X-chromosome. Here, “SNP alleles” are those which show polymorphisms relative to the reference allele. In an XX-XY system, SNP alleles should be alleles associated with the Y-chromosome and, as a corollary, should be on or near the male-determining region if the allele is tightly Y-specific. Higher degrees of recombination may be expected between sex chromosomes in endothermic species with homomorphic sex chromosomes [14, 25, 38]. In an XX-XY system if the two sex chromosomes recombine, SNP alleles should occasionally appear on the X chromosome, particularly in more distal regions of the sex chromosomes [25, 39]. If this occurs, then some males might be homozygous for SNP alleles at particular loci. Because of this, females could also be heterozygous, exhibiting a copy of the SNP allele. The probability of a female being homozygous for a SNP allele is low.

A previous study on inheritance patterns of an allozyme suggested an XX-XY system in R. clamitans [40]. While we assumed an XX-XY system, we also searched for loci showing a ZZ-ZW pattern. For evaluating loci associated with an XX-XY system, we targeted SNP loci that maximized female homozygosity at the reference allele, maximized male heterozygosity, and which minimized homozygosity for the SNP allele in both sexes, but particularly for females. Specifically, we kept loci where frequencies of female homozygosity for the reference allele were at least 80 %, homozygosity at the SNP allele were at most 10 %, and heterozygosity were no more than 20 %. For males, we kept loci where frequencies of homozygosity for the reference allele were at most 10 % and heterozygosity were at least 75 %. This allowed sex-linked markers to show higher degrees of SNP allele homozygosity for males if recombination were occurring. We also only included SNP loci which were sequenced for at least 90 % of each sex. For PA markers we selected loci which had restriction fragments sequenced in at least 90 % of one sex and not sequenced in at least 90 % of the other sex. For both SNP and PA markers, we analyzed tadpole and adults data separately and compared results between these groups afterwards. We took a similar but opposite approach for targeting loci with a ZZ-ZW system.

To assess the genetic association between each locus and phenotypic sex, we performed Cochran-Armitage tests on the adult samples with the “independence_test” function in the R package ‘coin’. We also used Hamming Distance matrices with a custom R script to estimate how well the combined markers perform in predicting phenotypic sex. Hamming Distance calculates the number of pairwise differences between all individuals across all loci. For all sex-linked loci which met our criteria and which had statistically significant associations with phenotypic sex, we used NCBI BLAST [41] to discover homologies of sex-linked SNPs markers. Specifically, we searched the “nr” database with the default settings in NCBI’s blastn suite. We searched using both the REF and SNP alleles and only report blastn hits related to sexual development or sex determination and those with e-values of 1.5 or below.

Estimating random sex-linkage

Here P is the probability that a given locus, i, is sex-linked. The probability that either a female is homozygous at a given locus or a male or heterozygous is 0.5. The number of individuals sequenced at the locus is n. After sequencing, we multiplied P by the number of high-quality SNPs produced to estimate how many SNPs were expected to exhibit a sex-linked pattern by chance. Using this calculation, we estimated the estimated number of expected sex-linked markers produced by chance across a range of sample sizes and loci for comparisons with other studies. We also estimated the number of expected random sex-linked loci for both the tadpole samples as well as the adult samples used here.