Mitogenomes reveal the timing and distribution of divergence events among trans-Beringian birds

Taxonomic selection

We selected thirty-nine lineage pairs based on specimen availability and distribution across Beringia. The lineages selected were from six avian orders and 14 families and broadly representative of the Beringian avifauna. For example, our strongest representation is in the three most diverse Beringian orders, with substantial percentages of each of these groups’ available lineages included: Anseriformes 42% (∼13/31), Charadriiformes 17% (∼12/69), and Passeriformes 19% (∼11/59; Table 1 and Winker et al., 2023). For species, we used the framework developed for Beringian birds given in Winker et al. (2023), in which clearly closely related congeners on each side of the Bering and Chukchi seas were considered sisters unless species-level phylogenies (which do not exist for all these genera) showed otherwise. One lineage (Phylloscopus spp.) has both subspecific and species-level divergence across Beringia. Lineages show variable biogeographic patterns, including taxa of predominantly Old World (OW; n = 202 individuals), predominantly New World (NW; n = 230 individuals), and Holarctic distributions, with a limited number of intra-Beringian comparisons, primarily involving endemic taxa. Individuals were categorized as OW or NW based on phenotype except for population-level contrasts, which were based on location. Migrant OW birds often occur in the western Aleutian Islands (NW), and some of these were used to increase sample sizes of OW specimens available (Table S1). Lineages varied in divergence level as reflected in current taxonomy, with 12 species-level, 11 subspecies-level, and 16 population-level divergences (Chesser et al., 2023; Table 1, Table S1).

Mitogenome dataset generation

We chose to use mtDNA sequence data for this study to provide a single-locus comparative framework for this system. This framework enables direct contrasts with historic data, the use of good published divergence calibrations, and a computationally tractable dataset for our multi-lineage comparative focus. We have studied some of these same taxa using ultraconserved elements (UCEs) and have found the UCE approach is unlikely (at this time) to provide robust answers to our focal questions (e.g., Winker, Glenn & Faircloth, 2018; McLaughlin et al., 2020; Spaulding et al., 2022; Oliver Brown, Glenn & Winker, 2025). Here we focus only on the mitogenome. Most taxonomic pairs had a minimum of five individuals of each population (or subspecies or species) per pair. We extracted DNA from high-quality frozen muscle tissues from the collections of the University of Alaska Museum and University of Washington Burke Museum (Table S1) using a DNeasy Tissue Kit, following the manufacturer’s protocol (Qiagen, Valencia, California, USA). Associated voucher specimens are archived at these institutions. We then prepared dual-indexed DNA libraries for each sample using the methods in Glenn et al. (2017), and calibrated them using a Qubit fluorimeter (Invitrogen Inc., Carlsbad, CA, USA). For each lineage in which a pairwise comparison was being made, libraries for all associated specimens were combined into an equimolar 500 ng pool, which was then enriched for ultraconserved element (UCE) loci following Winker, Glenn & Faircloth (2018). Fragment size distributions of each lineage-specific pool were quantified using a Bioanalyzer (Agilent Inc., Santa Clara, CA, USA), and pools were amplified and quantified with qPCR using a commercial kit (Kapa Biosystems, Wilmington, MA, USA). Pools were then combined at equimolar ratios and sequenced on one lane of a paired-end 150 bp (PE150) Illumina HiSeq 2500 (Illumina Inc., San Diego, CA, USA; UCLA Neuroscience Genomics Core). This procedure produces large amounts of mitogenomic data as a by-product, and that is our focus here.

