The Genomic Architecture of a Rapid Island Radiation: Recombination Rate Variation, Chromosome Structure, and Genome Assembly of the Hawaiian Cricket Laupala

Samples

We generated three F₂ interspecies hybrid families to estimate genetic maps. Multiple F₁ male and sibling females were intercrossed to generate F₂ mapping populations for the following species crosses: (1) a L. kohalensis female and L. paranigra male [“ParKoh”; 178 genotyped F₂ hybrid offspring; previously reported in Shaw et al. (2007)], (2) a L. kohlanesis female and a L. pruna male (“PruKoh”; 193 genotyped F₂ hybrid offspring), and (3) a L. paranigra female and a L. kona male (“KonPar”; 263 genotyped F₂ hybrid offspring). These four species are part of a recently radiated clade showing conspicuous mating song divergence (Mendelson and Shaw 2005). Approximate geographic distributions of the species, phylogenetic relationships, and parent collection localities are shown in Figure 1 and in Supplemental Material, Table S1. Crickets used in crosses were a combination of laboratory stock and outbred individuals [L. kohalensis (for ParKoh) and L. paranigra (for ParKoh and KonPar) were both reared in the laboratory for 3–15 generations; L. kohalensis (for PruKoh), L. pruna, and L. kona were wild caught]. All parental and hybrid generations were reared in a temperature-controlled room (20°) on Purina cricket chow and were provided with water ad libitum.

Genotyping

DNA was extracted from whole adults using the DNeasy Blood and Tissue Kits. Genotype-by-sequencing library preparation and sequencing were done at the Genomic Diversity Facility at Cornell University in 2014 following Elshire et al. (2011). The PstI restriction enzyme was used and DNA was sequenced (1 × 100 bp) on the Illumina HiSequation 2000 platform.

