Introduction

One of the central questions in sex chromosome research is whether the sex chromosomes found in different species are homologous (Ohno, 1967; Miura, 2007; Graves and Peichel, 2010; Pease and Hahn, 2012; O’Meally et al., 2012). Despite recent progress (see, for example, Gamble et al., 2014; Gamble and Zarkower, 2014; Chen et al., 2014; Nanda et al., 2014; Srikulnath et al., 2014), identifying the homomorphic sex chromosomes found in many cold-blooded vertebrates, invertebrates and some plants remains a major challenge. Researchers are currently interested in the factors that maintain homomorphy in these sex chromosome systems. Two key hypotheses have been put forward: either sex chromosome turnovers occur too frequently to allow substantial degeneration of the nonrecombining Y (or W) chromosome, or a low level of recombination is maintained between the X and Y (or Z and W) chromosomes in these systems (Perrin, 2009; Beukeboom and Perrin, 2014). Distinguishing between these two hypotheses requires us to quantify rates of sex chromosome turnover across the phylogeny, and this can only be done by identifying the sex chromosomes in a large number of species. Here, we propose a novel approach to identify sex chromosomes through high-density linkage mapping. Importantly, the siblings used to prepare this linkage map do not need to be sexed.

Shortly after the onset of reduced recombination between a proto-X and proto-Y chromosome (Bachtrog, 2013), we expect divergence in allele frequencies between the two gametologs. In this case, an X and Y chromosome sampled from a population are expected to have more differences between them than two X chromosomes. In a full-sib mapping cross, more markers will be heterozygous and therefore map informative in the XY male parent than in the XX female. In contrast, the number of heterozygous markers on each autosome is expected to be similar between the two parents. In a comparison of the two sex-specific maps, if one linkage group has a strong difference in marker number between the male and female maps, that linkage group is likely to correspond to the sex chromosome.

We demonstrate this method in the European tree frog Hyla arborea. Previous work on this species has shown that its sex determination system is at least 5 million years old, but that occasional X–Y recombination prevents degeneration of the Y chromosome (Stöck et al., 2011). The sex chromosome is largely homologous to chromosome 1 in Xenopus tropicalis, and contains the candidate sex determination gene Dmrt1 (Brelsford et al., 2013). Males and females differ strongly in recombination rate (Berset-Brändli et al., 2008; Dufresnes et al., 2014a, b), but X–Y recombination is higher and X–Y differentiation is lower in Balkan glacial refugia than in recently colonized northwestern Europe (Dufresnes et al., 2014b). In this study, we validate our method for identifying sex chromosomes by showing that only linkage group 1, homologous to X. tropicalis chromosome 1, has a strongly different number of markers between paternal and maternal linkage maps. This approach also enables us to localize sex differences in recombination to specific parts of the genome, and extends to the entire genome our previous work on conserved synteny in anurans that until now was largely limited to a single chromosome.

Materials and methods

Samples

We constructed a linkage map based on 90 offspring from a single H. arborea cross collected in Krk, Croatia (family K10 in Dufresnes et al., 2014b). DNA was collected from both parents using buccal swabs. Fertilized eggs were kept until shortly after hatching, when tadpoles were preserved in ethanol. DNA from parents and offspring was extracted using a BioSprint robotic workstation (Qiagen, Venlo, Netherlands).

Genotyping by sequencing

We collected genotyping-by-sequencing data using a novel molecular protocol adapted from Parchman et al. (2012), Peterson et al. (2012) and Purcell et al. (2014). The complete protocol is available in Supplementary Materials. Briefly, we digested DNA with restriction enzymes SbfI and MseI, ligated barcoded adapters, amplified each individual sample in four separate PCR reactions, pooled all PCR products and selected fragments between 400 and 500 bp using agarose gel. To increase sequence depth for the two parental DNA samples, PCR reactions for these samples were run at triple the volume of offspring samples. The library of 92 samples was sequenced on a single Illumina (San Diego, CA, USA) HiSeq 2000 lane by the Lausanne Genetic Technology Facility (Lausanne, Switzerland), with single-end 100 bp reads.

Genotype calling and filtering

We demultiplexed the raw reads using the process_radtags module of Stacks (Catchen et al., 2013), removing low-quality reads, reads with uncalled bases and reads that failed the Illumina ‘chastity’ filter. Adapters were trimmed from demultiplexed reads using a custom shell script. Reads were then aligned to the H. arborea low-coverage draft genome (Dryad, doi: 10.5061/dryad.n856c) using Bowtie2 (Langmead and Salzberg, 2012). Variants were called using the mpileup command of Samtools (Li et al., 2009), using only reads with mapping quality of at least 20. We used VCFtools (Danecek et al., 2011) to filter the resulting raw variants. First, genotypes with a quality score <20 were converted to missing data. Loci with >20% missing data, offspring heterozygosity >0.8 or a minor allele frequency <0.15 were removed.

