Comparative phylogeographic inference with genome-wide data from aggregated population pairs

Abstract

Comparing divergences across multiple sister population pairs has been a focus in phylogeography since its inception. Initial approaches used organelle genetic data and involved qualitative comparisons of phylogenetic patterns to evaluate hypotheses of shared and variable evolutionary responses. This endeavor has progressed with coalescent model-based statistical techniques and advances in next-generation sequencing, yet there remains a need for methods that can exploit aggregated genomic-scale data within a unified analytical framework. To this end, we introduce the aggregate joint site frequency spectrum (ajSFS) by validating its use within a hierarchical Bayesian framework through several in silico experiments. Subsequently, we applied our method against two published restriction site–associated DNA marker datasets consisting of eight local replicates of a lamprey species pair and six co-distributed passerine taxon pairs, respectively, with the aim of inferring variability in co-divergence and co-migration histories. We found that the lamprey population pairs exhibited temporal synchrony in both co-divergence and collective secondary contact times, yet an idiosyncratic pattern in secondary migration intensities. In contrast, the bird population pairs displayed thoroughly asynchronous co-divergence histories. Our results demonstrate that the ajSFS can be exploited for complex and flexible co-demographic inference, opening up new possibilities for comparative phylogeography and population genomic studies.

Aggregate site frequency spectrum (aSFS), approximate Bayesian computation (ABC), hierarchical co-demographic modeling, Multi-DICE, population genomic methods, random forest

A long-standing focus in comparative phylogeography is to infer the degree of shared response among replicate sister taxon pairs to emerging or existing landscape barriers for the sake of elucidating processes underlying allopatric or parapatric divergence (Avise et al. 1987; Bermingham and Moritz 1998; Bernatchez and Wilson 1998; Avise 2000; Arbogast and Kenagy 2001; Soltis et al. 2006; Hickerson et al. 2010). It is expected that geographic and climatic changes of major impact could result in similar demographic histories across co-distributed taxa, whereas trait differences (Papadopoulou and Knowles 2016) or even random chance (Rosindell et al. 2012) may cause idiosyncratic demographic trajectories. The interplay between these deterministic and stochastic forces leads to alternative predictions of community assemblage with respect to emergent histories of co-diversification, co-expansion, or co-dispersal. Hence, uncovering the distribution of divergence times across a set of population pairs is of interest to many researchers investigating broader evolutionary and ecological questions (Avise et al. 2016). This endeavor has traditionally relied on qualitatively observing collective patterns of pairwise divergences to circumstantially make inferences about multi-taxa histories (Hoffmann and Baker 2003; Lessios et al. 2003; Alvarado Bremer et al. 2005; Lourie et al. 2005; Goldstien et al. 2006; Richards et al. 2007). Likewise, it is commonplace to approximate whether geographic boundaries that divide structured populations coincide among multiple species or taxa, with the assumption that general similarities suggest concordant demography in joint historical response to environmental features (Taberlet et al. 1998; Zink 2002; Pastorini et al. 2003; Barber et al. 2006; Maliouchenko et al. 2007).

However, the advancement of inferential methods has allowed researchers to investigate historical processes using probabilistic models that offer parameter estimation. For comparative inference, such analyses can be conducted independently per constituent dataset, with results compared post hoc to qualitatively assess historical congruence (Carstens et al. 2004; Fedorov et al. 2008; Bacon et al. 2015; Edwards et al. 2016). Alternatively, multiple demographic histories can be modeled together simultaneously within a single unified analysis while accounting for individual variance across taxa in nonfocal nuisance parameters (Edwards and Beerli 2000; Carstens et al. 2016; Satler and Carstens 2016). One such approach involves hierarchical co-demographic modeling, which allows explicit hypothesis testing among taxa (Hickerson et al. 2007; Chan et al. 2014; Xue and Hickerson 2015, 2017) and increases accuracy due to the statistical pooling power or borrowing strength derived from globally modeling exchangeable datasets (Congdon 2001; Gelman et al. 2003; Beaumont 2010). Consequently, uncertainty measurements from this tactic consider the stochasticity for each constituent dataset in conjunction, in contrast to post hoc qualitative comparisons on individual results.

Hierarchical co-demographic modeling for sister taxon pairs is available in the software programs msBayes (Hickerson et al. 2007, 2014; Huang et al. 2011) and dpp-msbayes (Oaks 2014), which use an approximate Bayesian computation (ABC) inferential framework. In brief, ABC employs Monte Carlo simulations to construct a reference table of random parameter draws from prior distributions coupled with associated summary statistic vectors, which is subsequently applied to empirical data via a rejection algorithm that minimizes Euclidean distance of summary statistics to approximate posterior distributions of tested models and/or parameters (Beaumont et al. 2002; Bertorelle et al. 2010; Csilléry et al. 2010). The variant of ABC deployed by msBayes and dpp-msbayes assumes a hierarchical co-demographic model (hABC), and exploits distribution summaries of standard population genetic summary statistics calculated individually per population pair, to infer hyperparameters and parameter summaries that describe the variation in divergence times among taxon pairs. However, these methods have restricted modeling capabilities with respect to specifying historical events (e.g., complex co-migration trajectories) as well as an inflexible hyperparameterization strategy (e.g., cannot estimate proportion of taxa membership within a hypothesized single synchronous co-divergence pulse, or congruence in parameters besides co-divergence timing). More importantly, this approach is computationally constricted to single-locus data such as mitochondrial barcodes, or a small number of nuclear loci (Carnaval et al. 2009; Stone et al. 2012; Smith et al. 2014), such that statistical power and accuracy are strongly compromised (Overcast et al. 2017).

With the use of next-generation sequencing becoming more commonplace as a result of reducing costs, there is a growing need and utility for extending hierarchical co-demographic modeling to genomic-scale applications (Xue and Hickerson 2015; Oaks 2018). Such implementation could also improve the ability to explore more complex models, including incorporation of various migration and size change parameters as well as quantifying degree of congruence among multiple demographic parameters. The recently developed aggregate site frequency spectrum (aSFS), which permits comparative size change models of independent single populations coupled with genome-wide data and can be used in the R package Multi-DICE (Xue and Hickerson 2015, 2017; Prates et al. 2016; Mastretta-Yanes et al. 2018; Potter et al. 2018), opens the possibility of a large-scale data summarization that contains signal of co-demography across multiple population pairs. Notably, previous implementations of Multi-DICE supplemented hABC results with random forest, a machine learning algorithm that is experiencing increased attention in population genetics due to its ability to predict model assignment for empirical data from training on Monte Carlo simulations (Pudlo et al. 2016; Schrider and Kern 2018).

To this end, we describe, evaluate, and deploy a hierarchical co-demographic model of population pairs, as well as a generalized approach of comprehensively modeling individual population pair divergence and migration scenarios that can accompany our hierarchical model. We also introduce the data summarization associated with our co-demographic model, the aggregate joint site frequency spectrum (ajSFS). The ajSFS is similar in design to the aSFS with respect to its efficient aggregation of independent genomic-scale datasets, which can be effectively exploited by Monte Carlo simulations to achieve statistical inference of hyperparameters that quantify the degree of congruence in co-demographic histories. However, an important distinction from the aSFS is that the ajSFS is specifically intended for hierarchical models of aggregate population pairs rather than single populations, similar to msBayes and dpp-msbayes yet with more flexibility. To establish proof-of-concept, we conduct simulation experiments to assess performance using both random forest regression given our hierarchical model (hRF) as well as hABC, and then reanalyze two previously published restriction site–associated DNA marker (RAD) datasets (Oswald et al. 2017; Rougemont et al. 2017).

Both studies associated with the empirical data employed model selection of various divergence with migration scenarios to address a central inquiry of a comparative nature. The first comparative dataset consists of a species pair between river lamprey (Lampetra fluviatilis), which are anadromous and parasitic, and its freshwater-exclusive and nonparasitic counterpart brook lamprey (Lampetra planeri), replicated in eight sampling localities throughout Western France (Rougemont et al. 2017). For the original investigation, Rougemont et al. (2017) focused on how environmental connectivity affected the lamprey species pair's demographic history, particularly with respect to migration, across these putatively natural replicates. Specifically, four population pairs considered to be connected (CPPs) all had secondary contact as the model with highest likelihood, whereas ancient migration was mostly inferred for the remaining four population pairs deemed disconnected (DPPs). The coarse trends of several relative estimates (i.e., proportion of two parameters within a population pair) likewise contrasted CPPs from DPPs. The second comparative dataset to which we apply our method includes six passerine taxon pairs co-distributed along a barrier separating the Marañón Valley and Tumbesian region in Peru (Oswald et al. 2017). In the initial paper, Oswald et al. (2017) focused on how a shared phylogeographic break broadly affected bird pairs co-distributed in seasonally dry forests. Their results suggested asymmetric gene flow for every taxon pair following an overall asynchronous co-divergence. For our approach, we are interested in extending the findings for these two datasets by directly evaluating temporal synchrony among co-divergence times as well as variability within other demographic parameters related to postdivergence migration and population size change.

Materials and Methods

HIERARCHICAL CO-DEMOGRAPHIC MODEL OF POPULATION PAIRS