Reads were trimmed and demultiplexed using Flexbar version 2.5 (Dodt et al. 2012) and then mapped to the L. kohalensis de novo draft genome using Bowtie2 version 2.2.6 (Langmead and Salzberg 2012) with default parameters. We then called SNPs using two different pipelines: the Genome Analysis Toolkit (GATK) version 3.6.0 (DePristo et al. 2011; Van der Auwera et al. 2013) and FreeBayes version 0.9.13 (Garrison and Marth 2012). For GATK, we used individual BAM files to generate gVCF files using “HaplotypeCaller” followed by the joint genotyping step “GenotypeGVCF.” We then evaluated variation in SNP quality across all genotypes using custom R scripts (R Development Core Team 2016) to determine appropriate settings for hard filtering. The following metrics (based on the recommendations for hard filtering at https://gatkforums.broadinstitute.org/gatk/discussion/6925/understanding-and-adapting-the-generic-hard-filtering-recommendations) were used: quality-by-depth, Phred-scaled P-value using Fisher’s exact test to detect strand bias, root mean square of the mapping quality of the reads, u-based z-approximation from the Mann–Whitney rank sum test for mapping qualities, and u-based z-approximation from the Mann–Whitney rank sum test for the distance from the end of the read for reads with the alternate allele. For FreeBayes, we called variants from a merged BAM file using standard filters. After variant calling, we filtered the SNPs using “vcffilter,” a Perl-library part of the VCFtools package (Danecek et al. 2011), based on the following metrics: quality (>30), depth of coverage (>10), and strand bias for the alternative and reference alleles (SAP and SRP, both >0.0001). Finally, the variant files from the GATK pipeline and the FreeBayes pipeline were filtered to only contain biallelic SNPs with <10% missing genotypes using VCFtools version 0.1.15.

We retained two final variant sets: a high-confidence set, including only SNPs with identical genotype calls between the two variant discovery pipelines; and the full set of SNPs, which included all variants called using FreeBayes but limited to positions that were shared among the GATK and FreeBayes pipelines.

Linkage mapping

The genotype information from the parental lines was used to assign ancestry to the SNP loci. The parents of the crosses were heterogeneously heterozygous and only ancestry-informative loci were retained, i.e., all loci for which one or both parents were heterozygous were discarded. We were unable to obtain sequence data from the parents for PruKoh but we used sequence data from a single, nonparental L. pruna female and three available L. kohalensis females, which were all from the same populations as the parents. Ancestry was inferred if all three L. kohalensis individuals were homozygous for one allele and the L. pruna individual was homozygous for the alternative allele. All other loci were discarded. The loci were then further filtered based on genotype similarity and segregation distortion (see below for details).

The linkage maps deriving from the three species crosses were generated independently and by taking a three-step approach, employing both the regression mapping and the maximum likelihood (ML) mapping functions in JoinMap 4.0 (Van Ooijen 2006) as well as the three-point, error-corrected ML mapping function in MapMaker 3.0 (Lander et al. 1987; Lincoln et al. 1993). We used an intercross (F₂) scheme in both programs because we limited our SNP markers to those that are fixed between the parents (i.e., the only genotypes are A, B, and H).

In the first step, we estimated “initial” maps that are relatively low resolution (5 cM) but with high marker-order certainty. For initial maps, we first grouped (3.0 ≤ LOD score ≤ 5.0) and then ordered the high-confidence markers that showed no segregation distortion (markers with chi-square-associated P-values for deviation from Mendelian inheritance of <0.05 were discarded) and for which no marker had >95% similarity in genotypes across individuals compared to other markers (otherwise, one of each pair was excluded). When excluding similar loci, we favored those marker loci shared among the three mapping populations over markers unique to one or two crosses. We then checked for concordance among the three mapping algorithms. In most cases, the maps had highly concordant marker orders; any (rare) marker-order discrepancies between maps for the same cross were resolved by choosing the order that resulted in the highest map likelihood, shortest total map length, and the best fit of markers based on pairwise recombination frequencies with neighboring markers (see also JoinMap 4.0 manual).

These initial maps were then filled out using MapMaker with marker loci passing slightly more lenient criteria: markers were drawn from the full set of SNPs if they had false discovery rate (FDR)-corrected (Benjamini and Hochberg 1995) P-values for chi-square tests of deviation from Mendelian inheritance of ≤0.05 and if <99% of their genotypes were in common with other markers. First, more informative markers (no missing genotypes, >2.0-cM distance from other markers) were added, satisfying a log-likelihood threshold of 4.0 for the positioning of the marker (i.e., assigned marker position is 10,000 times more likely than any other position in the map). Remaining markers were added at the same threshold, followed by a second round for all markers at a log-likelihood threshold of 3.0. We then used the ripple algorithm on five-marker windows and explored alternative orders.

In the second step, “comprehensive” maps were obtained in MapMaker by sequentially adding markers from the full set of SNPs that met the more lenient criteria described above to the initial map. Markers were added if they satisfied a log-likelihood threshold of 2.0 for the marker positions, followed by a second round with a log-likelihood threshold of 1.0. We then used the ripple algorithm again on five-marker windows and explored alternative orders. Typically, MapMaker successfully juxtaposes SNP markers from the same scaffold. However, in marker-dense regions with low recombination rates, the likelihoods of alternative marker orders coalesce. In such regions, when multiple markers from the same genomic scaffold were interspersed by markers from a different scaffold, we repositioned the former markers by forcing them in the map together. If the log likelihood of the map decreased by >3.0 (factor 1000), only one of the markers from that scaffold was used in the map. The comprehensive maps provide a balance between marker density and confidence in marker ordering and spacing.

The third step was to create “dense” maps. We added all remaining markers that were not yet incorporated in step two, first at a log-likelihood threshold of 0.5, followed by another round at a log-likelihood threshold of 0.1. We then used the ripple command as described above. The dense maps are useful for anchoring of scaffolds and for obtaining the highest possible resolution of variation in recombination rates, but with the caveat that there is some uncertainty in marker order. Uncertainty is expected to be higher toward the centers of the linkage groups (LGs) where crossing-over events between adjacent markers become substantially less frequent (see Results).

Collinear linkage maps

We obtained 815,109,126; 522,378,849; and 311,558,401 reads after demultiplexing for ParKoh, KonPar, and PruKoh, respectively. Average sequencing depth ± SD across all individuals in the F₂ mapping population after filtering, (before and) after marker selection based on segregation distortion and ancestry information for linkage mapping was (62.4× ± 162.5 SD) 52.2× ± 31.4 SD, (44.3× ± 58.5 SD) 38.1× ± 23.8 SD, and (56.1× ± 105.7 SD) 41.8× ± 29.3 SD, respectively.

In the initial maps, 158 (ParKoh), 170 (KonPar), and 138 (PruKoh) markers were grouped into eight LGs at a LOD score of 5.0, corresponding to the seven autosomes and one X chromosome in Laupala. Markers group in eight LGs for a range of LOD scores: 5.0 < LOD score < 15.0 for ParKoh (no markers unlinked), 5.0 < LOD score < 42.0 for KonPar (no markers unlinked), and 4.0 < LOD score < 14.0 for PruKoh (between four and seven markers are unlinked depending on the LOD score threshold). The corresponding marker spacing was 5.14, 4.85, and 7.33 cM. The comprehensive maps contained 526, 650, and 325 markers with an average marker spacing of 1.91, 1.37, and 3.25 cM, respectively; and on the dense maps we placed 608, 823, and 383 markers with on average 1.69, 1.17, and 3.81 cM, respectively, between markers.

The recent divergence times and the limited levels of postzygotic isolation observed in this system led us to hypothesize that the linkage maps would show a high degree of collinearity. The visual comparison of marker positioning showed that the relative locations of shared scaffolds were similar across the linkage maps in both the initial and the comprehensive maps (Figure 2 and Figure S1). However, we also observe substantial variation in the total genetic length of homologous LGs, indicating recombination rate variation (Figure 2 and Figure S1). This variation may in part result from chromosomal rearrangements. However, we can only reliably detect rearrangements in our maps if they are alternatively fixed between L. pruna and L. kohalensis on the one side and L. paranigra and L. kona on the other side. In that specific scenario, the inverted marker order is visible when contrasting the PruKoh and KonPar maps, while the ParKoh map would show reduced recombination in that area (Figure 1C). Despite the apparent variation in recombination rates among homologous LGs, Spearman’s rank correlation of pairwise linkage-group comparisons was high (ρ varied between 0.91 and 1.00) and similar to values seen in comparisons of intraspecific linkage maps (e.g., Poursarebani et al. 2013); the quantitative measure of collinearity was largely consistent across LGs and across cross types (Table 2). Finally, merging the maps into a consensus pseudomolecule assembly allowed us to measure the error between individual maps and the merged assembly. Correlations between linkage maps and the pseudomolecule assembly were generally high (>0.95), indicating substantial synteny (Figure S2).

Open in a separate window

Figure 2

Initial linkage maps. Bars represent LGs for ParKoh, KonPar, and PruKoh. Lines within the bars indicate marker positions. The scale on the left measures marker spacing in cM. Blue lines connect markers on the same scaffold between the different maps. The map for ParKoh is shown twice to facilitate comparison across all three maps. See Figure S1 for comprehensive maps.

Limited heterogeneity in segregation distortion

We expected genetic incompatibilities to be more likely to occur in the ParKoh cross than in the KonPar and PruKoh cross, because L. kohalensis and L. paranigra are more distantly related than any of the other species pairs (Figure 1). We tested this hypothesis by examining the degree of segregation distortion in markers within 10-cM sliding windows across the linkage maps. Overall, segregation distortion was limited and average genotype frequencies were close to Mendelian expectations (Figure 3). However, LG3 showed a bias against L. kohalensis homozygotes in the ParKoh cross but not in any of the other crosses. Additionally, there was significant variation in the frequency of heterozygotes across the LGs [linear model: frequency(heterozygotes) ∼ LG × cross: R² = 0.21, F_20,1547 = 20.7, P < 0.0001]. The post hoc Tukey honest significant differences test revealed that LG7 had the lowest abundance of heterozygotes overall and within each of the intercrosses, and that levels of heterozygosity on LG7 were similar across the maps (Table S2). Together, these results show that for some LGs and in some crosses, certain genotype combinations were less common than expected, potentially as a result of genetic incompatibilities or meiotic drive.

Open in a separate window

Figure 3

Segregation distortion. For each of the seven autosomal LGs within the three comprehensive linkage maps (from top to bottom: ParKoh, KonPar, and PruKoh), a sliding window of the negative log-transformed P-values for the χ2-square test for deviation from a 1:2:1 segregation ratio is shown across markers with black lines in the top panels. In the panel below, the trace of the frequency of heterozygote genotypes (blue lines) and homozygote genotypes for both parental alleles (black and red lines, respectively) is shown. For each intercross, dashed gray lines indicate P = 0.01 (top panels) or expected allele frequencies based on 1:2:1 inheritance (bottom panels).

Variable recombination rates across the genome

We anchored a total of 1054 scaffolds covering 720 million bp, a little below half the current genome assembly. The pseudomolecule length ranges from 25,387,579 bp (LG6) to 137,279,277 bp (LG3), with N50 values between 569,059 bp (LG6) and 1,249,941 bp (LGX) (see Table S3 for scaffold number, N50, and assembly size per LG; and Figure S3 for coverage variation across the LGs). This gives us enough power to make inferences about broadscale recombination rate variation, but not about the existence of small-scale recombination hotspots. Average recombination rates varied from between 0.75 (KonPar) and 0.93 cM/Mb (ParKoh) on the X chromosome to between 3.12 (KonPar) and 4.24 cM/Mb (PruKoh) on LG6 (Table 3). We note that the recombination rate for LG6 might be artificially inflated because of lower assembly quality (expressed as N50) of this LG relative to the other LGs in all linkage maps and in the pseudochromosomes (Table S3). Both including and excluding the sex chromosome, there is a significant linear relationship between chromosome size and genetic length (linear mixed effect model with cross as random variable; with X: β = 0.62, F_1,23 = 14.95, P = 0.0008; without X: β = 0.69, F_1,20 = 29.7, P < 0.0001) and between chromosome size and broadscale recombination rate (with X: β = −34.1, F_1,23 = 29.1, P < 0.0001); without X: β = −24.0, F_1,23 = 63.7, P < 0.0001).

