A dense linkage map for a large repetitive genome: discovery of the sex-determining region in hybridizing fire-bellied toads (Bombina bombina and Bombina variegata)

Abstract

Genomic analysis of hybrid zones offers unique insights into emerging reproductive isolation and the dynamics of introgression. Because hybrid genomes consist of blocks inherited from one or the other parental taxon, linkage information is essential. In most cases, the spectrum of local ancestry tracts can be efficiently uncovered from dense linkage maps. Here, we report the development of such a map for the hybridizing toads, Bombina bombina and Bombina variegata (Anura: Bombinatoridae). Faced with the challenge of a large (7-10 Gb), repetitive genome, we set out to identify a large number of Mendelian markers in the nonrepetitive portion of the genome that report B. bombina vs B. variegata ancestry with appropriately quantified statistical support. Bait sequences for targeted enrichment were selected from a draft genome assembly, after filtering highly repetitive sequences. We developed a novel approach to infer the most likely diplotype per sample and locus from the raw read mapping data, which is robust to over-merging and obviates arbitrary filtering thresholds. Validation of the resulting map with 4755 markers underscored the large-scale synteny between Bombina and Xenopus tropicalis. By assessing the sex of late-stage F2 tadpoles from histological sections, we identified the sex-determining region in the Bombina genome to 7 cM on LG5, which is homologous to X. tropicalis chromosome 5, and inferred male heterogamety. Interestingly, chromosome 5 has been repeatedly recruited as a sex chromosome in anurans with XY sex determination.

Keywords: anurans; genome assembly; hybrid zone; large-scale synteny; linkage map; population pileups diplotyping; segregation distortion; sex-determining region; targeted capture.

Figure 1

Overview of B. variegata genome assemblies. The large circle represents the total read set. Assemblies are numbered 1–4. hi K, highly over-represented kmers (at least 100x average frequency); Tr, B. v. variegata transcriptome; repeat(−) reads, repeat-subtracted portion of the total read set. CLC contigs with a BLAST+ hit (query: transcriptome) are represented by a blue rectangle.

Figure 2

Polarization of the raw read coverage. The top two plots show the raw read coverage along the reference sequence of locus 332172 for F0 B. variegata (A) and F0 B. bombina (B). The B. variegata individual is homozygous for the reference state (R) at all sequence positions, whereas the B. bombina individual has a number of variants. Two homozygous (156 and 343) and two heterozygous (110 and 224) variant positions are highlighted. For these four, the matrix entries are listed below the plots. A polarized matrix, M_p, is computed from these read counts (see text, C), in which sequence states associated with B. variegata have positive entries and sequence states associated with B. bombina have negative entries. For each sample, raw read counts are then multiplied by M_p. Average positive entries and average negative entries result in a B. bombina score and a B. variegata score, respectively, and when plotted in a coordinate system (D), samples can be assigned to three clusters representing BbHOM, HET, and BvHOM. Note that the heterozygous variants (B) do not interfere with the clustering into three diplotypes.

Figure 3

The distribution of repeat types. We show the 200 REPdenovo contigs with the highest copy number. Transposable element orders represented by more than 10 contigs in this set are identified by color. The classification follows Wicker et al. (2007). Contigs without a match in Repbase (blastn and tblastx) are labeled as “no match” and ordered separately. LTR, long terminal repeat retrotransposon; DIRS, Dictyostelium intermediate repeat sequence; TIR, terminal inverted repeat DNA transposon.

Figure 4

The Bombina linkage map. The linkage map was visualized with LinkageMapView (v. 2.1.2) (Ouellette et al. 2018). Horizontal bars represent marker loci. Colors indicate marker density in cM/locus from 0.2 (red) to 2.1 (blue). This figure is best viewed at maximum magnification.

Figure 5

Segregation distortion, χ², by family and linkage map position. Dashed horizontal lines are significance thresholds: the lower line is the Bonferroni correction based on the number of chromosome arms, and the upper line is the critical value for the Benjamini and Hochberg false discovery rate (the experiment-wise α is 0.05 for both). For each significant spike, which is indicated with an arrowhead, the genotype showing the strongest deviation is noted along with a (+) or (−) label, where (+) = excess and (−) = deficit. Groups of spikes with the same deviation are separated by vertical lines above the plot. For clarity, 22 observations from 20 loci with χ² > 20 are excluded from the plot.

Figure 6

Synteny between B. variegata and X. tropicalis. Circos (v0.69-6) (Krzywinski et al. 2009) plot of 737 B. variegata target sequences from the 12 LGs (Bv, unit is cM) aligned against the X. tropicalis genome assembly (Xt, unit is Mb) with BLAST+ (v. 2.9.0) (Camacho et al. 2009).

Figure 7

Estimated frequency of two sex-diplotype combinations among homozygous F2 individuals, b. See text for the definition of b. The global minimum on LG5 indicates the sex-determining region. The blue line represents the null hypothesis of b = 0.5.

Figure 8

Raw read coverage at locus 5568 in the F0 generation. Plots follow the format of Figure 2: raw read coverage (y-axis) is shown for each reference position (x-axis) and sequence state (R = reference state) for F0 B. variegata (A) and F0 B. bombina (B). Nonreference states are highlighted in color. Variants in the sex-linked haplotype (A) are connected with a dashed line.