We present a hierarchical co-demographic model extending that which is described in Xue and Hickerson (2017) to accommodate population pairs instead of single populations, multiple historical events per individual taxonomic unit, and various co-demographic hyperparameters (Fig. 1; Table S1). As with Xue and Hickerson (2017), one level of the hierarchical model consists of n individual population pairs that are each treated as an independent iteration with its own demographic parameters (therefore no assumptions about shared ancestry are enforced), whereas on another level are co-demographic hyperparameters that explain variability among analogous demographic parameters from the constituent population pairs. Population indices correspond along all n number of incorporated population pairs (e.g., population “0” in one population pair is associated with population “0” from all other population pairs). Demographic parameter symbols notate the corresponding set of taxon-specific demographic parameters across population pairs, with all of the constituent elements, which allude arbitrarily to one of the n individual taxa, identified by numeric subscripts, for example, ${\tau }^1=\{{\tau }_1^1,{\tau }_2^1,\dots ,{\tau }_n^1\}$⁠. To estimate variability in taxon-specific responses, we defer to the hyperparameters defined in Xue and Hickerson (2017), specifically ζ_T, S_T, and ψ. In general, these could be applied to any set of analogous demographic parameter values among taxa, thus parameter symbols are subscripted to distinguish co-demographic hyperparameters (e.g., ζ_T,τ1 and ψ_τ1 for τ¹). Accordingly, a demographic parameter symbol with an s subscript denotes a vector of values that coincide with ψ shared pulses of co-demography, with the vector elements identified by numeric subscripts that reference the specific pulse to which they are assigned (e.g., ${\tau }_s^1=\{{\tau }_{s,1}^1,{\tau }_{s,2}^1,\dots ,\tau {{}_{s,\psi }^1}_{,}\}$ for synchronous co-divergence timings). Likewise, an i subscript represents the vector of idiosyncratic draws that are each unique in value, with numeric subscripts arbitrarily pertaining to individual population pairs (e.g., ${\tau }_i^1=\{{\tau }_{i,n-S_T+1}^1,{\tau }_{i,n-S_T+2}^1,\dots ,{\tau }_{i,n}^1\}$ ).

Figure 1

Hierarchical co-demographic model. An example model with eight population pairs and half of these pairs within a single temporally synchronous co-divergence pulse. The pulse of synchronous co-divergence occurred at ${\tau }_s^1$⁠, whereas the remaining four population pairs experienced idiosyncratic co-divergence at independent times of ${\tau }_{i,5}^1$⁠, ${\tau }_{i,6}^1$⁠, ${\tau }_{i,7}^1$⁠, and ${\tau }_{i,8}^1$⁠, respectively. For visual simplicity, the individual model used here does not have migration nor size change.

For this study, we focused on evaluating hyperparameters, in particular the total proportion of the whole set of taxon-specific demographic parameters that are congruent in value, ζ_T (Table S1). The specific demographic parameter values per population pair (Table 1) were therefore regarded as nuisance parameters, which are not of direct interest yet important to account for noise and uncertainty stemming from population differences. In particular, we emphasized ζ_T with respect to co-divergence time, which is formally designated as ζ_T,τ1 but hereby referred to as ζ for convenience. Within our simulation-based inferential approach, we treated ζ as a freely varying hyperparameter that governs variation in τ¹ values across population pairs a priori. Conversely, to efficiently derive values of ζ_T for other demographic parameters, we deployed a novel post hoc protocol that instead calculates these hyperparameter values from freely varying taxon-specific parameter values (e.g., ζ_T,τ2 is determined from random draws of τ² per population pair). This was accomplished by counting how many individual values were within a predefined distance from the median for the total set of draws per simulation, with zero counts assigned to ζ_T = 1/n because this value represents complete idiosyncrasy and thus is inclusive of cases whereby ζ_T = 0.0. Notably, the predefined distance from the median approximates our use of β to condition the ${\tau }_i^1$ prior distribution, for example, assuming this value is set to 0.05 natural log units, then the degree of clumping up to 0.05 natural log units around the central tendency within the range of taxon-specific draws determines the value of post hoc ζ_T, similar to β = 0.05 natural log units. Values for a post hoc ζ_T hyperparameter then can be combined across simulations to construct hABC hyperprior and hyperposterior distributions. This simple post hoc conversion process demonstrates that hyperparameters, which provide a convenient interpretation of parameter variability across population pairs, can still be extracted when not incorporated into the simulation design. Such a possibility allows further exploration of the hierarchical co-demographic model in a quick and efficient manner, potentially eliminating the need for heavier initialization of hyperparameters, which must be accommodated by more simulations and thus increasingly expensive computation.

Table 1

Index of parameters for individual population pair model

Symbol	Description
τ1	Time of divergence
τ2	Time of instantaneous change in migration rates and effective population sizes
M  1-0:1	Migration rate from tip population “0” to “1” forward-in-time during recent phase
M  1-1:0	Migration rate from tip population “1” to “0” forward-in-time during recent phase
M  2-0:1	Migration rate from tip population “0” to “1” forward-in-time during ancient phase
M  2-1:0	Migration rate from tip population “1” to “0” forward-in-time during ancient phase
ε0	Size change magnitude for tip population “0”
ε1	Size change magnitude for tip population “0”
N  0	Present-day effective population size for tip population “0”
N  1	Present-day effective population size for tip population “1”
u	Mutation rate for entire population pair

Additionally, we focused on estimation of the following parameter summaries, which describe the distribution of parameter values across population pairs and can supplement interpretation of inferred hyperparameter values: (1) timing for a single pulse of synchronous co-divergence (formally designated as ${\tau }_{s,1}^1$⁠, but hereby referred to as ${\tau }_s^1$ for convenience since ${\tau }_s^1=\{{\tau }_{s,1}^1\}$ for only a single pulse); (2) mean absolute deviation from the median (Δ; e.g., $\mathrm{\Delta }({\tau }^1)=\frac{1}{n}{\sum }_{j=1}^n|{\tau }_j^1-\mathrm{median}({\tau }^1)|$ ); (3) dispersion index (Ω; i.e., variance/arithmetic mean); and (4) arithmetic mean (Table S1). Notably, although a priori governing hyperparameters are prescriptive determinants for the full distribution of parameter values, these are distribution summaries that are descriptive in nature, and in fact derived similarly from the same set of simulated values as our post hoc ζ_T hyperparameters. Importantly, our usage of Δ is a new implementation within this context of co-demographic modeling. As a deviation statistic, it provides a standard measure of statistical dispersion that benefits from conveying information in an intuitive manner, specifically the average difference from the central tendency within a set of values such that a range of 2 × Δ contains approximately half of the population pairs, as well as being more robust to low sampling than Ω due to relying less on normality assumptions (Herrey 1965; Pham-Gia and Hung 2001; Leys et al. 2013).

Importantly, we elaborate upon two caveats of our hierarchical co-demographic model with respect to biological interpretation and our assumption of population pair independence within the Supporting Information.

DEMOGRAPHIC MODEL OF INDIVIDUAL POPULATION PAIRS

To build our hierarchical co-demographic model, we used a generalized parameterization for individual population pairs that can efficiently explore the data independently for each population pair and can be further informed to construct the individual taxon level of our hierarchical scheme. Several commonly tested divergence scenarios are nested within our comprehensive individual-level model (Fig. 2A), including isolation with continuous migration through to the present (IM; Fig. 2B), isolation with continuous migration that ceases before present (i.e., ancient migration [AM], often associated with parapatric divergence; Fig. 2C), strict isolation with a period of no migration followed by continuous migration through to the present (i.e., secondary contact [SC]; Fig. 2D), and strict isolation (SI) with absolutely no migration (often associated with allopatric divergence; Fig. 2E). Such isolation/migration syndromes have traditionally been treated as competing discretized models, which are tested by performing model selection through comparing relative likelihoods or model posterior probabilities (Pelletier and Carstens 2014; Roux et al. 2016; Jackson et al. 2017; Oswald et al. 2017; Rougemont et al. 2017). Our global model instead captures these syndromes under different regions of parameter space (Table 1), an approach we favor because a wide span of overall divergence scenarios can be quickly assessed independently across multiple taxa under a single analysis. Importantly, our generalized population pair model can discover trajectories that are intermediate between the aforementioned traditional population pair scenarios (e.g., mostly consistent with SC but with a small yet detectable amount of early phase migration, which would not be captured if selecting among discrete models). Furthermore, in cases where similar patterns are elucidated throughout the individual taxa, the flexible parameterizations can be easily modified to increase model specificity, thus increasing identifiability, while still allowing population pair datasets and associated individual parameterizations to be treated as exchangeable within our downstream hierarchical co-demographic model construction (Schlüter et al. 1997; Gelman et al. 2003; Teh et al. 2005). As an example, preliminary analyses with our generalized parameterizations could reveal that several individual datasets resemble IM while the remaining are indicative of AM, in which case it would be desirable to restrict the migration intensity parameters such that only the model space inclusive of either IM or AM is permitted within a corresponding hierarchical co-demographic approach.

Figure 2

Individual population pair model. (A) Generalized model with two phases of asymmetric migration, as well as instantaneous size changes occurring independently between tip populations and simultaneously with shift in migration rates. This model subsumes other commonly used divergence models, such as (B), (C), (D), and (E), under specific parameter combinations. Demographic parameter symbols are annotated onto this model here and described in Table 1. (B) IM model, which the generalized model converges toward when approximating M^1-0:1 = M^2-0:1 and M^1-1:0 = M^2-1:0 (or alternatively τ¹ = τ²). (C) AM model, which the generalized model converges toward when approaching M^1-0:1 = 0.0 and M^1-1:0 = 0.0. (D) SC model, which the generalized model converges toward when near M^2-0:1 = 0.0 and M^2-1:0 = 0.0. (E) SI model, which the generalized model converges toward when resembling M^1-0:1 = 0.0, M^1-1:0 = 0.0, M^2-0:1 = 0.0, and M^2-1:0 = 0.0. Moreover, this model does not contain size change to further illustrate that the generalized model becomes consistent with stable size if ε parameters have values close to 1.0 (or alternatively τ¹ = τ²).