Most LGs showed wide regions of strongly reduced recombination rates in the center of the LGs (Figure 4). The general pattern of peripheral peaks in recombination rates juxtaposing large recombination “deserts” was observed in all three intercrosses, but some additional cross-specific peaks in recombination rates were observed on almost all LGs (Figure 4).

Open in a separate window

Figure 4

Recombination and Marey maps. Grayscale symbols and lines indicate the relationship between the physical distance (scaffold midposition) in Mb on the x-axis and the genetic distance in cM for each of the eight LGs on the left y-axis. ○’s represent the dense ParKoh linkage map, ▵’s and ⋄’s represent that of the KonPar and PruKoh cross, respectively. The corresponding lines (ParKoh: solid; KonPar: dashed; PruKoh: dotted) indicate the fitted smoothing spline (10 d.f.). The red lines (same stroke style) show the first order derivative of the fitted splines and represent the variation in recombination rate as a function of physical distance (cM/Mb, on the right y-axis). Gray bars indicate the approximate location of male song rhythm QTL peaks. The yellow ☆ in the LG1 panel highlights the QTL peak that colocalizes with a female preference QTL peak (Shaw and Lesnick 2009).

Collinearity of genetic maps

Based on the recent divergence (Mendelson and Shaw 2005) and strong premating isolation of Laupala species in the absence of conspicuous morphological and ecological differences (Otte 1994; Shaw 1996; Mendelson and Shaw 2005), we expected limited variation in chromosome structure and few signatures of genetic incompatibility between the species. In line with these expectations, we found that LGs are collinear across interspecies crosses (Figure 2). This was true both for comparisons of nonindependent species pairs (between ParKoh and the other two crosses), as well as for the independent contrast of the PruKoh map vs. the KonPar map. We saw some instances where markers occupied regions that may have been translocated or inverted. However, these instances were rare (Figure 2 and Figure S1) and recombination rates were similar among homologous LGs across the hybrid families (Figure 4). Moreover, the quantitative measures of collinearity, Spearman’s rank correlation among maps and between maps and the pseudomolecule assembly, as well as limited segregation distortion (but see discussion of one exception below) both supported the collinearity hypothesis; if there were larger or many small inversions between the species pairs, the order of markers in the different maps would be more variable, in which case significantly lower values for Spearman’s rank correlation would be expected than the uniformly high correlations (most values are between 0.95 and 0.99) we observe here.

