Ultraconserved elements (UCEs) illuminate the population genomics of a recent, high-latitude avian speciation event

Introduction

Among non-model organisms, population genetic studies have used a diverse set of markers, tending to concentrate on those with sufficiently high substitution rates to provide useful data at shallow levels of evolutionary divergence, for example, from the populations-to-species levels (Avise, 1994; Hillis, Moritz & Mable, 1996; Pearse & Crandall, 2004). This approach usually provides answers to the specific questions asked by researchers, but the historic focus on markers with high substitution rates has produced studies that include relatively few loci and often have little to no overlap with loci used for other taxa. This lack of consistency in the loci used across different studies compromises our ability to make direct comparisons of population genetic parameters among taxa (e.g., in divergence statistics and in estimates of gene flow and effective population sizes). Improvements in sequencing platforms and genomic data collection approaches are changing this general pattern by enabling us to efficiently collect much larger samples of the genome, up to and including whole-genome sequences (Ellegren, 2014). However, the sheer quantity of data obtained from whole-genome sequencing can require excessively long computation times, and may be overkill for many questions. The parallel difficulties of collecting a moderate sample of the genome from identical loci across diverse species argue for a sequence data collection approach that (a) subsamples the genome to (b) obtain orthologous markers across a broad taxonomic scope. This type of approach would provide a tractable number of loci for analyses while improving among-study comparisons and larger-scale comparative meta-analyses. Ultraconserved elements (UCEs) are one class of genome-wide marker that might provide a solution to these problems.

Ultraconserved elements are conserved sequences shared among divergent animal genomes (Bejerano et al., 2004; Siepel et al., 2005; Stephen et al., 2008; Janes et al., 2011), and many UCE loci are likely to be involved in controlling gene expression (Marcovitz, Jia & Bejerano, 2016). UCEs in vertebrates show little overlap with most types of paralogous genes, and, as a marker class, UCE loci are broadly distributed across the genome and are typically transposon-free (Derti et al., 2006; Simons et al., 2006; McCormack et al., 2011; Harvey et al., 2016). We focus on the set of UCEs previously defined for tetrapods and now in widespread use (McCormack et al., 2011; Faircloth et al., 2012). Outside of their functional relevance, UCE loci have demonstrated utility for recovering deeper-level phylogenetic relationships (McCormack et al., 2013; Faircloth et al., 2015; Gilbert et al., 2015) and shallower-level genus and population relationships (Smith et al., 2014; Harvey & Brumfield, 2015; Leaché et al., 2015; Harvey et al., 2016; Manthey et al., 2016; Oswald et al., 2016; Mason et al., 2018). Although UCEs are highly conserved at their core, which enables universal capture of loci across diverse groups of organisms (Faircloth et al., 2012, 2015; Faircloth, 2013; Starret et al., 2017), lower levels of purifying selection away from the core allow substitutions to accumulate in the flanking regions. Using human genome data, Faircloth et al. (2012) demonstrated that the increased variation in UCE flanking sequence might be adequate to make these loci useful for questions at shallow levels of divergence, and this hypothesis has been supported by subsequent empirical studies (Smith et al., 2014; Harvey & Brumfield, 2015; Harvey et al., 2016; Oswald et al., 2016; Mason et al., 2018). However, the utility of UCE loci for studying population genetics, population divergence, and/or incipient speciation is only beginning to be tested, and both the value and challenges of using UCEs at these shallow levels remain underexplored.

Here, we examine the utility of UCEs for studying the population genomics of divergence between two bird species, McKay’s bunting (Plectrophenax hyperboreus) and snow bunting (P. nivalis). McKay’s buntings breed on remote islands in the Bering Sea (where our samples are from) and are the highest-latitude endemic songbirds; their range is restricted to the North Pacific region. Snow buntings breed throughout the rest of the high-latitude Holarctic (our samples are from the southern edge of the Being sea on the Alaska Peninsula and an Aleutian island; Table S1). McKay’s bunting is thought to have arisen ∼18–74 Kya during the last glacial maximum through divergence from snow buntings, and previous work suggests gene flow between the two may be ongoing (Maley & Winker, 2010). These species are interesting to study using UCEs because prior work (Maley & Winker, 2010) enables us to compare population genetic statistics derived from UCEs versus traditional population genetic markers (mtDNA sequence and amplified fragment length polymorphisms, AFLPs). These species also allow us to test the utility of UCEs for studying very shallow divergences between sister lineages where gene flow may be ongoing.