To this end, our generalized population pair model (Fig. 2A; Table 1) treats the two “tip” lineages as sister panmictic populations derived from their most recent common ancestral population. Deep population structure for the ancestral or either of the descendent populations is assumed to be absent, with tip populations indexed arbitrarily as “0” and “1” to distinguish parameters. Specifically, population “0” and population “1” are subject to respective present-day effective population sizes (N_E), N⁰ and N¹, in number of haploids. Divergence of the ancestral population occurs at time τ¹ and the descendent tip populations experience two phases of asymmetric continuous migration as well as instantaneous size changes. Both shifts in migration rates as well as size changes are simultaneously associated with time τ². These two time parameters are scaled by number of generations. The four different asymmetric migration intensity parameters (M) are notated as M^x-y:z, with x symbolizing recent (1) or ancient (2) phases and y:z designating source and sink population indices, respectively, from a forward-in-time perspective (i.e., M^1-0:1, M^1-1:0, M^2-0:1, and M^2-1:0). Migration intensities are rates of effective haploid migrants per generation, that is, scaled by N_E × m, where m is the generational probability of per haploid individual migration. The two distinct population size change magnitude parameters (ε) are affixed with the aforementioned population indices (i.e., ε⁰ and ε¹), with population-specific size change magnitudes scaled by the ratio of N_E before size change to N_E after size change from a forward-in-time context (i.e., past:future), such that values <1 denote expansion and >1 signify contraction. Importantly, mutation rate (μ), in number of substitutions per site per generation, can also be incorporated to exploit sequence length information (l) for calibration of N_E parameters (via the population genetics parameter θ = 4N_Eμl), which then inherently calibrates time parameters as assumed by coalescent theory. If μ is omitted, monomorphic sites are implicitly ignored and thus parameter estimates are conditioned on given N_E, which may be informed by a priori information. In such case, N_E ought to still be governed by a prior distribution to account for uncertainty as well as allow for heterogeneity in N_E among population pairs in the downstream hierarchical co-demographic model.

AGGREGATE JOINT SITE FREQUENCY SPECTRUM

To enable inference under supervised machine learning or ABC, we introduce a summarization appropriate for genome-wide single nucleotide polymorphism (SNP) data that extends the aSFS summary statistic vector described in Xue and Hickerson (2015), the ajSFS (Fig. 3). To construct the ajSFS, a joint site frequency spectrum (joint-SFS) with relative SNP frequencies (i.e., proportion of total SNPs rather than SNP count) is calculated for every population pair. As aforementioned, monomorphic sites can also be optionally included to better calibrate absolute parameter estimates as long as the μ parameter (Table 1) is used during simulation, although this count needs to be converted into the proportion to SNPs within the joint-SFS due to the use of relative frequencies. Importantly, all constituent joint-SFS are projected to equivalent sampling for number of individuals, after which they are ideally comparable in sampling of genomic sites as well. With respect to number of sites, although slight deviations are negligible, major differences in sampling among joint-SFS can be addressed by conservatively assuming the lowest amount for the entire aggregate dataset during simulation. After being independently built, the multiple joint-SFS are collated and bins are rearranged in descending order per joint allele frequency class, thus forming a new summary statistic vector, the ajSFS. As with the aSFS, the intuition behind this construction is to efficiently combine independent datasets such that important aspects of the aggregate data are directly extracted, in this case the relative distribution of each joint allele frequency class across population pairs, whereby similarities and deviations between constituent datasets are captured by hyperparameters and associated parameter summaries. To clarify, this exchangeability assumption is the reason that relative SNP frequencies and equivalent sampling for the multiple joint-SFS are necessary to derive the ajSFS. Importantly, although information about the identity of individual population pairs is discarded, the ajSFS heavily increases computational tractability for simulation-based approaches while exploiting valuable multi-taxa insight, allowing assemblage-level hypotheses to be broadly tested while complementary methods could be used to ascertain specific individual responses. To better illustrate our ajSFS protocol, we provide an explicit example within the Supporting Information.

Figure 3

Construction protocol for ajSFS. (A) Example of eight independent joint-SFS with relative SNP frequency values, visualized here as heatmaps, projected to the same haploid sampling. (B) The first joint allele frequency class across the eight individual datasets is combined, such that the eight separate cells are now brought together. (C) The eight bins are reordered in descending fashion. (D) The process repeats for all joint allele frequency classes. (E) The final ajSFS, consisting of values from all eight of the original joint-SFS and visualized here as a one-dimensional vector, is constructed.

SIMULATING REFERENCE TABLES AND TEST DATA FROM PRIOR AND HYPERPRIOR DISTRIBUTIONS

To demonstrate our hierarchical co-demographic modeling approach with the ajSFS, we conducted numerous simulations (Excoffier et al. 2013; Xue and Hickerson 2017) to be exploited by in silico experiments as well as empirical analyses of two publicly available datasets, specifically reduced representation genomic-level SNPs from eight replicates of a lamprey species pair and six co-distributed taxon pairs of passerine birds, respectively, within a Bayesian inferential framework. Specifically, simulations were organized into reference tables, which consist of independent entries of a summary statistic vector generated from random parameter draws from user-specified prior distributions, and can be employed for machine learning or rejection sampling to discover the most likely parameter values that fit an observed dataset. Our first group of simulations was produced to test the performance of our generalized population pair modeling scheme within the hierarchical framework. To this end, we constructed three ajSFS reference tables that were governed by this generalized model, an IM scenario, and an SI scenario, respectively. Accordingly, pseudo-observed datasets (PODs) were simulated under nine different co-demographic scenarios, including those three applied for the reference tables, to explore model misspecification (Table S2). For our other batch of ajSFS simulations, we produced reference tables that matched the sampling specifications of the two empirical datasets, allowing both in silico experiments that further acted as power analyses relevant to the empirical hyperparameter and parameter summary inferences, as well as statistical inference on the empirical data. Furthermore, to better inform the underlying hierarchical models, we created joint-SFS reference tables under the generalized population pair model to conduct preliminary empirical analyses independently across the individual taxa (Tables S3 and S4).

For the ajSFS simulations associated with the empirical analyses, we employed alternative hierarchical modeling strategies (Tables S3 and S4): (1) weakly informed parameterization—wide taxon-specific prior distribution ranges that matched those in the preliminary analyses on the individual population pairs; (2) SC parameterization (lamprey analyses only)—enforced SC scenario with more constraint in other taxon-specific prior distributions, as indicated by the preliminary results; (3) fixed ${\tau }_s^1$ parameterization (bird analyses only)—${\tau }_s^1$ fixed in value by leveraging a previous hABC result, thus inference was conditional on a hypothesized synchronous time of co-divergence while still allowing ${\tau }_i^1$ to freely vary; and (4) SC and fixed ${\tau }_s^1$ parameterization—SC parameterization, which was supported by the preliminary results among the birds as well, with simultaneous fixing of ${\tau }_s^1$ .

To explore certain properties of the ajSFS, we subsequently transformed the initial reference tables into alternative sampling specifications. For the lamprey simulations, the ajSFS with 20 haploids per tip population constituted a large composite vector, therefore we down-projected these data to 10 haploids per tip population to discern effects from this sampling difference between 10 versus 20 individuals. The down-projection followed the same procedure applied by the Python module δaδi, whereby the number of individuals within the joint-SFS can be reduced through the hypergeometric distribution, which integrates over every subsampling combination (Gutenkunst et al. 2009). For the bird simulations, we deployed three variants of the ajSFS to assess the effects of incorporating sequence length information: (1) ajSFS configuration omitting monomorphic sites, such that the ajSFS was scaled to relative proportions of SNPs only, as with the lamprey ajSFS; (2) ajSFS configuration including monomorphic sites, with the monomorphic frequency class scaled to the proportion of invariants to SNPs per population pair (remaining ajSFS bins still in relative SNP proportions), as well as subsequent normalization of all bins to account for the orders of magnitude difference in scale and sampling variance between monomorphic and polymorphic bins; and (3) ajSFS configuration involving the same normalization procedure as (2) but with monomorphic bins removed to isolate any effect from this normalization. Normalizing ajSFS simulations entailed a standard protocol of subtracting the mean value and then dividing by the variance for each individual element, with mean and variance calculated per respective element from the whole reference table.

Every ajSFS reference table across in silico experiments consisted of 1.2 million ajSFS simulations, and collectively totaled to 18 separate sets of simulations from: three modeling schemes, which were individually tested against nine sets of 100 ajSFS PODs with varying parameterizations, for the first experiment (Table S2); three parameterizations (Table S3) under two sampling projections for number of individuals (i.e., six reference tables) in the case of the lamprey experiment; and three parameterizations (Table S4) with three data constructions involving monomorphic sites and normalization (i.e., nine reference tables) under the bird experiment. Additionally, each of the two individual population pair reference tables for our preliminary empirical ABC analyses, that is, one for the lampreys (Table S3) and one for the birds (Table S4), were built with 500,000 joint-SFS vectors according to the same simulation specifications (i.e., simulation model, number of individuals, and amount of genetic data) as the corresponding ajSFS reference tables for the lamprey and bird data, respectively. To clarify, ABC reference tables were constructed only given 20 haploid samples for the lamprey data, and only without monomorphic sites nor normalization for the bird data, because these both ostensibly offered the best inferential resolution as determined by the results of the simulation experiments.