Variation in the organization and structure of chromosomes can contribute to postzygotic reproductive isolation after speciation as well as to the speciation process directly (Noor et al. 2001; Rieseberg 2001). We conclude that, at least for the Laupala group that radiated on the Big Island of Hawaii in the last 500,000 years, structural rearrangements have not played a major role in the evolution of reproductive isolation. This is because, similar to two hybridizing Heliconius species (Davey et al. 2017), we observed that chromosome-wide recombination rates are relatively conserved and large chromosomal rearrangements are absent. We hypothesize that, for Laupala on the Big Island, premating isolation combined with (partial) geographic separation (i.e., low migration rates) provides a sufficiently strong barrier to gene flow between sister species. Indeed, as has been shown in recent models of the role of inversions in speciation, genomic rearrangements can only invade and spread in diverging populations if levels of gene flow and the contribution of structural variation to isolation (by linking adaptive alleles or incompatibilities) is high relative to the strength of assortative mating (Feder et al. 2014; Dagilis and Kirkpatrick 2016).

We acknowledge that the power to detect rearrangements and recombination rate variation between the genomes of the species studied here is limited by the resolution of the maps, the sample size, and the error rate of the sequencing platform. The average spacing of markers is between 1.37 and 3.25 cM; the resolution at which we can detect rearrangements is on the order of 10⁵ and 10⁶ bp. Therefore, it is difficult to attribute subtle variation in marker order and genetic distance (<1–2 cM) between the maps to genomic rearrangements, as opposed to mapping errors and sampling variance, because variation in recombination frequencies between closely related species might be detected only at finer scales (Stevison et al. 2017). However, the overall collinearity of the maps supports our conclusion that chromosomal architectures are broadly conserved between the species studied here.