We then generated separate lists of paternal-informative and maternal-informative markers, excluding loci that were heterozygous in both parents or in neither parent. For example, a male-informative marker with two alleles, A and G, could have genotypes AG in the male and GG in the female, with offspring expected to show the AG and GG genotypes in equal numbers. A female-informative marker could have genotypes GG in the male and AG in the female, offspring showing the same pattern as in the previous example. A marker with AG genotypes in both parents would be excluded, as in heterozygous offspring we would be unable to determine which parent provided the A allele and which provided the G allele. Because the H. arborea reference genome is highly fragmented (N50, 1.2 kbp), we retained only the marker with the lowest proportion of missing data from each genome scaffold for each of the two sex-specific data sets, assuming no recombination within the short genome fragments. Finally, we searched for and corrected Mendelian segregation errors. For example, if 45 offspring were homozygous for the C allele at a locus, 44 heterozygous for C and T and one homozygous for T, the rare homozygous genotype was corrected to a heterozygous call. Markers with >5% of offspring having the homozygous genotype unexpected based on parent genotypes were removed. For markers with unexpected homozygous genotypes at a frequency <5%, these were corrected to heterozygous genotypes. Genotypes were then exported in 012 format using VCFtools.

Linkage mapping

We constructed sex-specific genetic maps using MSTmap (Wu et al., 2008), using cross type ‘DH.’ Following Gadau et al. (2001), to account for the unknown phase of markers in each parent we added each marker to the data set twice, once with each possible phase. Each linkage group then appeared in the output twice with opposite marker phases; duplicate linkage groups were removed manually.

Assessment of conserved synteny with Xenopus tropicalis