We discuss in greater detail the data specifications and model parameterizations of our simulations within the Supporting Information.

EXPERIMENTS ON CO-DEMOGRAPHIC INFERENCE USING SIMULATION REFERENCE TABLES

To perform in silico experiments for our hierarchical co-demographic methodology, we explored bias and accuracy under hRF hyperparameter inference of ζ (Liaw and Wiener 2002) as well as hABC hyperparameter and parameter summary estimation (Table S1) under simple rejection sampling with mean, median, and mode point estimates of the posterior. Specifically, we tested the application of our generalized modeling scheme for population pairs, parameterization options that sequentially increased the degree of a priori information, and ajSFS construction schemes that varied in sampling of individuals, inclusion of monomorphic sites, and normalization of ajSFS bins. To accomplish this, in our first in silico experiment that tested our generalized population pair scenario against other models, every single POD underwent separate inferences against each of the three simulated reference tables. Here, hABC inference included hyperparameter estimation of ζ and parameter summary estimation of Δ(τ¹), Ω(τ¹), and $\overline{x}$ (τ¹). For the other ajSFS reference tables that were based on the two empirical specifications, we conducted “leave-one-out” cross-validation, which entails iteratively treating a single random simulation as a POD and performing inference with the remaining reference table (Csilléry et al. 2012). In this case, hABC inference included estimates of ζ, other ζ_T hyperparameters derived post hoc to the taxon-specific draws, ${\tau }_s^1$⁠, and parameter summaries (i.e., Δ, Ω, and $\overline{x}$ ) for τ² and several migration parameters. Additionally, we examined robustness to error in assuming SC or the fixed value of ${\tau }_s^1$ within the context of our empirical applications, by selecting a single set of PODs from the second reference table in sequence (i.e., generated with the second-in-sequence parameterization option) but using the third and final reference table in sequence (i.e., simulated with the SC and fixed ${\tau }_s^1$ parameterization for both the lampreys and the birds) for hRF and hABC inference of ζ. Notably, the difference in the second parameterization option between the lamprey and bird analyses allowed us to isolate the effect of using solely an SC parameterization and fixed ${\tau }_s^1$ parameterization, respectively, therefore allowing an independent assessment of each tactic. Across in silico experiments, estimated values were compared against true values by calculating Pearson's r correlation and root mean squared error (RMSE). Additionally, hRF confusion matrices, which required rounding estimates to the nearest discrete interval, and the mean of the hABC hyperposterior distributions per true value were produced for ζ hyperparameter estimation. We further describe this inferential procedure for the simulation experiments within the Supporting Information.

COMPARATIVE INFERENCE ON EMPIRICAL AGGREGATED DATASETS OF LAMPREY AND BIRD POPULATION PAIRS USING SIMULATION REFERENCE TABLES

To empirically demonstrate our hierarchical modeling approach with the ajSFS, we applied our method against two publicly available genome-wide SNP datasets to detect co-demographic pulses of divergence with migration. We selected these two comparative datasets because they starkly contrast in terms of underlying taxonomic relationship, biogeographic question, and sampling, thus span a broad spectrum of empirical applications. The first was collected from a lamprey species pair, consisting of the anadromous, parasitic river lamprey (Lampetra fluviatilis) and freshwater brook, nonparasitic lamprey (Lampetra planeri), across eight sampling localities throughout western France (Rougemont et al. 2017). We adhered to the population and pair assignments as well as groupings of localities based on environmental connectivity of Rougemont et al. (2017), specifically L. fluviatilis and L. planeri populations were designated population “0” and population “1,” respectively, following the structure of the data files, and the connected population pairs (CPPs) included the population pairs AA, BET, OIR, and RIS (from sampling locations Aa, Bethune, Oir, and Risle River, respectively), whereas BRE, CEN, ODO, and SAU (from sampling locations Bresle, Loire-Cens, Odon River, and Garonne-Saucats, respectively) were part of the disconnected population pairs (DPPs). The second dataset was derived from ddRAD sequences for six taxon pairs of dry forest birds from Peru: Arremon abeillei nigriceps—A. a. abeillei, Campylorhynchus fasciatus fasciatus—C. f. pallescens, Melanopareia maranonica—M. elegans, Mimus longicaudatus maranonicus—M. l. longicaudatus, Saltator striatipectus peruvianus—S. s. immaculatus, and Thamnophilus bernardi shumbae—T. b. bernardi (Oswald et al. 2017). The first taxon listed for each of these bird population pairs was sampled within the Marañón Valley, whereas the second listed taxon was located within the Tumbesian region, and these taxa were respectively allocated to population “0” and population “1” following the structure of the data files from Oswald et al. (2017).

We first conducted preliminary analyses on the individual population pair datasets independently to explore appropriate parameterization for co-demographic modeling. This was accomplished via ABC under the simple rejection algorithm (Csilléry et al. 2012), which allowed uncertainty measures to be easily obtained. Moreover, a single reference table could be repeatedly applied per joint-SFS across an entire comparative dataset; hence, we have only one set of simulations each for the lampreys and birds yet separate posteriors for every respective population pair. This exploratory step then can be efficiently reiterated many times over to discover broadly applicable prior distributions among constituent population pairs. Subsequently, hierarchical co-demographic estimates were performed (Xue and Hickerson 2017) with both hRF (Liaw and Wiener 2002) and hABC (Csilléry et al. 2012), following the inferential protocol for the corresponding simulation experiments. Notably, within this context of empirical analysis, our manner of informing prior distributions from a previous analysis on the data constitutes an empirical Bayes approach (Morris 1983; Chaudhuri et al. 2018), and our succession of parameterization options that update prior distributions resembles a simplified sequential ABC method (Sisson et al. 2007). To heuristically and visually assess goodness-of-fit of our hierarchical model, we used principal component analysis (PCA; Cornuet et al. 2010; Berlin et al. 2011; Sjödin et al. 2012) to execute both a predictive check sampled from the prior and a predictive check sampled from the posterior (Rubin and Stern 1994).

We present additional details of our empirical analyses within the Supporting Information.

Results

SIMULATION EXPERIMENTS OF HIERARCHICAL CO-DEMOGRAPHIC INFERENCE

Our first simulation experiment, which investigated using our generalized population pair approach against other migration scenarios within our hierarchical co-demographic model (Tables 2 and S5–S7), found that inference under our generalized model has substantially greater accuracy as the true model underlying the PODS increases in migration complexity. Conversely, employing an SI scenario for inference was better in cases whereby migration was minimal for the true model, whereas applying an IM scenario for inference was seemingly an intermediate between these two extremes. Unsurprisingly then, performance was optimized when deploying the correct model for inference (i.e., model generating the reference table was the same as that simulating the PODS; r = 0.652–0.781 for hRF inference of ζ). Similarly, slight model misspecifications were well tolerated, for example, estimation on PODs simulated under AM given the SI reference table (r = 0.684 for hRF inference of ζ), or estimation on PODs generated under strict SC versus SC with minor ancient phase migration (r = 0.623 and 0.685, respectively, for hRF inference of ζ given the reference table under the generalized model). However, major model misspecifications drastically reduced accuracy, such as the SI reference table unable to handle the six POD sets that were not simulated under SI or the two AM scenarios, or the IM reference table behaving more consistently but still suffering against PODs generated with more dynamic migration trajectories as well as SI PODs.

Table 2

Hierarchical co-demographic model misspecification of three reference tables against nine pseudo-observed datasets testing inference of ζ with random forest regression

	GM reference table		IM reference table		SI reference table
	r	RMSE	r	RMSE	r	RMSE
GM	0.652	0.213	0.612	0.231	0.465	0.247
IM	0.490	0.262	0.708	0.213	0.344	0.285
SI	0.531	0.261	0.703	0.226	0.781	0.188
AM	0.581	0.249	0.521	0.262	0.684	0.225
AX	0.517	0.249	0.434	0.268	0.503	0.252
SC	0.623	0.239	0.557	0.261	0.388	0.277
SX	0.685	0.212	0.553	0.243	0.314	0.281
XM	0.640	0.226	0.526	0.249	0.415	0.260
XX	0.511	0.241	0.453	0.249	0.399	0.263

Order of models follows the same order they were introduced in the Supporting Information; GM—population pair model with generalized parameterization of migration incorporated within hierarchical co-demographic model; IM—isolation with migration population pair model incorporated within hierarchical co-demographic model; SI—strict isolation population pair model incorporated within hierarchical co-demographic model; AM—ancient migration population pair model incorporated within hierarchical co-demographic model; AX—AM with minor assumption violation in recent phase migration; SC—secondary contact population pair model incorporated within hierarchical co-demographic model; SX—SC with minor assumption violation in ancient phase migration; XM—hierarchical co-demographic model that mixes IM, AM, and SC; XX—hierarchical co-demographic model that mixes GM, AX, and SX.