Genetic incompatibilities

We expected genetic incompatibilities to be more likely between genomes of more distantly related species. Accordingly, we discovered a single region covering approximately half of LG3 with high segregation distortion in ParKoh; we found no such deviations in the other two crosses (Figure 3). Inspection of genotype frequencies indicated that fewer individuals than expected were homozygous for L. kohalensis alleles at loci in this region. In a controlled cross, segregation distortion can be caused by prezygotic effects such as meiotic drive of selfish genetic elements and distorter genes [e.g., like sd in Drosophila melanogaster (Larracuente and Presgraves 2012)], and by postzygotic genetic incompatibilities (Dobzhansky 1937; Muller 1942; Burt and Trivers 2006; Presgraves 2010). Genotypic errors may produce superficially similar patterns but are unlikely to distort segregation over large genomic regions and with consistent bias toward the same genotypes. Although meiotic drive is a possible alternative to genetic incompatibilities, we do not see the same effect in the other cross involving L. paranigra, where selfish genetic elements or segregation distorters ought to have a similar effect. Overall, the large region on LG3 reveals a potential local postmating barrier to gene flow that could contribute to strengthening existing prezygotic barriers in secondary contact zones or following episodes of migration.

Recombination landscape

Chromosomal rearrangements influence genomic divergence by locally altering recombination rates within and among species. Felsenstein (1974, 1981) illuminated the role of intraspecific recombination in purging deleterious alleles and the role of interspecific recombination in decoupling coadapted alleles. In recent years, the role of recombination and its interaction with divergent selection and adaptation on genomic scales have received considerable attention (e.g., Yeaman and Whitlock 2011; Feder et al. 2012; Samuk et al. 2017) and technological advances are shifting focus toward characterizing recombination rates across genomes (Butlin 2005; Slatkin 2008; Noor and Bennett 2009; Barb et al. 2014; Burri et al. 2015). Comparisons of genomic recombination landscapes illuminate the role of recombination in speciation, but focus has been on systems where environmental and ecological pressures are driving divergence between locally adapted populations.