We used Illumiprocessor (v.2.0.6; Faircloth, 2013) with all raw fastq files to trim low-quality (phred score <30; i.e., <99.9% accuracy) bases and clip adapter contamination. For each lineage, we then obtained a full annotated mitochondrial reference sequence from GenBank, from a conspecific if available. If not, we used a reference from a closely related species, typically a congener (Table S2). Phyluce (v.1.5.0; Faircloth, 2016) and several of its dependencies were used to call SNPs from reads in individual fasta files against these reference mitogenomic sequences, lineage by lineage. We next aligned raw reads for each individual to the reference and converted them to .bam files using bwa-mem (v.0.7.7; Li & Durbin, 2009), SAMtools (v.0.1.19; Li et al., 2009), and PICARD (v.1.106; http://broadinstitute.github.io/picard). We then cleaned all sequences to remove unmapped reads and beyond end-of-reference alignments, tagged them with read groups corresponding to their specimen ID number, and marked PCR and sequencing duplicates, also using PICARD. We then obtained coverage summary statistics by taxon (Tables S3, S4) from the sorted .bam files using SAMtools.

Using the Genome Analysis Toolkit (GATK, v.3.3.0; McKenna et al., 2010), we indexed each .bam file and called high-quality SNPs against the reference following an established population genomic pipeline (Harvey et al., 2016). We then generated individual mtDNA fasta sequences for each study individual, also using GATK.

Gene isolation and divergence estimates

We first aligned full mitogenomes in MEGA X (Stecher, Tamura & Kumar, 2020) using Muscle (v.3.8.1551; Edgar, 2004) on default settings. Due to downstream computational limitations (see ‘MTML-msbayes analyses’), we then created a data subset by trimming the alignments according to the coordinates of three mitochondrial genes—NADH dehydrogenase 2 (ND2), cytochrome b (cytb), and cytochrome c oxidase 1 (CO1) in the annotated reference sequences, producing alignments for each of these genes.

We then partitioned each aligned dataset into groups of OW and NW individuals based on taxonomic identity or specimen location, as applicable. To visualize specimen relatedness, we created haplotype networks and neighbor-joining (NJ) trees for each taxon. Haplotype networks were generated for the ND2 dataset with the median-joining algorithm in the program Network v.4.6.1.1 (Bandelt, Forster & Röhl, 1999). NJ trees for the ND2 dataset were generated in MEGA X. We then labeled and color-coded these trees in FigTree (v1.4.4; Rambaut, 2018).

We calculated between-group genetic distances using JC-distances with a gamma distribution in MEGA X, to match the model used in Campbell, Braile & Winker (2016) (K Campbell, pers. comm., 2022). We used Arlequin 3.5 (Excoffier & Lischer, 2010) to calculate net nucleotide divergence (D_A; Nei, 1983). Both metrics are common in the historic literature and thus enable direct comparisons with prior work.

For a subset of seven taxon pairs for which historic Sanger-sequencing data existed (Humphries & Winker, 2011), we also used MEGA X to calculate nucleotide diversity (π, Nei & Li, 1979) values in both the historic Sanger data and the current Illumina data. We used these values to inform comparisons between these datasets.

Finally, to estimate the proportion of variance in divergence attributable to taxonomic depth, we performed a linear regression of taxonomic level (population, subspecies, and species) against full mitogenomic JC-distance, as calculated above.

Dating divergence events

We estimated times of divergence (t) for pairwise comparisons for the JC-distance values of the whole mitogenome and ND2 datasets. For whole-mitogenome divergence date estimates, we used calibrations from Quillfeldt (2017) and Arcones, Ponti & Vieites (2021). For our ND2 divergence rate estimates, we used the calibration for shallow avian divergence events from Lerner et al. (2011). As Lerner et al. (2011) used only passerines to determine their rate, we adjusted this base value using ordinal rates from Arcones, Ponti & Vieites (2021) to the avian orders present in this study using the equation (Lerner_scaled = (Arcones_{Desired_lineage}/Arcones_{Passeriformes}) * Lerner_unscaled). This enabled us to make direct comparisons among taxa. Other mtDNA calibration approaches might be used (e.g., Nabholz, Lanfear & Fuchs, 2016), and can be applied directly to our pairwise divergence calculations, although as our results suggest, we consider that new, shallow calibrations are needed for data of the sort we generated. We made pairwise estimates of F_ST using DnaSP v.5 (Librado & Rozas, 2009).

MTML-msBayes analyses

We assessed the simultaneous versus non-simultaneous divergence (estimating the variation in divergence times) of all taxon pairs using a hierarchical approximate Bayesian computational (hABC) approach implemented in MTML-msBayes (Huang et al., 2011). hABC is a likelihood-free Bayesian inference framework that uses stochastic simulations and summary statistics to approximate posterior distributions of model parameters, incorporating hierarchical structure of parameters. In this case, a coalescent demographic model is used, with demographic parameters drawn from user-specified priors to simulate patterns of sequence variation. Summary statistics from the simulated datasets are then compared to those from the observed data, and posterior parameter distributions are estimated. MTML-msBayes estimates the number of divergence times across all taxon pairs (Ψ), the mean divergence time in coalescent generations (E(τ)), and the dispersion index of the divergence times (Ω =Var(τ)/E(τ)), defined as the variance-to-mean ratio of divergence times. We used the dispersion index, Ω, to test the hypothesis of simultaneous versus nonsimultaneous divergence of all taxon pairs (i.e., whether all taxon pairs diverged at the same time). The degree to which Ω values deviate from 0.0 (tested using Bayes factors) is interpreted as greater support for non-simultaneous divergence. For further details see Huang et al. (2011).

Due to the computational intractability of analyzing the entire mitogenome, datasets were assembled for each lineage by realigning the previously generated full mitogenomes to extract (i) the ND2 gene and (ii) a ∼4,000 bp subsampled dataset comprising the ND2, cytb, and CO1 genes (hereafter three-gene). These markers were selected based on their prior use in the literature. We further restricted the dataset to taxa with clear evidence of divergence, reducing the initial 39 lineages to 19. We considered lineages to show evidence of divergence if they possessed a pairwise mitogenomic F_ST significantly different from 0 (13/19), or if haplotype sharing (via, for example, introgression, incomplete lineage sorting, or recent divergence) across well-defined OW/NW clades was suspected to have caused a reduced F_ST (6/19). In the six lineages where the latter was the case (Gallinago spp., Numenius phaeopus, Calidris alpina, Larus canus, Pinicola enucleator, and Calcarius lapponicus), this assessment was made based on ND2 haplotype networks, NJ trees, and prior evidence for putative gene flow in these taxa (Zink et al., 1995; Humphries & Winker, 2011; McLaughlin et al., 2020; Spaulding et al., 2022). The cytb and CO1 genes were tested in separate runs of 3 million replicates under identical parameters to the ND2 run.

Following Humphries & Winker (2011) and Topp et al. (2013), we defined a set of Bayesian priors to include as model bounds. We set the upper bound of migration rate, UpperMig, to 10, allowing simulations to assume a maximum of 10 migration events per generation. We kept UpperRec, denoting the upper bound of recombination rate in the study locus, at 0, as recombination is not expected in the avian mitogenome (Berlin, Smith & Ellegren, 2004). The upper bound of ancestral θ (upperAncPopSize), denoting relative ancestral population size relative to current size, was set to 0.5, erring on the side of wider tolerance. Coalescent theory, implemented in msBayes’ hierarchical simulation method, would predict that individuals that contributed ancestry to the current population should be at most one half of the current population. As we do not have estimates of effective Pleistocene population sizes for most avian taxa, we allowed this value to vary as needed. Additionally, the wide variety of ecologically diverse taxa in this study makes it unlikely that a single, uniform effective population size pattern holds true across the dataset (e.g., Maley & Winker, 2010; McLaughlin et al., 2020); thus, an unbounded prior is likely better.

During initial data exploration runs consisting of 500,000 replications, τ values were left unbounded. However, after several exploratory runs, all probable τ values were consistently concentrated between 0 and ∼0.15. For our final analysis we bounded upper τ to 0.2 for all subsequent runs, providing increased coverage within the sampling space. We note that our upper bound prior of migration rate (UpperMig) of 10, allowing simulations to assume a maximum of 10 migration (gene flow) events per generation, reflects a level at which populations are converging (Winker, 2021). We erred on the side of a wider tolerance here because seasonal migration and dispersal vary widely across our lineages.

We performed three million replicates for both datasets. Simple rejection methods (Huang et al., 2011) were used to construct the posterior distributions.

Divergence levels and timing

Pairwise sequence divergence (JC-distance) between lineages for full mitogenome data ranged from 0.0 to 3.82% across all taxa, with an average of 1.39 ± 1.19% (Table 2). Species pairs had a significantly higher average genetic divergence (2.32 ± 1.09%) than population (0.47 ± 0.51%; p < 0.00001), but not subspecies pairs (1.66 ± 1.05; p = 0.16). Subspecies and population pair divergences also differed (p = 0.0003). ND2 divergence depths were generally higher than full mitogenomes, ranging from 0.0 to 4.96% for all taxa, 2.76 ± 1.42% for species, 1.83 ± 1.17% for subspecies, 0.52 ± 0.71% for populations, and an overall average of 1.60 ± 1.46% (Table 2, Table S5). As with full mitogenome data, ND2 divergence levels were significantly different between species and populations and between populations and subspecies (Table 2, Table S5).

Estimated full mitogenome divergence dates between taxon pairs ranged from 69,000 years to 7.96 Mya, with an overall average of 1.95 Mya (Table 3, Table S6). ND2 divergence date estimates, corresponding to the higher divergence rate estimate of Lerner et al. (2011), were lower, ranging from 0 (Gavia stellata, which did not have any ND2-specific mutations) to 1.84 Mya, with an overall average of 0.59 Mya (Table 3, Table S6). Divergence estimates increased only loosely with taxonomic depth (Fig. 2, Table 2; Adjusted R² = 0.4), although taxonomic level was highly significant in predicting divergence (p < 0.001).

Thirteen of 39 taxon pairs were significantly differentiated by pairwise F_ST, of which eight are currently classified as species, two as subspecies, and three as populations (Table 4). Of these, all 13 remained significant after false discovery rate adjustment (Benjamini & Hochberg, 1995; Strimmer, 2008) for both the full mitogenome and ND2 datasets.

Additionally, both datasets in ten pairwise contrasts (Anas crecca, Mareca spp., Calidris alpina, Gallinago spp., Corvus corax, Larus canus, Numenius phaeopus, Anthus rubescens, Pinicola enucleator and Calcarius lapponicus) showed strong NW/OW (Fig. 3) clades with varying amounts of haplotype sharing. Prior work in these taxa also showed this pattern (Zink et al., 1995; Humphries & Winker, 2011; McLaughlin et al., 2020; Spaulding et al., 2022; Doniol-Valcroze et al., 2023). However, in all cases, we were unable to distinguish between incomplete lineage sorting and putative gene flow. While distinguishing between incomplete lineage sorting and gene flow is possible with large-scale transposable element data (e.g., Doronina et al., 2015; Suh, Smeds & Ellegren, 2015), these data were not available to us for any of the 39 taxon pairs investigated.

As F_ST is directly affected by gene flow, the ten taxa that showed strong NW/OW haplotype clades were considered to be substantially differentiated and were included in further analyses with the 13 taxa differentiated by F_ST. This produced a total subset of 19 taxa, as Anas crecca, Mareca spp., Corvus corax, and Anthus rubescens were already significantly differentiated by pairwise F_ST before considering haplotype clades.

Divergence events

Our initial ND2 dataset unambiguously showed the strongest posterior probability support for a mean number of vicariant events (Ψ) of one (i.e., simultaneous divergence), with ∼20.1% of the posterior probability on that value (Fig. 4). However, the data also showed a long, right-hand tail with consistent, relatively high support (∼5%) for all Ψ values greater than one –up to the maximum Ψ of 18 (Table S7). When additional data were added to form the 3-gene dataset (Fig. 4), support for one vicariance event increased to 80.4% of the total, with a corresponding decrease in support for the tail of the distribution (Table S7). However, because estimates of the number of divergence events (Ψ) do not take the temporal spacing of these events into account, it is an unreliable estimator of the number of true divergence events (Overcast, Bagley & Hickerson, 2017).

The dispersion index of divergence times (Ω), a direct test for the simultaneous or nonsimultaneous divergence of all comparisons, where nonsimultaneous divergence is indicated by a Ω significantly different from zero, also differed between these datasets. The ND2 dataset (Fig. 5) showed Ω values significantly different from zero (posterior mean of Ω = 0.0246, Bayes factor = 0 when comparing the probability distribution of Ω ≤ 0.01 and Ω > 0.01), indicating moderate support for more than one (i.e., nonsimultaneous) vicariance event. The 3-gene dataset, in contrast, failed to produce a nonzero Ω-estimate, indicating strong support for simultaneous divergence.

Distribution of avian divergence events in Beringia

Our data indicate that, on average, avian lineages in Beringia differ in divergence depth based on taxonomic level, with species generally having diverged deeper in time than subspecies, and both, on average, having deeper divergences than populations (Fig. 2, Table 2). However, there is also substantial variation within categories, with taxonomic level explaining 45% of the variation in divergence timing. Other taxonomies are available in which different species limits are given than those used here for some taxa (e.g., Withrow, Gibson & Winker, 2025; AviList Core Team, 2025). But these affect a minority of cases in our study and do not materially affect the unsurprising outcome that while mitochondrial divergence levels increase between taxonomic categories (Fig. 2, Table 2), they have considerable variation and overlap between categories and are not reliable indicators of the completion of speciation (cf., Winker, 2021). We make no taxonomic recommendations based on these results.

Our divergence date estimates were more recent than expected for lineages that had been previously studied (Table 1, Fig. 2; e.g., Zink et al., 1995; Humphries & Winker, 2011). Divergence dates show some apparent clumping of events at around 1.3–1.8 Mya, but from 1.0 Mya toward the present divergences appear to be relatively continuous when all taxonomic levels are overlaid (Table 1, Fig. 2). The cluster of recent, population-level divergences from ∼0.0–0.3 Mya likely reflects a pulse of recent, and possibly ongoing divergences among these taxa, albeit with variable levels of gene flow across taxa. This is supported by pairwise F_ST values, which, with the exception of Uria aalge, were nonsignificant for all taxa in this cluster (Table 4). The two population-level comparisons outside of this cluster (Aythya marila, ∼629 Kya and Uria lomvia, ∼1 Mya; Table 1, Fig. 2) showed significant differentiation in F_ST for only one lineage (A. marila).

Divergence events here are strongly biased towards the more recent half of the Pleistocene, and none occur earlier than ∼1.9 Mya (Table 1, Fig. 2). The apparent lack of divergence events during the earliest third of the Pleistocene, followed by what seem to be one or a small number of events in the more recent two thirds, is unexpected. At face value, this result would suggest that avian speciation was only affected by (or at least, more strongly affected by) climatic cycles in the middle-to-late Pleistocene. Climate records (Ruddiman, Raymo & McIntyre, 1986) do not show marked changes in intensity or length of glacial cycles that would explain this apparent response. Also, prior research has suggested that divergences did occur in the early Pleistocene (Zink & Klicka, 2000; Lovette, 2005; Humphries & Winker, 2011; Winker et al., 2023).

One possible explanation for this apparent discrepancy is in the type of sequencing data used here to estimate divergence times versus the historic data and divergence calibrations using such data. The data in Humphries & Winker (2011) and Peters et al. (2014), both generated for similar taxa with equivalent markers, were developed using 1–2x Sanger sequencing, which generates a single consensus sequence from one or two sequencing reads, each deriving from many thousands of amplicons from one set of primers. Alternatively, Illumina sequences like those generated here, are the consensus of ∼100 individual reads, each derived from a single library molecule, which generates the cluster of sequenced molecules. Both methods can produce high-quality data, but the patterns and rates of errors are different (Fridman et al., 2021). High-coverage Illumina data can reduce the incidence of miscalled bases, reducing the frequency of errors in the data and, hence, decreasing overall mitochondrial divergence values versus data with more errors. This was empirically supported in a direct comparison between the prior Sanger and current, high-coverage Illumina ND2 data from the same specimens in the seven shared taxa between this study and that of Humphries & Winker (2011). Nucleotide diversity (π), used in this case as a rough proxy for error rate when considering the difference between two datasets, showed that divergence values calculated from Sanger sequence data were on average 1.67 times higher than values from the Illumina data (Table 5). Although per-sequence error rates are generally higher with Illumina instruments than traditional Sanger sequencing, our filter of discarding base calls with quality scores of <30, combined with our high coverage, creates datasets of higher quality than typical Sanger sequence data (Illumina, 2011). Variation in data quality among sequencing chemistries, instrumentation, and experiments is widely recognized (e.g., Stoler & Nekrutenko, 2021), but we are unable to make direct comparisons between our data and corresponding historic Sanger datasets beyond pointing out that our filter produces data of nearly an order of magnitude higher quality than Sanger data (Q30 vs. ∼Q20; Illumina, 2011) and demonstrating higher variation in the latter (Table 5). Higher-quality data with fewer variants called likely had a major impact on our divergence date estimates (e.g., Table 5).

Additionally, divergence time estimates here were likely impacted by the lack of an appropriate calibration, as the substitution rates from Lerner et al. (2011) were estimated from lower-coverage Roche 454 sequence data and would probably also result in an underestimate of divergence times when applied to data from high-coverage sequences (Luo et al., 2012). Thus, it is likely that the consistently earlier divergence dates generated in this study are an artifact of a mismatched high-coverage dataset (with its correspondingly lower variation than found in data from using older methods, e.g., Table 5) and a calibration rate estimate based on lower-coverage data (Lerner et al., 2011), which would push our divergence estimates closer to the present.

Current methods of dealing with the problem of incorrectly calibrated avian molecular clocks, such as the calibration rates of Arcones, Ponti & Vieites (2021), are not applicable here. Arcones, Ponti & Vieites (2021) calibrated mtDNA divergence rates at the order level by using 622 full mitogenomes representing 33 extant families of birds, with fossil calibration points at depths of 11 to 53 Mya. These rates are suitable for estimating divergence times in deep evolutionary time within Aves, where saturation is high, reciprocal monophyly is assured, and the discrepancy between gene and population divergence is minimal (Arbogast et al., 2002). In our shallow-divergence dataset, Arcones, Ponti & Vieites’s (2021) rates produced divergence dates an order of magnitude older than previous analyses of Beringian birds (Maley & Winker, 2010; Humphries & Winker, 2011) and dates generated here using the more shallowly calibrated (1.9–5.7 Mya), passerine-specific rates estimated by Lerner et al. (2011; Table S6).

Number of avian divergence events in Beringia

For the 18 taxa found to be significantly differentiated (Table 4), one vicariance event would represent complete synchrony and 18 would indicate complete asynchrony of divergence. Here, we find the most support for a single vicariance event, which is biologically surprising, given the staggered taxonomic depths and divergence times of these lineages. Dating calibration issues do not affect this outcome.

One explanation for this might be that successive cycles of isolation and gene flow, most likely driven by glacial-interglacial cycles, have muddied individual vicariance events, such that only the last one or two events are easily detectable. The exact number of vicariance events affecting these lineages is unclear due to the relative difficulty that MTML-msBayes has in estimating Ψ (number of divergence events) among variably temporally clustered divergence events (Overcast, Bagley & Hickerson, 2017). Thus, when considering the frequency of glacial-interglacial cycles (Lisiecki & Raymo, 2005) and the blurring effects on divergence that gene flow has likely had (Oliver Brown, Glenn & Winker, 2025), it might not be possible to sharpen our understanding of divergence clustering in this system.