For the leave-one-out cross-validations under the empirical data specifications, our experiments strongly suggest that our weakly informed parameterization is a challenging modeling approach, whereas the priors informed through empirical Bayes led to much greater accuracy. Specifically, the SC parameterization achieved powerful inference under the lamprey specifications (Fig. S1; Tables 3 and S8–S10), and the fixed ${\tau }_s^1$ parameterization likewise offered substantial, albeit less so, gains in sensitivity for the bird analyses (Fig. S2; Tables 4 and S11–S13). Moreover, combining the SC parameterization with a fixed ${\tau }_s^1$ value yielded additional increase in performance under both the lamprey and the bird sampling configurations, illustrating that either setting M² parameters to 0.0 or fixing ${\tau }_s^1$ to a single value separately boosts inferential resolution. Furthermore, we do not see any ostensible changes in the efficacy of the SC and fixed ${\tau }_s^1$ parameterization when comparing these results to our misspecification test that sampled PODs instead from the preceding reference table generated under the less informed parameterization, that is, under the lamprey SC parameterization (Tables S14-S16) and bird fixed ${\tau }_s^1$ parameterization (Tables S17–S19), respectively. However, despite making the SC and fixed ${\tau }_s^1$ parameterization equivalent, we witness quite greater accuracy for the lamprey-inspired cross-validation (e.g., r = 0.965 for hRF estimation of ζ given the 10 haploids data sampling). This is expected given that the lamprey-based simulations contained more population pairs and were much more informed with respect to the prior distributions of τ² and ε parameters than for the bird-related simulations.

Table 3

Cross-validation of hierarchical co-demographic modeling under lamprey data specifications given 10 haploids for hRF and hABC mean point estimates

	Weakly informed parameterization		SC parameterization		SC and fixed ${\tau }_s^1$ parameterization
	r	RMSE	r	RMSE	r	RMSE
	hRF hyperparameter estimation
ζ	0.491	0.255	0.824	0.163	0.965	0.076
	hABC hyperparameter estimation/model selection
ζ	0.492	0.254	0.773	0.182	0.933	0.109
${\zeta }_{{\tau }_2}$	0.302	0.183	0.852	0.110	0.868	0.102
ζM  1-0:1	0.737	0.196	0.424	0.091	0.601	0.088
	hABC parameter summary estimation
${\tau }_s^1$	0.301	291,790	0.541	237,101	N/A	N/A
Δ(τ1)	0.342	94,872	0.872	47,614	0.922	38,797
Ω(τ1)	0.361	111,858	0.867	54,368	0.923	36,098
$\overline{x}$ (τ1)	0.668	131,999	0.962	47,946	0.902	43,535
Δ(τ2)	0.595	50,942	0.872	14,828	0.923	11,469
Ω(τ2)	0.270	104,351	0.845	21,952	0.901	13,866
$\overline{x}$ (τ2)	0.709	75,374	0.960	17,183	0.926	12,609
Δ(M1-0:1)	0.903	0.252	0.729	0.343	0.733	0.295
Ω(M1-0:1)	0.788	1.487	0.849	0.287	0.851	0.267
$\overline{x}$ (M1-0:1)	0.916	0.244	0.915	0.373	0.945	0.300

Table 4

Cross-validation of hierarchical co-demographic modeling under bird data specifications given no monomorphic sites nor normalization for hRF and hABC mean point estimates

	Weakly informed parameterization		Fixed ${\tau }_s^1$ parameterization		SC and fixed ${\tau }_s^1$ parameterization
	r	RMSE	r	RMSE	r	RMSE
	hRF hyperparameter estimation
ζ	0.536	0.244	0.673	0.214	0.736	0.199
	hABC hyperparameter estimation/model selection
ζ	0.435	0.259	0.631	0.236	0.668	0.227
ζM  1-0:1	0.551	0.275	0.335	0.268	0.576	0.240
ζM  1-1:0	0.661	0.238	0.674	0.222	0.786	0.182
	hABC parameter summary estimation
${\tau }_s^1$	0.494	725,994	N/A	N/A	N/A	N/A
Δ(τ1)	0.563	270,060	0.721	262,022	0.542	266,144
Ω(τ1)	0.518	313,422	0.750	284,508	0.493	281,559
$\overline{x}$ (τ1)	0.867	334,913	0.652	287,501	0.715	290,100
Δ(M1-0:1)	0.744	0.119	0.660	0.110	0.730	0.119
Ω(M1-0:1)	0.610	0.398	0.441	0.418	0.691	0.351
$\overline{x}$ (M1-0:1)	0.752	0.123	0.681	0.111	0.748	0.119
Δ(M1-1:0)	0.817	0.044	0.857	0.042	0.849	0.044
Ω(M1-1:0)	0.747	0.164	0.769	0.174	0.818	0.147
$\overline{x}$ (M1-1:0)	0.817	0.045	0.859	0.043	0.872	0.044

Our simulation experiments demonstrated nearly identical sensitivities among the two lamprey sampling projections, suggesting that sampling more individuals throughout this range yields negligible benefit, if any. In contrast, the analyses based on the bird data show a surprising decrease in hABC power when monomorphic sites were incorporated. Notably, there was greater error in including monomorphic bins, which we coupled with the normalization of ajSFS bins, than normalizing the ajSFS containing only SNPs, corroborating that monomorphic sites specifically impacted hABC performance. However, using monomorphic sites slightly improved ζ inference with hRF (r = 0.725 with monomorphic sites; r = 0.717 with normalization only; r = 0.673 without either), which inherently weighs the importance of the various bins through its machine learning algorithm. Together, this indicates that although there is additional signal to be gained from sequence length, exploiting this information is not straightforward (e.g., must optimize weights between monomorphic and SNP bins, which can be several orders of magnitude different in number) and may only marginally improve accuracy.

Finally, it is strongly evident from the hRF confusion matrices (Tables S6, S9, S12, S15, and S18) and mean of the hABC hyperposterior distributions per true value (Tables S7, S10, S13, S16, and S19) across all the simulation experiments, as well as the leave-one-out cross-validation plots (Figs. S1 and S2), that the signal for hyperparameter estimation was substantially stronger when the true value of ζ was 1.000 (i.e., completely synchronous co-divergence times), whereas there was a bit of upward bias when ζ approached idiosyncrasy.

EMPIRICAL ANALYSES OF REPLICATE LAMPREY SPECIES PAIR

There was a range of 7870–11,899 SNPs produced for the eight lamprey joint-SFS vectors when down-projected to 20 haploids, with a minimal decrease to 6263–9580 SNPs after further down-projection to 10 haploids (Table S20). Individual ABC demographic analyses unveiled a universal pattern across all eight population pairs of ancient M² estimates near zero in value (mode M^2-0:1 = 0.004–0.020; mode M^2-1:0 = 0.006–0.018), indicative of an SC scenario that followed a period of isolation (Tables 5 and S21) and hence our SC parameterization for the downstream hierarchical co-demographic model. Migration rates from the parasitic anadromous river populations of L. fluviatilis to nonparasitic freshwater brook populations of L. planeri remained minimal to present day, with mode estimates of less than 0.01 for five of the population pairs (mode M^1-1:0 = 0.000–0.017). Notably, although this set of migration results would support fixing M^1-1:0 to a value of zero in the SC parameterization as we did for the M² parameters, we chose to allow M^1-1:0 vary within a highly constrained range at lower values to remain more consistent with the original modeling from Rougemont et al. (2017) (Table S3). In contrast, migration rates in the opposite direction from the start of SC to present day were high, with M^1-0:1 mode estimates >1.0 for all but one population pair and on the order of several haploid migrants per generation for half of the population pairs (mode M^1-0:1 = 0.953–7.225).

Table 5

Lamprey individual population pair ABC mode point estimates

	Connected population pairs (CPP)				Disconnected population pairs (DPP)
	AA	BET	OIR	RIS	BRE	CEN	ODO	SAU
τ1	523,114	558,859	882,020	819,962	854,589	884,503	845,225	861,784
τ2	142,646	121,057	353,965	365,931	320,608	314,055	456,163	305,911
M  1-0:1	2.934	3.959	7.225	4.081	1.354	1.334	0.953	1.359
M  1-1:0	0.003	0.002	0.017	0.000	0.000	0.012	0.005	0.017
M  2-0:1	0.015	0.013	0.004	0.006	0.005	0.019	0.020	0.017
M  2-1:0	0.007	0.006	0.007	0.016	0.008	0.018	0.013	0.018
ε  0	0.093	0.187	0.299	0.359	0.307	0.493	0.307	0.420
ε  1	0.395	0.383	0.460	0.547	0.466	0.507	0.603	0.399
N  0	760,517	737,342	674,257	809,991	642,219	619,047	450,862	607,578
N  1	12,241	13,092	4572	12,159	5604	5596	3922	5712

Expansion magnitude estimates were low to moderate for L. fluviatilis populations (mode ε⁰ = 2.0 × to 10.8 ×) and consistently low for L. planeri populations (mode ε¹ = 1.6 × to 2.4 ×), with the latter result motivating us to model twofold growth in the L. planeri populations (i.e., population “1”) among all population pairs for the SC parameterization. The ranges for the two historical event timings of divergence and migration/size change were moderate (mode τ¹ = 523,114–884,503, with six of the eight pairs clumping ca. 850,000; mode τ² = 121,057–456,163, with five of the eight pairs clumping ca. 350,000), making the degree to which times were in synchrony ambiguous based on these individual population pair analyses alone. Notably, τ¹ and τ² were highly correlated across population pairs (r = 0.868 given mode estimates), although this could largely be an effect of the τ² prior being conditional on τ¹. Nonetheless, this suggests that regardless of the actual timing values, which coalescent theory dictates are inherently scaled by N_E and thus absolute time estimates can only be achieved by assuming a value for N_E or including information about monomorphic sites while incorporating a mutation rate parameter, there is a degree of congruence in the temporal length of isolation throughout lamprey population pairs.