Here, we show for species separated by conspicuous sexual barriers that there is limited variation in recombination rates across the maps of three interspecific crosses (Figure 4), but strong heterogeneity in recombination rates across the genome. Genome-wide average interspecific recombination rate varied between 1.3 and 1.5 cM/Mb (Table 3), similar to intraspecific rates observed in dipterans and substantially lower than social hymenopterans and lepidopterans (Wilfert et al. 2007). We note that our estimates are derived from interspecific maps, which may lead to somewhat lower estimates compared to intraspecific maps (e.g., Beukeboom et al. 2010), because genetic incompatibilities and rearrangements may reduce rates of crossing over; however, differences between intra- and interspecific recombination might be negligible if rearrangements are rare (e.g., Davey et al. 2017). Moreover, we anchored ∼50% of the nucleotides in the draft assembly to LGs, and there remain many scaffolds not mapped to a genomic position. These “missing” scaffolds are expected to add to the physical length of the chromosomes more so than to the genetic length of the chromosomes, thus lowering the overall recombination rate and increasing heterogeneity in recombination rates among the crosses. However, our study emphasized relative patterns of recombination, which should not be affected by our sampling. While we can only approximate intraspecific recombination rates at this point, we note that recent divergence of the study species and collinearity of the linkage maps support conservation of recombination landscapes across intraspecific and interspecific comparisons.

Interestingly, for all three species pairs we document high variability in interspecific recombination across genomic regions. We found a stereotypical U-shaped pattern: large regions of low recombination in the center, with recombination rates well below 1 cM/Mb and occasionally approaching 0, flanked by steep inclines reaching rates up to 6 cM/Mb (Figure 4). This pattern is consistent with earlier findings in plants (Anderson et al. 2003), invertebrates (Rockman and Kruglyak 2009; Niehuis et al. 2010), and vertebrates (Backström et al. 2010; Roesti et al. 2013; Singhal et al. 2015) including humans (Kong et al. 2002); but differs from observations in, for example, Drosophila (Kulathinal et al. 2008), which show heterogeneity in recombination rates, but not necessarily much higher rates on the periphery of the chromosomes. Commonly invoked drivers of local recombination suppression, such as selection against recombination due to negative epistasis or the maintenance of linkage disequilibrium between mutually beneficial alleles (Smukowski and Noor 2011; Stevison et al. 2011; Smukowski Heil et al. 2015; Ortiz-Barrientos et al. 2016), are not likely to leave chromosome-wide signatures. Rather, the observed U-shaped pattern is more likely attributable to structural properties of chromosomes, such as the location of the centromere and heterochromatin-rich regions (Copenhaver et al. 1999; Haupt et al. 2001). Roesti et al. (2013) observe similar recombination landscapes in stickleback fish and suggest it might be due to peripheral clustering during meiosis prophase I to facilitate homolog pairing (Harper et al. 2004; Brown et al. 2005). Regardless of the mechanism, the observed genomic architecture will drive substantial heterogeneity in the propensity of favorable and/or maladaptive alleles to come together, break apart, and introgress in heterospecific backgrounds.

We thank Stephen Chenoweth and two anonymous reviewers for helpful comments that strongly improved the quality of this manuscript. We further thank the Shaw laboratory, in particular Mingzi Xu, as well as Michael Sheehan and other members from Cornell’s Department of Neurobiology and Behavior for input that contributed to the interpretation of the findings. This work was supported by the National Science Foundation (DEB 1241060, IOS 1257682, and IOS 0843528).