We searched the X. tropicalis genome sequence (available at ftp://ftp.xenbase.org/pub/Genomics/JGI/Xentr7.1/) for each H. arborea scaffold that contained a marker on the linkage map using blastn. Blast hits were retained only if the e-value of the best hit was five orders of magnitude better than that of the second hit (Purcell et al., 2014). Homology between the X. tropicalis genome and H. arborea linkage map was visualized using an Oxford grid plot.

Results

Genotyping by sequencing

We obtained 153 million demultiplexed, quality-filtered Illumina reads, with number of reads per individual ranging from 1.1 to 1.9 M in offspring, and 5.8 and 6.5 M in the parents. Approximately 21% of these reads mapped uniquely to the H. arborea genome draft; this low rate is likely because of the incompleteness of the draft genome assembly (1.2 Gbp of an expected 5 Gbp). After genotype calling and filtering, we retained 1608 male-informative and 1461 female-informative single-nucleotide polymorphisms or small indel markers, with an average sequence depth of 16.5.

Linkage mapping

We recovered 12 linkage groups for both sex-specific maps, consistent with the species karyotype (both maps available in Supplementary Materials). The number of informative markers per linkage group was similar between male and female maps for 11 of the 12 linkage groups. In contrast, linkage group 1, the sex chromosome pair, contained nearly three times as many male-informative as female-informative markers (Figure 1).

Figure 1
figure 1

Linkage group 1, the sex chromosome in H. arborea, contains nearly three times as many informative markers on the male-specific map compared with the female-specific map. The 11 autosomal linkage groups fall close to the dashed line, illustrating equal marker numbers between the sexes.

Full size image

A subset (10.7%) of mapped markers could be assigned positions in the X. tropicalis genome sequence (Table 1). These markers revealed strongly conserved synteny between Hyla and Xenopus, with most Hyla linkage groups matching a single Xenopus chromosome (Figure 2). Because no data are available to associate linkage groups with the Hyla karyotype, we named the Hyla linkage groups based on their homologous chromosomes in X. tropicalis. In several cases, even the ordering of markers along the chromosome is largely conserved across the entire chromosome (for example, Hyla linkage groups 5, 6 and 9) or large chromosomal regions (for example, Hyla linkage groups 2 and 3). Marker order is conserved on the linkage group 1, with the exception of a small region. Brelsford et al. (2013) noted that a region of X. tropicalis chromosome 1 corresponding to avian chromosome 26 was not found on the same chromosome in H. arborea; in this study, we map that rearranged region to linkage group 2.

Table 1 Number of mapped markers, and number available for comparisons between sex-specific Hyla arborea maps and between H. arborea and Xenopus tropicalis
Full size table
Figure 2
figure 2

Oxford grid comparison of H. arborea linkage map and X. tropicalis genome shows blocks of conserved synteny, often extending across entire chromosomes. Marker order is largely conserved on linkage group 1 that is sex linked in several frog species but autosomal in X. tropicalis.

Full size image

Most male linkage groups contained a cluster of markers near the center of the linkage group with an identical segregation pattern, suggesting a large region of suppressed recombination (Figure 3a). This pattern was most pronounced in the larger linkage groups. In contrast, the number of markers per cM was close to uniform across the female map (Figure 3b).

Figure 3
figure 3

Marker density is very high at the centers of male linkage groups (a), and approximately uniform across female linkage groups (b). Linkage groups containing more markers also show a larger difference in map length between the male and female maps. Taken together, this suggests that recombination in males occurs only near the ends of large chromosomes.

Full size image

Discussion

Linkage maps provide an opportunity to investigate recombination and serve as a resource for research on genome structure, dynamics and rearrangements. Our female- and male-specific maps respectively contain 1450 and 1598 markers, with an average intermarker distance of 1.65 and 0.72 cM. This marker density compares favorably with other published amphibian maps; to our knowledge, only the model species X. tropicalis has a higher marker density (2886 markers, intermarker distance 0.38 cM; Wells et al., 2011). Availability of a high-density linkage map will enable future research on how population genetic patterns vary across different regions of the genome, essential for understanding the genomics of sex chromosome evolution and speciation. Immediately, it also demonstrates the potential of linkage mapping for identification of sex chromosomes without knowledge of offspring sexes. The sex-specific linkage maps we present here show that the previously documented absence of recombination in males is restricted to the central parts of the chromosome. In contrast, chromosome ends have similar rates of recombination between the two sexes. In addition, we identify a high level of synteny between the whole genome of Hyla and the distantly related model frog species X. tropicalis.

A new method for identifying sex chromosomes

Using high-density linkage maps, we demonstrate that we can identify the sex chromosome from families with sex-unknown offspring. When paternal and maternal chromosomes are drawn from the same population, the number of heterozygous single-nucleotide polymorphisms per chromosome pair is expected to be equal for the father and mother. If the X and Y chromosomes show significant genetic differentiation, heterozygosity will be higher in the male than the female for that chromosome only. By plotting male-informative markers against female-informative markers for each of 12 linkage groups, we show that linkage group 1 is an outlier with nearly three times as many male-information as female-informative markers, consistent with previous studies (Berset-Brändli et al., 2008; Brelsford et al., 2013) identifying this as the sex chromosome.

Our approach will serve as a powerful tool to identify the sex chromosome of species with large numbers of offspring available, but for which juvenile sex identification is difficult or impossible. It does, however, have several limitations. First, it requires significant differentiation in allele frequencies across a large region of the sex chromosome. The important parameter for detecting the sex chromosome is the ratio of differentiation between the X and Y chromosomes (DXY) to polymorphism on the X (πX). In the parents of this cross, our estimated DXY/πX is 2.9. Among the 11 autosomal linkage groups, LG 9 has the largest deviation from equal marker number between the sexes, with 1.7 × more female-informative than male-informative markers. The method’s power will be highly dependent on polymorphism and linkage disequilibrium in the parental source population as well as the level of differentiation between the X and Y. However, to the extent that other species are similar to this H. arborea population, we suggest that a DXY/πX ratio of >2 will be necessary to detect the sex chromosome. Second, it requires identifying X-linked and Y-linked markers as different alleles of the same locus, and not different loci. If a Y chromosome is missing large regions that are present on the X, or if sequence divergence is too high to recognize orthologous X and Y sequence reads, the number of male-informative markers on the sex chromosome will decline, reducing the method’s power. Finally, a large autosomal supergene with two differentiated haplogroups (see, for example, Purcell et al., 2014) could produce a false positive signal, if one parent of the mapping cross is heterozygous for the inversion and the other is homozygous. Further testing of our method will be necessary to determine how often its requirements are met in other species.

Synteny between Hyla and Xenopus

Chromosomes are largely conserved between H. arborea and X. tropicalis, despite divergence time estimates ranging from 165 to 305 million years (Aris-Brosou and Yang, 2002; Wiens, 2007). A previous study demonstrated conserved synteny across H. arborea, Rana temporaria, Bufo siculus and X. tropicalis at chromosome 1 (Brelsford et al., 2013), with one small segment of the X. tropicalis chromosome found on a different chromosome in H. arborea. In H. arborea, Dufresnes et al. (2014a) inferred two other linkage groups that each contained two markers mapped to the same X. tropicalis chromosome (chromosomes 2 and 8), hinting at conserved synteny in at least two other chromosomes. More recently, Sun et al. (2015) reported conservation of large synteny blocks between X. tropicalis and the Tibetan frog Nanorana parkeri, although the Nanorana genome assembly did not allow comparison of entire chromosomes. We extend these previously published results by assigning each Hyla linkage group to one Xenopus chromosome (or in the case of LG 4A7A, two chromosomes), and by identifying large blocks of conserved marker order on several chromosomes (Figure 2a). Marker order is similar in four of the 12 Hyla chromosomes (1, 5, 6 and 9), and a single large-scale inversion is apparent on two others (2 and 3). Of the 10 X. tropicalis chromosomes, three (4, 7 and 8) are split in two, with 4A and 7A reunited into a new Hyla chromosome.

A similar level of synteny has been observed in birds based on comparison of chicken and zebra finch (Stapley et al., 2008), representing close to 100 Myr of divergence. In contrast, some systems have substantial levels of interchromosomal rearrangement even with much shorter divergence times, as in ants (Purcell et al., 2014). Interestingly, the Z chromosome in birds has undergone more intrachromosomal rearrangement than autosomes (Griffin et al., 2006; Stapley et al., 2008). The X. tropicalis chromosome 1 is partially homologous to the avian Z and is recurrently used as the sex chromosome in multiple frog species (Brelsford et al., 2013). Thus, it is surprising to find little structural variation in chromosome 1 between Xenopus and Hyla.

Differences in male and female recombination rate

Heterochiasmy, or difference between sexes in recombination rate, is a well-documented phenomenon in many taxa. Although Haldane (1922) and Huxley (1928) noted that recombination tends to be lower in the heterogametic sex, many exceptions have been observed and comprehensive explanations for this phenomenon are still lacking (Mank, 2009). In frogs, previous low-density linkage maps documented a strong asymmetry in recombination between males and females, with males having greatly reduced recombination (see, for example, Nishioka and Sumida, 1994; Sumida and Nishioka, 1994; Berset-Brändli et al., 2008; Cano et al., 2011; Rodrigues et al., 2013). Our high-density linkage map allows us to localize the recombination events along each linkage group. The differences in map length between the male and female are particularly striking. In females, the map length is roughly proportional to the marker number, suggesting that recombination rate scales linearly with chromosome length. On the other hand, the male map lengths of all linkage groups are close to 100 cM, implying approximately one crossover per chromosome regardless of the physical size of the chromosome. It should be noted that these calculations are based on single-nucleotide polymorphism density per cM. This measure is expected to be only a rough proxy for the physical length per cM, as genotyping-by-sequencing markers and heterozygosity may not be uniformly distributed across the genome. Our finding that male recombination is restricted to chromosome ends, whereas female recombination is relatively constant across the full length of each linkage group, suggests that chromosome ends experience more recombination and genes in these regions may evolve differently (Cutter and Payseur, 2013). Our localization of heterochiasmy is fully consistent with previous cytogenetic findings. In female frogs, chiasmata have been shown to occur in proportion to chromosome size (being more frequent in large bivalents than in small ones) and uniformly along the chromosomes. Male frogs, in contrast, typically show two chiasmata per meiosis (independent of bivalent size), and always localized in terminal regions, conferring a characteristic ring shape to chromosomes during diakinesis (see, for example, Morescalchi and Galgano, 1973).

Most frogs with known systems of sex determination are male heterogametic. In Hyla, Bufo and Rana, recombination is reduced in males (Berset-Brändli et al., 2008; Stöck et al., 2013; Rodrigues et al., 2013), consistent with the Haldane–Huxley rule. However, we still do not know whether one of these patterns is the cause of the other. Insights into the basis of these patterns will likely emerge from exceptional systems. For instance, investigating recombination rates in frogs with ZW sex determination will help us to disentangle the causal relationships between sex determination and heterochiasmy.

Conclusions

Overall, we have used our high-density linkage maps to achieve three goals. First, we have evaluated a novel method to infer the sex chromosome from the comparison of sex-specific maps. Second, we have demonstrated conserved synteny between the whole genomes of Hyla and Xenopus. Finally, we have identified the relative positions and frequency of recombination events in male and female H. arborea. These maps will also serve as vital tools for future population genomic work.

Data archiving

Demultiplexed genotyping-by-sequencing reads and whole-genome reads are available on NCBI Sequence Read Archive (BioProjects PRJNA292823 and PRJNA292815). The draft low-coverage genome assembly is available on Dryad (DOI 10.5061/dryad.n856c).