Our hierarchical co-demographic approach using hRF and hABC (i.e., posterior distributions for ζ, ${\tau }_s^1$⁠, Δ(τ¹), Ω(τ¹), and $\overline{x}$ (τ¹); see Table S1 for descriptions of parameter symbols) supports a high degree of temporally synchronous co-divergence during the late Pleistocene (i.e., at τ¹) of the freshwater brook populations from their anadromous river complements (Figs. 4 and S3; Tables 6 and S22). Additionally, our hierarchical analyses elucidated key questions regarding variability in SC timing as well as strength of SC co-migration. In fact, the collective SC times (i.e., at τ²) among the eight lamprey population pairs were inferred to be similarly synchronous based on posterior distributions for ζ_τ2, Δ(τ²), and Ω(τ²), a perhaps unexpected insight into the co-migration dynamics between the nonparasitic L. fluviatilis and parasitic L. planeri populations. Notably, the hierarchical model constrained the start of SC to coincide with potential size change, which could have contributed to this co-demographic signal. For co-migration intensities, ζ_M1-0:1, Δ(M^1-0:1), and Ω(M^1-0:1) posterior distributions were shifted toward idiosyncrasy with respect to the corresponding priors, representative of incongruent contemporary migration rates from L. fluviatilis freshwater localities to L. planeri anadramous habitats. Point estimates of Δ(M^1-0:1) further imply wide variation of migration rates across population pairs, indicating a total range of at least four haploid migrants per generation, assuming the minimum total variability equals 2 × Δ. Moreover, our PCA-based predictive checks constructed against both simulations retained from the prior as well as re-simulations from the posterior, with subsequent projections of the empirical data, demonstrate that our hierarchical co-demographic model could reasonably generate the observed lamprey data (Fig. S4).

Figure 4

Lamprey hABC prior and posterior distributions. Prior histograms in light blue and posterior histograms in dark green of, from left to right and then top to bottom: ζ, Δ(τ¹), Ω(τ¹), ${\tau }_s^1$⁠, ${\zeta }_{{\tau }_2}$⁠, Δ(τ²), Ω(τ²), $\overline{x}$ (τ²), ζ_M1-0:1, Δ(M^1-0:1), Ω(M^1-0:1), and $\overline{x}$ (M^1-0:1). Distributions are shown for the SC parameterization as well as SC and fixed ${\tau }_s^1$ parameterization, given 10 haploids. For each distribution, the y-axis is in units of the probability mass, and the x-axis is in units of the corresponding hyperparameter or parameter summary. Vertical lines in green, blue, red, and dotted black represent mean, median, and mode point estimates as well as 50% credibility intervals, respectively.

Table 6

Lamprey hierarchical co-demographic hRF as well as hABC mean point estimates and 50% credibility intervals given 10 haploids

	SC parameterization			SC and fixed ${\tau }_s^1$ parameterization
	hRF hyperparameter estimation
ζ	0.841			0.894
	Mean point estimate	50% credibility intervals		Mean point estimate	50% credibility intervals
	hABC hyperparameter estimation/model selection
ζ	0.785	0.625	1.000	0.861	0.750	1.000
${\zeta }_{{\tau }_2}$	0.857	0.750	1.000	0.898	0.875	1.000
ζM  1-0:1	0.146	0.125	0.125	0.137	0.125	0.125
	hABC parameter summary estimation
${\tau }_s^1$	452,288	357,693	510,347	N/A	N/A	N/A
Δ(τ1)	85,625	0	150,088	62,568	0	104,898
Ω(τ1)	64,478	0	112,096	55,439	0	101,084
$\overline{x}$ (τ1)	469,010	417,493	516,582	419,315	361,117	460,631
Δ(τ2)	37,361	16,686	53,530	29,559	12,833	40,861
Ω(τ2)	23,026	3179	38,287	19,819	2244	34,052
$\overline{x}$ (τ2)	151,655	134,974	165,902	137,221	121,016	148,784
Δ(M1-0:1)	2.100	1.761	2.405	2.155	1.816	2.464
Ω(M1-0:1)	1.947	1.542	2.363	2.000	1.601	2.401
$\overline{x}$ (M1-0:1)	3.772	3.477	4.089	3.754	3.442	4.063

EMPIRICAL ANALYSES OF CO-DISTRIBUTED BIRD TAXON PAIRS

In the case of the six dry forest bird joint-SFS vectors, there was an imputed total of 260,637–849,951 sites given the extracted information about invariants as well as number of SNPs ranging from 5446 to 12,455 (Table S23). Unlike for the lampreys, individual-level ABC estimates for migration parameters were generally quite low, with 18 of the 24 mode M estimates less than 0.01, or one effective haploid migrant every 100 generations (maximum mode M = 1.533; Tables 7 and S24). In fact, both ancient phase M² migration parameters were approximately zero in value (i.e., <0.01) for all six bird population pairs, whereas both recent phase M¹ migration parameters had much more variance. Notably, although this supports an SC parameterization for the downstream hierarchical co-demographic model, to remain more consistent with the original modeling from Oswald et al. (2017), we chose to continue allowing the M² parameters vary within a highly constrained range at lower values until the final parameterization option, thus we first explored the effect of solely fixing ${\tau }_s^1$ (Table S4).

Table 7

Bird individual population pair ABC mode point estimates

	Arremon	Campylorhynchus	Melanopareia	Mimus	Saltator	Thamnophilus
τ1	553,978	1,743,372	1,210,602	264,790	627,849	479,811
τ2	188,695	869,525	700,885	135,352	378,118	172,539
M  1-0:1	0.005	1.369	0.003	1.533	0.018	0.448
M  1-1:0	0.013	0.003	0.001	0.005	0.003	0.446
M  2-0:1	0.006	0.002	0.005	0.004	0.006	0.001
M  2-1:0	0.000	0.001	0.001	0.000	0.001	0.000
ε  0	0.091	0.063	0.068	0.082	0.085	0.052
ε  1	0.072	0.065	0.071	0.062	0.091	0.234
N  0	444,477	406,780	357,884	453,885	460,533	273,187
N  1	826,113	652,565	940,162	556,931	615,656	601,824
U	1.324 × 10–9	1.609 × 10–9	3.869 × 10–9	1.317 × 10–9	1.322 × 10–9	1.131 × 10–9

In contrast to migration, inference of size change magnitude was consistently moderate to high (mode ε⁰ = 10.9 × to 19.2 × expansion; mode ε¹ = 4.3 × to 16.1 × expansion). For the two historical event timings, there was a much greater dispersion of estimates among population pairs (mode τ¹ = 264,790–1,743,372; median τ² = 135,352–869,525) in comparison to the lamprey inferences, thus providing a clearer test case whereby the expectation is asynchrony. Interestingly, however, τ¹ and τ² were still highly correlated (r = 0.976 given mode estimates); as with the lamprey analyses, this may be due to the more recent time being inherently conditional on the earlier time, yet implies that the temporal length of isolation is shared across the bird population pairs.

Corroborating the individual-based approach, hRF estimates and hABC posterior distributions for ζ, Δ(τ¹), and Ω(τ¹) (see Table S1 for descriptions of parameter symbols) suggest temporal idiosyncrasy in co-divergence deep within the past between the Marañón Valley and Tumbesian region (Fig. 5; Tables 8 and S25). Although these inferences allow some degree of intermediate synchrony, there is little room for interpretation of a fully shared, singular co-divergence event. Furthermore, we found similar asynchronous patterns for present-day migration rates in both directions, with the post hoc hyperparameter and parameter summary estimates consistent with a range of migration intensities from nearly nonexistent to moderate among the six taxon pairs. Together, these results from our hierarchical analyses indicate high variability in co-demographic dynamics, portraying a scenario of many independent population pair trajectories. Furthermore, our PCA-based prior and posterior predictive checks exhibit reasonable overall goodness-of-fit for our parameterization options. Notably, although the posterior predictive check for the fixed ${\tau }_s^1$ parameterization suggests a small degree of model mismatch, this could be due to an insufficient number of simulations for adequately sampling the hyperparameter space, which is implied by the satisfactory model fit for the more informed SC and fixed ${\tau }_s^1$ parameterization option that succeeded and is in fact nested within this fixed ${\tau }_s^1$ parameterization (Fig. S5).

Figure 5

Bird hABC prior and posterior distributions. Prior histograms in light blue and posterior histograms in dark green of, from left to right and then top to bottom: ζ, Δ(τ¹), Ω(τ¹), ${\tau }_s^1$⁠, ζ_M1-0:1, Δ(M^1-0:1), Ω(M^1-0:1), $\overline{x}$ (M^1-0:1), ζ_M1-1:0, Δ(M^1-1:0), Ω(M^1-1:0), and $\overline{x}$ (M^1-1:0). Distributions are shown for the fixed ${\tau }_s^1$ parameterization as well as SC and fixed ${\tau }_s^1$ parameterization. For each distribution, the y-axis is in units of the probability mass, and the x-axis is in units of the corresponding hyperparameter or parameter summary. Vertical lines in green, blue, red, and dotted black represent mean, median, and mode point estimates as well as 50% credibility intervals, respectively.

Table 8

Bird hierarchical co-demographic hRF as well as hABC mean point estimates and 50% credibility intervals

	Fixed ${\tau }_s^1$ parameterization			SC and fixed ${\tau }_s^1$ parameterization
	hRF hyperparameter estimation
ζ	0.535			0.547
	Mean point estimate	50% credibility intervals		Mean point estimate	50% credibility intervals
	hABC hyperparameter estimation/model selection