Methods

Data analysis

We calculated coverage depths, SNP positions within loci, and SNP-specific and locus-specific F_ST values on the complete VCF file using VCFtools (ver. 0.1.12b; Danecek et al., 2011). After thinning the VCF file to 1 SNP/locus (which is required in demographic analyses when unlinked variation is important) and converting the VCF file to STRUCTURE format using PGDSpider (ver. 2.1.0.3; Lischer & Excoffier, 2012), we performed tests of Hardy–Weinberg equilibrium and computed observed and expected heterozygosities, homogeneity of variance, population structure (population F_ST, including a 10,000-replicate G-test; see Goudet et al. (1996)), and the probabilities of each individual’s assignment to a particular population using discriminant analysis of principal components (DAPC) in adegenet (ver. 2.0.1; Jombart & Ahmed, 2011). To calculate nucleotide diversity, average distance between taxa (d_xy), and net average distance (d_A), we created a concatenated FASTA file of all individual sequences using catfasta2phyml by Johan Nylander (https://github.com/nylander/catfasta2phyml), and we analyzed this file in MEGA (ver. 6; Tamura et al., 2013) using the maximum composite likelihood method.

We estimated recombination using the four-gametes test as implemented in IMgc (Woerner, Cox & Hammer, 2007), which also produces sequence datasets from which the effects of recombination have been removed. Resulting sequences were used for IMa2p (ver. 1.0; Sethuraman & Hey, 2015) analyses to attempt to estimate demographic parameters, but those analyses did not converge under a variety of full- and sub-sampling schemes and are not reported. We nevertheless include the results of IMgc because accounting for recombination is a critical part of workflows using full sequences (i.e., not just SNPs), and these results provide needed insight into the levels of recombination found in UCE loci for studies of this type.

We used Diffusion Approximations for Demographic Inference (δaδi; ver. 1.7.0; Gutenkunst et al., 2009) to infer demographic parameters from the data under a variety of divergence scenarios (models) after excluding Z-linked loci (for δaδi analyses only). Z-linked loci in birds are on the sex chromosome, have a different inheritance scalar from autosomal loci, and sample population sex ratios affect allele frequency estimates (Jorde et al., 2000; Garrigan et al., 2007). We identified Z-linked loci in our data using BLASTn (ver. 2.3.1; Zhang et al., 2000), by aligning the reference set of UCE loci against the zebra finch (Taeniopygia guttata) genome (NCBI Annotation Release 103). We excluded UCE loci that strongly matched (E-values ∼0.0) the zebra finch Z chromosome. After removing Z-linked loci from our complete VCF file, we converted this reduced dataset to biallelic format (which dropped one locus with >2 alleles at a SNP site) and thinned the data to one SNP per locus using VCFtools. Then we converted the resulting VCF file to the joint site frequency spectrum (SFS) format required by δaδi using a PERL script by Kun Wang (https://groups.google.com/forum/#!msg/dadi-user/p1WvTKRI9_0/1yQtcKqamPcJ). Because we lacked an outgroup, we used a folded SFS in our analyses (Gutenkunst et al., 2009), which lacks polarization of SNPs (Fig. S1).

Because we had prior evidence that these species represent two genetic populations (based on taxonomy and results from Maley & Winker, 2010), we used δaδi to infer what general two-population divergence model best fit the data. We then used that model to estimate demographic parameters (i.e., effective population sizes, split time, and migration). We ran six different models spanning the standard possible demographic histories of two populations, five basic and one derivative: (1) neutral (no divergence, or still strongly mixing), (2) split with migration, (3) split with no migration, (4) isolation with bidirectional migration and population growth, (5) isolation with population growth and no migration, and (6) a custom split-bidirectional-migration model (a simple derivative of split-migration; Fig. 1). The neutral, split-with-migration, and isolation-with-migration-and-population-growth models are provided in the δaδi file Demographics2D.py as snm, split_mig, and IM, respectively. The no-migration models (3 and 5 above) use the split_mig and IM models, respectively, with migration parameters set to zero. The split-bidirectional-migration model (figshare https://doi.org/10.6084/m9.figshare.6453125.v1) adds bidirectional migration to the split-migration model to examine potential asymmetry in gene flow.

We performed a series of optimization runs (8–47 each) of each basic model, adjusting parameters (grid points, upper and lower bounds) to identify high log composite likelihoods. We then ran each model repeatedly, varying parameters within bounds that yielded the highest likelihood during optimization, until three runs yielded the highest observed likelihood value. Our reasoning was that this level of repeatability indicated a best-fit neighborhood for each model. We report this highest observed likelihood, except for poorer models, which yielded variable likelihoods, in which case we averaged and report the highest five values. After identifying the best-fit model based on likelihood values over successive runs, we ran the best-fit model 10 times each with jackknifed datasets to estimate the 95% confidence interval (CI) for each parameter.

We estimated the average per-site substitution rate by BLASTing the FASTA file containing all confidently scored loci (those meeting our quality filters as described in the Bioinformatics section, above) for all individuals (3,431 loci) against the budgerigar (Melopsittacus undulatus) genome (NCBI release 102) and the rifleman (Acanthisitta chloris) genome (NCBI release 100), using time to most recent common ancestor (TMRCA) date estimates of 60.5 Ma (budgerigar) and 53 Ma (rifleman) (Claramunt & Cracraft, 2015). These taxa were chosen as the nearest relatives with complete genomes and fossil-dated nodes available. We imported BLAST results (as a hit table, csv) into a spreadsheet, removed duplicate lower-affinity hits, then summed total length of base pairs, total substitutions, and calculated substitutions per site. This value (substitutions per site) was annualized by multiplying it by 2 TMRCA (e.g., 121 Ma for the total time along the branches of divergence between buntings and budgerigars). To account for the uncertainty associated with using divergence time estimates from distant relatives, we averaged the resulting substitutions per site per year rates (6.83 × 10⁻¹⁰ and 6.67 × 10⁻¹⁰, respectively), and we used the average rate (6.75 × 10⁻¹⁰) to convert parameter estimates obtained from δaδi analyses into biologically relevant estimates of effective population size(s) and split time(s). We converted substitution rates to substitutions/site/generation using a generation time of 2.7 years for snow buntings. We estimated generation time using survival and breeding data from Smith (1994) and the method of Sæther et al. (2005), in which generation time (G) is calculated as G = α + (s/(1 − s)), where α is age of first breeding (1 year) and s is annual adult survival. Finally, because our δaδi analyses used filtered SNPs, our demographic estimates used an adjusted surveyed sequence length of: (total sites surveyed, including invariant sites) × (proportion of SNPs used after thinning to 1 SNP/locus) = 1,103,715 bp.

Results

Assembly produced 632,401 contigs (min = 224 bp, max = 17,453 bp) with a mean length of 396.6 bp (±0.27 bp 95% CI) for a total of 250,802,355 bp. Fully 9,194 contigs were over one Kb in length. After identifying UCE loci and removing potential paralogs, we recovered 4,018 UCE loci. After filtering UCE loci for quality, calling SNPs, phasing (reconstructing haplotypes), and applying additional quality filters, we identified 2,635 loci that contained data for all individuals and were variable. This complete matrix of variable loci included a total of 9,449 SNPs (averaging 3.6 sites per locus). Per-site sequencing depth for these SNPs averaged 26.3 reads (±16.9 SD). An additional 587 loci exhibited variation but the data were not of sufficient quality (i.e., GQ < 10) among all individuals to confidently call both alleles. There were 796 high-quality invariant loci (loci with invariant data, rather than an absence of data), providing a full dataset of 3,431 loci with mean length of 1153.6 bp (±4.95 bp 95% CI). The shortest locus was 228 bp, the longest 2,543 bp, and 2,482 loci were longer than one Kb (Fig. S2). The total length of these loci was 3,957,876 bp. The distribution of SNP variation among loci confidently called for all individuals is given in Fig. 2. Nucleotide diversity (π) was 0.000519 overall, 0.000523 for snow buntings, and 0.000493 for McKay’s buntings.

No alleles showed fixed differences (F_ST = 1.0) between the two populations, and few alleles showed strong segregation. No variable sites had an F_ST value above 0.9, and there were only three each at 0.86 and 0.72 (Fig. S3; two of these sites were on the same locus). One of the five loci with the highest F_ST values was Z-linked; all of the others were on different chromosomes (figshare https://doi.org/10.6084/m9.figshare.6453125.v1). There were 128 Z-linked loci among the 2,635 variable loci. As noted, only one showed high F_ST between the two species. The two populations had an overall F_ST = 0.034, which was significant (P = 0.03). The average distance between taxa (d_xy) was 5.3 × 10⁻⁴, and the net average distance (d_A) was 2.0 × 10⁻⁵. DAPC in adegenet assigned all individuals to their correct taxon of origin (retaining the first four PCs), with 100% probabilities for each, indicating a high level of genomic diagnosability (Fig. S4).

Fully 2,510 loci were in Hardy–Weinberg equilibrium; 124 were not (one was triallelic). McKay’s buntings had fewer unique alleles (4,238) than snow buntings (4,389), concordant with the smaller population size of McKay’s buntings. Bartlett’s test rejected homogeneity of variance between observed heterozygosity (H_o = 0.18, 0.19) and expected heterozygosity (H_e = 0.20, 0.22), but H_o did not differ from H_e (t = −3.1653, df = 2,633, P = 1.0).

The four-gametes test suggested that recombination occurred in hundreds of loci. For 405 loci, locus lengths were shortened by IMgc to meet the four-gametes test, and for 252 loci one or more individuals were removed to meet the same criteria (a few of these loci had both done; IMgc automatically performs one or the other or both operations to obtain non-recombinant sequence data). There were thus 15.4–24.9% of variable loci exhibiting patterns indicative of recombination. As noted in the Methods, these sequence data, together with all other unchanged sequences, were not used further; we used only SNP data for further analyses.

In testing our six, two-population models with δaδi, the highest maximum log composite likelihood values were obtained for the split-with-migration model (−112.76), which made it the best-fitting model for these data (model 2 in Fig. 1). We obtained successively lower likelihood values for the neutral (−588.45), isolation with bidirectional migration and population growth (−803.30), and isolation with population growth and no migration (−2026.93) models. The final model tested, split-bidirectional-migration, had an intermediate likelihood of −286.49. The split-with-no-migration model was unstable under all conditions tried, and we could not get it to run to convergence. We provide jackknifed estimates and CIs for the best-fitting, split-with-migration model in Table 1.

Discussion

Our data provided sufficient variation to answer fundamental questions about these two recently diverged taxa, despite a lack of fixed genetic differences and evidence for moderate levels of gene flow. Thus, our study adds to evidence showing the utility of UCEs for illuminating key evolutionary attributes among populations with shallow levels of divergence (Table 2). These data also provide a direct comparison to markers previously used to investigate recent divergence (i.e., mtDNA, AFLPs) for these same taxa (see below). As UCEs are used more frequently for population genomics, in addition to systematics, new actions become desirable (Table 2). Some of the key approaches are: sequencing at increased depth, genotyping individuals (determining both alleles of a locus), implementing GQ filters, accounting for recombination, improving mutation rate estimates, and implementing population genomics analytical pipelines rather than those oriented more typically toward systematics. Questions often differ at population levels, but researchers are successfully applying a variety of approaches that demonstrate the utility of UCEs in population genomics (Table 2).

In considering UCEs as a class of markers that subsamples the genome, it is useful to note that our estimated substitution rates (mean of 6.75 × 10⁻¹⁰ substitutions per site per year) are roughly an order of magnitude slower than the mutation rate estimated across the entire genome of three generations of Ficedula flycatchers (Smeds, Qvarnström & Ellegren, 2016). This is perhaps not surprising given the conserved nature of these loci. Using a very different method to estimate substitution rates (scaling UCE results to an mtDNA molecular clock), Harvey et al. (2017) estimated rates of 1.74 × 10⁻¹² to 2.32 × 10⁻¹¹ substitutions per site per year for UCE loci, one to two orders of magnitude slower than comparable RAD-seq data from the same animals and also slower than the rates we estimated here for buntings. In addition to differences in methodology, Harvey et al. (2017) had shorter loci on average (mean locus length 604 bp) than loci in our study. Nevertheless, more study of substitution rates in loci with UCEs is warranted because these estimates are important when converting modeled demographic parameters into biological units. The effect of our estimated substitution rates on our demographic estimates (if, e.g., our substitution rate estimates are wrong) is that for some, they are positively correlated; lower substitution rates would drive effective population sizes and split times lower (N_e: nu1, nu2, and Nref and T in Table 1). Migration rate (m) estimates in Table 1 are unaffected by substitution rates.

The nucleotide diversity levels that we observed are approximately an order of magnitude lower than typical levels across the avian genome (0.0011–0.005; Ellegren, 2013). This is likely the result of purifying selection acting on UCE loci, effecting an apparent lower substitution rate. Our values are more similar to the values for Z-linked loci in other bird species (Balakrishnan & Edwards, 2009; Huynh, Maney & Thomas, 2010; Lavretsky et al., 2016).

When applied to this relatively recently derived pair of taxa, UCE results raise the question of whether McKay’s bunting is a full biological species. Although McKay’s bunting is taxonomically recognized as a species, this dataset shows substantial levels of gene flow (see Wright, 1943; Cabe & Alstad, 1994; Winker, 2010), and the lack of fixed alleles is surprising given that we sampled thousands of loci and four individuals from each of two putative species; we would expect several fixed differences to occur by chance through neutral processes. There are some noted plumage differences between the two taxa (Maley & Winker, 2007). But while our results enabled 100% diagnosability (which might decline with broader sampling; it was 96.5% using AFLPs in Maley & Winker, 2010), they also suggest widespread genomic similarity between McKay’s and snow buntings (e.g., relatively low F_ST). Given phenotypic differences between the taxa, it seems likely that there are fixed allelic differences in portions of the genome not included in our data that could be detected by more extensive surveys of each species’ genome. The status of the taxa as biological species, however, is more likely to hinge on gene flow (i.e., the geographic partitioning of traits that may be responding to adaptation is not equivalent to speciation).

There are reports of male McKay’s buntings present outside their breeding range and possible hybridization between McKay’s and snow buntings (Sealy, 1967, 1969). Snow buntings are also common on the breeding range of McKay’s buntings at St. Matthew Island prior to and during early portions of the breeding season, although most individuals leave before fledging (Winker et al., 2002). Just one pair of snow buntings has been recorded on the island during fledging (Winker et al., 2002). Observations thus suggest the possibility of hybridization; our data provide a confirmation and a quantification of it. The levels of gene flow that we found, 2.8–5.2 individuals per generation (Table 1), seem rather high for two putative biological species (Rice & Hostert, 1993; Hostert, 1997; Winker, 2010). Further study will be needed to determine species limits between these taxa, including larger sample sizes, broader genomic coverage, and proper caution for interpreting genomic results in terms of species delimitation (Robinson et al., 2014; Sukumaran & Knowles, 2017).

In comparing UCE-based estimates of demographic parameters with those based on mtDNA sequence (Maley & Winker, 2010), we find little overlap (Table 3), and our UCE-based split-time estimate is an order of magnitude earlier. Although effective population size estimates for McKay’s buntings are close (though non-overlapping), those for snow buntings are one-to-two orders of magnitude smaller, a difference that is only partially explained by differences in effective population size for autosomal and mtDNA estimates. These differences may also be driven by the different selection regimes operating on the two marker classes. For example, purifying selection on UCEs will result in background selection on linked variation in flanking regions, reducing (through hitchhiking) the effective population size (Charlesworth & Charlesworth, 2016). Previously, mtDNA results suggested that gene flow was highly asymmetric (Maley & Winker, 2010), concordant with what was likely a post-glacial introgression of McKay’s buntings into snow buntings during a snow bunting range expansion into Beringia. Our UCE-based estimates have much narrower confidence limits (and without the strong asymmetry found in mtDNA; Table 3), but they do suggest moderate levels of gene flow between the two species.

Conclusions

Although more work is needed to understand demographic estimates made using UCEs relative to those obtained using other markers, UCEs provide rich, high-quality data for population genomic studies (Table 4). They are thus an important new class of genomic marker that should provide broad comparative value among diverse population genomics studies, with ever-increasing value as additional studies using UCEs (or whole genomes from which UCEs can be obtained) are conducted.

Supplemental Information