ζ	0.634	0.333	0.833	0.635	0.500	0.833
${\zeta }_{{\tau }_2}$	0.347	0.125	0.667	0.376	0.167	0.667
ζM  1-0:1	0.211	0.125	0.125	0.238	0.167	0.167
	hABC parameter summary estimation
Δ(τ1)	372,473	136,303	587,086	364,315	130,951	573,636
Ω(τ1)	379,101	95,619	651,958	371,343	87,854	601,943
$\overline{x}$ (τ1)	1,238,151	1,089,975	1,366,331	1,221,745	1,043,583	1,332,173
Δ(τ2)	0.241	0.126	0.333	0.235	0.117	0.335
Ω(τ2)	0.757	0.393	1.095	0.747	0.344	1.108
$\overline{x}$ (τ2)	0.249	0.131	0.342	0.242	0.121	0.341
Δ(M1-0:1)	0.131	0.087	0.169	0.138	0.094	0.174
Ω(M1-0:1)	0.285	0.164	0.384	0.282	0.170	0.376
$\overline{x}$ (M1-0:1)	0.142	0.097	0.181	0.151	0.103	0.191

Discussion

COUPLING HIERARCHICAL CO-DEMOGRAPHIC MODELING OF POPULATION PAIRS WITH THE ajSFS

Bayesian inference of demographic history within the context of population genetics and phylogeography, especially in the genomic era, is often accomplished via approximation through Monte Carlo simulations (Beaumont et al. 2002; Beaumont and Rannala 2004; Beaumont 2010; Hoban et al. 2012; Hoban 2014). This allows flexibility to tractably test complex demographic scenarios assuming a large quantity of data, such as with our hierarchical co-demographic model using the ajSFS. However, as a central tenet of Bayesian statistics, it is best practice to incorporate sensible priors to inform a statistical model (Rannala 2015), which may be derived through empirical Bayes and/or converged toward with sequential ABC. Although model accuracy is then conditional on the assumptions imposed from that prior information, a model that is instead too permissive can become problematic due to identifiability issues and expensive efforts searching probability space. This is an essential consideration for our co-demographic model, which is rich in nuisance parameters from the individual population pairs. Indeed, this concern was highlighted in our weakly informed parameterization options, which may have had parameter draw combinations across different hyperparameter values resulting in similar ajSFS vectors, and/or a large proportion of wasted simulations from improper sampling of parameter space. Fortunately, using more efficient prior distributions that enforced an SC scenario and/or fixed ${\tau }_s^1$ resulted in drastic improvements in accuracy while maintaining utility of the model as there was a high degree of flexibility still permitted and robustness to misspecification of fixed parameters. Importantly, we illustrated through our preliminary empirical analyses that our generalized population pair model was able to efficiently identify SC as a mutually applicable migration scenario. However, in cases whereby co-demography cannot be compressed into a simplified scenario such as SC (e.g., when migration trajectories are dynamic or the demographic syndrome varies among constituent population pairs), we also demonstrated that incorporating our generalized approach within a hierarchical model is a valid choice. In such instances, parameter space can still be reduced through limiting certain taxon-specific prior distributions (e.g., for ε parameters), and/or likelihood space can be more thoroughly explored through an increase in simulations.

With respect to potential sensitivity to sampling configuration, there was negligible change in accuracy as well as actual estimates between the alternative lamprey sampling projections. This is perhaps unsurprising given similar observations with the aSFS (Xue and Hickerson 2015) and multidimensional SFS (Smith et al. 2017). Regardless, this constitutes a significant practical benefit to potential users because ajSFS vector size scales very quickly relative to haploid sampling, which can incur intense computational cost especially with a large number of simulations. In other words, sampling at a lower number of individuals within this range ostensibly has little direct effect on the statistical power or actual results from this type of analysis, at least assuming the specifications of the lamprey dataset. Importantly, however, sampling more individuals and subsequently down-projecting could perhaps indirectly increase accuracy by producing higher fidelity empirical data vectors (i.e., better representative of the sampled population), due to integrating over subsamples from a more comprehensive sampling scheme as well as minimizing the impact of errant individual calls. Nonetheless, our in silico experiments show promise for scaling down sampling efforts and sequencing costs when designing a project of this nature. However, because this study may not universally apply to other potential datasets, researchers ought to run their own cross-validation analyses to precisely assess optimal sampling. Conversely, more attention should be devoted to sampling more population pairs, as implied when comparing leave-one-out RMSEs given the lamprey dataset with eight pairs versus the bird dataset with six pairs, and as previously demonstrated with related comparative methods (Chan et al. 2014; Xue and Hickerson 2015, 2017).

In contrast, our effort to employ sequence length information under the bird data specifications did not yield substantial improvement. Performance actually worsened for hABC inference, and although hRF estimation indeed benefitted from incorporating monomorphic bins, suggesting that some information can be extracted through optimized weighing of bins, the gain was slight. In fact, hRF accuracy increased much more from normalizing bins than using invariant sites. Hence, rather than incorporating the total set of sites to calibrate absolute estimates, there is perhaps greater value in using informed prior distributions, which can be derived from a separate analysis better suited for exploiting sequence data.

Our leave-one-out experiments demonstrated that useful co-demographic hyperparameters, specifically the post hoc conversion of taxon-specific parameter draws into ζ_T estimates, can be extracted even when not initially used to construct the hierarchical model. Although hyperprior distributions would be extremely skewed under such post hoc conversion as a consequence of random parameter draws among independent units without hyperparameter governance, this strategy provides convenient exploration of co-demographic results without additional expensive computation. Post hoc converted hyperposterior distributions may also be used to guide and inform subsequent hierarchical co-demographic modeling efforts (e.g., fixing τ² to the same value for all population pairs based on a result of ζ_τ2 = 1.000), akin to our empirical Bayes tactic of fixing ${\tau }_s^1$⁠. Additionally, there was evident power to infer the parameter summary Δ (Table S1) during the simulation exercises. Considering its straightforward co-demographic interpretation, that is, the average difference per population pair from the median value, and lack of assumption in distributions such as normality, this suggests that Δ estimation will provide great utility for future co-demographic analyses.

We introduced an approach that significantly simplifies modeling multiple population pairs by ignoring shared ancestry and assuming exchangeability, yet has minimal information loss with respect to co-demography. Specifically, our findings validate how the ajSFS can be exploited by a hierarchical co-demographic model of population pairs for direct estimation of congruence in parameter values. Importantly, we showcase a wide range of sampling specification, including number of individuals from 8 to 20 haploids, inclusion versus omission of monomorphic sites, and amount of SNPs. Moreover, the organismal structure of the two biological systems is of two extremes, whereby multiple replicates of a single species pair are dispersed throughout different localities and may possibly have a shared recent history in the lamprey case, whereas highly divergent population pairs are co-distributed across a shared barrier and have less a priori expectations of congruence for the birds. As a result, this development provides a broadly applicable method that can beneficially supplement analyses on individual datasets.

COMPARATIVE PHYLOGEOGRAPHIC HISTORY OF RIVER AND BROOK LAMPREYS IN WESTERN FRANCE

Our hierarchical co-demographic model elucidated a strong signal of synchronous co-divergence for the lamprey population pairs, which challenges the suggestion of idiosyncrasy by the neighbor-joining tree from Rougemont et al. (2017). This latter, less parsimonious scenario implies independent regional occurrences of colonization into different freshwater drainages at asynchronous times with complex inter-locality gene flow to maintain each of the two constituent species, also known as the transporter hypothesis (Rougemont et al. 2017). Instead, our results are consistent with either (1) a continental-scale event simultaneously yet separately splitting all of the population pairs, which were already spatially dispersed previously into eight ancestral populations, or (2) spatial expansion into the present total geographic distribution from a singular localized isolation of L. fluviatilis and L. planeri, which perhaps resided at one of the sites that presently harbors one of the population pairs. In the latter case, all the current-day species pair replicates coalesce into this same divergence. Importantly, our analyses do not directly discriminate these two possibilities, although the general correlation between genetic diversity and directionality, specifically an increase in number of SNPs toward more Southern latitudes with a stark jump at the halfway point, corroborates an intermediate scenario of Northern colonization from established population pairs spanning several of the Southern sampling sites (presumably OIR, CEN, ODO, and SAU). This is also consistent with the trajectories of many other organisms that persisted during the last glacial maximum (LGM) within refugia located in Southern Europe and subsequently expanded into the North following global warming ca. 20 KYA (Hewitt 2000). Our hABC approach also indicates a temporally shared collective start of SC, which likely happened before a putative post-LGM Northern expansion given the τ² and $\overline{x}$ (τ²) estimates. In Southern Europe, another incident of wide impact could have concurrently initiated collective SC across distinct species pair replicates that may have been established within refugia, as with the simultaneous co-divergence. However, the inference of highly variable co-migration magnitudes suggests that environmental differences existing among the sampling locations elicited a range of migration histories, signifying that Northbound spatial expansion coincided with a change in gene flow for colonizing lamprey population pairs. Notably, our results of synchrony within co-divergence and collective SC timings along with idiosyncratic co-migration support that variability in co-migration intensities rather than durations likely caused the heterogeneity in population pair genomic differentiation observed by Rougemont et al. (2017), addressing an open question that they posed.

The question of whether colonization or vicariance incited the species pair split may be resolved by our individual population pair estimates. Namely, the ubiquitous combination of considerable migration from anadromous L. fluviatilis to freshwater L. planeri with extremely low reverse-direction migration, as well as the order of magnitude larger N_E for L. fluviatilis populations over corresponding L. planeri, may point to one or more peripatric founder events as the likely beginning for the nonparasitic brook lamprey populations (Mayr 1982). On a related note, the generally increased strengths of expansion in the L. fluviatilis samples (≈2× to 11×) compared to the weak expansions of the L. planeri counterparts (≈2×) could also perhaps be further consistent with patterns left by allele surfing (Travis et al. 2007), such that spatial expansion originating from a geographically constricted distribution in the source resulted in an overall signal of demographic N_E expansion for the founding river lamprey, whereas the paired brook lamprey retained this at a lower intensity due to being relegated to the front of the wave.

In agreement with the population pair membership described in Rougemont et al. (2017), our individual population pair ABC inferences recovered appreciably stronger unidirectional migration rates in connected population pairs (CPP: AA, BET, OIR, and RIS; mode M^1-0:1 > 2.934 haploids/generation) versus disconnected population pairs (DPP: BRE, CEN, ODO, and SAU; mode M^1-1:0 < 1.360 haploids/generation). Notably, the CPP migration intensities potentially encroach on panmixia (Wright 1931; Whitlock and McCauley 1999; Beerli and Palczewski 2010), which is to be expected considering the associated environmental connectivity. However, in contrast to the initial lamprey study, we did not find support for either an AM or IM scenario across DPP. To reconcile this, it is possible that both approaches may be detecting different aspects of the same complex co-demographic trajectory. For example, a scenario of oscillation between isolation and contact could cause an SC, AM, or even IM inference depending on the temporal scale captured. In fact, Rougemont et al. (2017) speculated upon this very case for their AM model choice.

Altogether, our results demonstrate a shared biogeographic history that influenced congruence in co-demographic responses yet with specific characteristics, potentially including those related to environmental connectivity (i.e., CPP vs. DPP), still shaping migration differences throughout localities.

COMPARATIVE PHYLOGEOGRAPHIC HISTORY OF CO-DISTRIBUTED PASSERINES THROUGHOUT SEASONALLY DRY FORESTS OF THE MARAÑÓN VALLEY AND TUMBESIAN REGION IN PERU

Unlike our treatment of the lamprey data, we detected asynchronous patterns across all facets of bird population pair co-demography, with both our hierarchical endeavor as well as individual analyses indicating variability across co-divergence times, migration rates, and size change magnitudes. Consistent with the initial investigation (Oswald et al. 2017), we rejected a single mechanism causing co-divergence, and in fact our findings suggest more widespread idiosyncrasy than those authors first discovered. Specifically, our wide range of τ² estimates suggests highly variable temporal dynamics for postdivergence migration, addressing an open question posed by Oswald et al. (2017). Moreover, our various estimates for the M parameters demonstrate a spectrum of admixture syndromes. For example, we found convergence toward SI in Melanopareia versus a resemblance to SC in the other five bird taxa. Of those five population pairs, Campylorhynchus, Mimus, and Saltator displayed stronger unidirectional intensity from the Marañón Valley to the Tumbesian region, whereas Arremon had stronger gene flow in the opposite direction, and Thamnophilus exhibited a more symmetric balance. Opposed to the lamprey case, this variation in directionality would imply that neither direction of postdivergence migration was generally more favorable than the other, which was also posited in the initial study. Notably, our results indicate that the greatest exchange of migrants regardless of direction was in Mimus, Campylorhynchus, and Thamnophilus, with the latter two poorly fitting a two cluster model under a population structure method according to Oswald et al. (2017). Additionally, Oswald et al. (2017) stated that Thamnophilus was the only population pair in the dataset previously shown to be in SC within the Marañón Valley. Curiously, Thamnophilus also was estimated to have the strongest Marañón Valley expansion magnitude coupled with the weakest Tumbesian expansion magnitude across the bird pairs.

Such support for idiosyncrasy would then be consistent with independent colonizations occurring throughout time, with exact timing for each individual population pair perhaps dictated by species-specific dispersal ability (Smith et al. 2014) or differences in pre-existing competitors (Hardin 1960). Pleistocene climatic cycles may have also acted on these differences to mediate wide variation in co-divergence times as well as co-migration behavior. Likewise, less stable conditions within the Tumbesian region due to El Niño events may have the consequence of such observations too, as suggested by Oswald et al. (2017). Regardless, we are in strong agreement with the original investigation that vicariant allopatry from the Andean uplift forming a barrier within a wider distribution is a very unlikely explanation due to our timing of isolation being too recent and asynchronous.

Our hierarchical approach concurring with the original findings of idiosyncratic co-demography validates that our model experiences minimal information loss from its simplification of population pairs as exchangeable units, which discards information about identity. Moreover, our generalized population pair model was able to extend the conclusions of Oswald et al. (2017), discovering even greater independence among evolutionary trajectories in the historical assembly of these co-distributed passerine taxon pairs.

Conclusion

Investigating co-demographic histories among population pairs is a central aim of comparative phylogeography that focuses on estimating the geographic distributions and temporal variability of multiple sister species divergences (Taberlet et al. 1998; Avise 2000; Arbogast and Kenagy 2001; Hickerson et al. 2010; Barratt et al. 2018; Peñalba et al. 2018; Rincon-Sandoval et al. 2019). Further comparative insight could be attained, although by exploring congruence in other demographic parameters such as for migration and size change, as these might reflect other important trait-based factors underlying multi-taxa responses to biogeographic and climatic events (Papadopoulou and Knowles 2016). For example, similar migration and expansion behaviors across populations coupled with synchronous co-divergence times could reflect pervasive effects from a common vicariant event with little contribution from species-specific traits, whereas asynchronous migration or expansion trajectories mixed with synchronous co-divergence may suggest that species-specific traits effected differential postdivergence responses from a shared allopatric event. Such inquiries can elucidate other contexts beyond co-vicariance (i.e., allopatry), including heterogeneous environmental gradients that are perhaps facilitated by climate (i.e., parapatry) (Riddle et al. 2000; Hewitt 2001; Moritz et al. 2009), founder/colonization events (i.e., peripatry) (Hickerson and Meyer 2008; Leache et al. 2009), and widespread, commonly experienced biotic interactions such as symbiosis or predation (Stone et al. 2012; Satler and Carstens 2016; Garrick et al. 2017). A classical strategy for addressing questions of this nature is to conduct individual inferential analyses on the constituent taxa within a comparative dataset, the results of which are concatenated to heuristically determine congruence. Although straightforward, this approach deploys an undefined and therefore nonreproducible protocol for interpreting conclusions as well as ignores the uncertainty underlying each single estimate.

To this end, we presented here a hierarchical model that assesses several aspects of co-demography among independent population pairs, including synchrony in divergence and SC timings as well as congruence in migration rates. Our method is targeted for increasingly available reduced representation SNP data (e.g., RAD-seq) through its application of the ajSFS, complementing other recent advances that can likewise treat genomic-scale data from multiple taxon pairs, yet within a multi-level likelihood framework (Bunnefeld et al. 2018; Oaks 2018). However, our approach is differentiated by its flexible co-divergence model, which is designed to easily scale in number of population pairs and can accommodate complex population pair demography. Furthermore, our unique modeling scheme joins population pairs within a unified analysis without assuming shared ancestry among them, and detects shared patterns in focal parameter values while integrating out unknown population-specific nuisance parameters, for example, differences in effective population size. This allows ecological assembly-level questions to be explicitly addressed, measurements of uncertainty that consider stochastic and inter pair variance, and increased accuracy from pooling datasets with shared parameterizations. This methodology could be applied to many classic systems and questions, with future empirical implementations that could include regions such as Isthmus of Panama, Andes uplift, and Hawaiian islands chain, as well as organisms such as African cichlids and Heliconius butterflies. With the increase of studies publishing genomic-scale data for nonmodel population pairs on a comparative scale (Roux et al. 2016; Harvey et al. 2017; Rougeux et al. 2017; Stange et al. 2017; Barratt et al. 2018; Peñalba et al. 2018; Rincon-Sandoval et al. 2019), we anticipate this to be a useful and relevant addition to the suite of tools used within population genetics, biogeography, evolutionary biology, and community ecology.

AUTHOR CONTRIBUTIONS

ATX performed the research. ATX and MJH designed the study and wrote the manuscript.

ACKNOWLEDGMENTS

We would like to thank E. Ma for assistance in designing the model diagram figures; E. A. Myers and I. Prates for access to computational resources; A. D. Kern for thoughtful comments on the manuscript and lab support; A. C. Carnaval, S. Boissinot, and B. M. Henn for thoughtful comments on the manuscript; and A. Siepel for lab support. This work was supported by grants from National Science Foundation (DEB-1253710 to MJH; DEB-1343578 to A. C. Carnaval, MJH, and K. C. McDonald; DEB-1457232 to M. K. Fujita), National Institutes of Health (1R15GM096267-01 to S. Boissinot and MJH), FAPESP (BIOTA, 2013/50297-0 to A. C. Carnaval, MJH, and K. C. McDonald), and NASA through the Dimensions of Biodiversity Program. This work would not have been possible without help from the City University of New York High Performance Computing Center, with support from the National Science Foundation (CNS-0855217 and CNS-0958379).

DATA ARCHIVING

Scripts and result files have been deposited into Dryad (https://doi.org/10.5061/dryad.xksn02vbw).

Associate Editor: J. Light

Handling Editor: M. Servedio