multi‐dice: r package for comparative population genomic inference under hierarchical co‐demographic models of independent single‐population size changes

Abstract Population genetic data from multiple taxa can address comparative phylogeographic questions about community‐scale response to environmental shifts, and a useful strategy to this end is to employ hierarchical co‐demographic models that directly test multi‐taxa hypotheses within a single, unified analysis. This approach has been applied to classical phylogeographic data sets such as mitochondrial barcodes as well as reduced‐genome polymorphism data sets that can yield 10,000s of SNPs, produced by emergent technologies such as RAD‐seq and GBS. A strategy for the latter had been accomplished by adapting the site frequency spectrum to a novel summarization of population genomic data across multiple taxa called the aggregate site frequency spectrum (aSFS), which potentially can be deployed under various inferential frameworks including approximate Bayesian computation, random forest and composite likelihood optimization. Here, we introduce the r package multi‐dice, a wrapper program that exploits existing simulation software for flexible execution of hierarchical model‐based inference using the aSFS, which is derived from reduced genome data, as well as mitochondrial data. We validate several novel software features such as applying alternative inferential frameworks, enforcing a minimal threshold of time surrounding co‐demographic pulses and specifying flexible hyperprior distributions. In sum, multi‐dice provides comparative analysis within the familiar R environment while allowing a high degree of user customization, and will thus serve as a tool for comparative phylogeography and population genomics.

Originally developed for single-locus DNA data sets easily collected from multiple taxa (Burbrink et al., 2016;Hickerson, Stahl, & Takebayashi, 2007;Ornelas et al., 2013), this methodology has been extended to accommodate SNP data sets derived via recently emerging technologies such as RAD-seq and GBS, thereby improving inferential resolution through vastly greater sampling of independent gene tree histories across genomes from multiple taxa (Xue & Hickerson, 2015). This has been accomplished by exploiting the aggregate site frequency spectrum (aSFS), which has been established to contain signal of variability in demographic histories across taxa. Producing an aSFS involves creating single-population site frequency spectra (SFS) independently across taxa and combining these according to a standardized re-ordering procedure based on relative proportions of total SNPs per allele frequency class. This protocol therefore does not require sites to be homologous between taxa and in turn allows data to be collected across distantly-related taxa (more details about data preparation given in Implementation section). Construction of the aSFS can then be applied to coalescent simulations produced under a hierarchical co-demographic model that treats taxa as independent, unidentified and exchangeable units, and coupled with a statistical framework such as approximate Bayesian computation (ABC) to make comparative multi-taxa inference (Prates et al., 2016;Xue & Hickerson, 2015). This simulation approach could potentially be modified with other techniques, including machine learning algorithms such as random forest (RF) (D ıaz-Uriarte & Alvarez de Andr es, 2006; Pudlo et al., 2016;Strobl, Boulesteix, Zeileis, & Hothorn, 2007;Svetnik et al., 2003) and partial least squares regression (PLS) (Boulesteix & Strimmer, 2007;Wegmann, Leuenberger, & Excoffier, 2009). To elaborate, RF involves constructing decision trees based on "training" simulations to form a classification or regression scheme that subsequently can be applied to observed data, and PLS entails maximizing the variance explained in response variables in a manner similar to principal component analysis, which can be employed as a transformation procedure to potentially mediate high dimensionality of correlated summary statistics, such that inherently exists among aSFS bins. Alternatively, the aSFS could be deployed within a composite likelihood optimization (CL) framework, a statistical approach commonly used for demographic inference based on SFS data (Bustamante, Wakeley, Sawyer, & Hartl, 2001;Excoffier, Dupanloup, Huerta-S anchez, Sousa, & Foll, 2013;Gutenkunst, Hernandez, Williamson, & Bustamante, 2009;Lukic & Hey, 2012;Sawyer & Hartl, 1992).
The aSFS enables researchers to exploit data produced by nextgeneration sequencing to explore a variety of hypotheses that relate climatic and landscape changes with the evolution and demographic histories of biotic assemblages through hierarchical co-demographic modelling. Here we make this analytical pipeline available as the R package MULTI-DICE (Multiple Taxa Demographic Inference of Congruency in Events). To demonstrate and explore implementation of MUL-TI-DICE, we conducted a series of simulation studies that summarize an expanded set of options within our aSFS approach, including: (i) employing RF as an additional inferential tool; (ii) enforcing a "buffer" on prior space such that co-demographic events have an a priori minimal difference in time from each other; (iii) truncating the hyperprior range for improved hyperparameter estimation.

| Hierarchical co-demographic model
Our hierarchical co-demographic model consists of n taxa, which refer to independent panmictic populations with no assumption of or requirement for recent shared ancestry (Mazet, Rodr ıguez, Grusea, Boitard, & Chikhi, 2016), randomly assigned to Ψ instantaneous expansion ( Figure 1a) or contraction (Figure 1b) times. Of the Ψ times, there are w times corresponding to synchronous pulse events that involve at least two taxa, and r times corresponding to idiosyncratic events ungrouped from any pulses with only a single taxon, such that Ψ = w + r (Table 1). The proportion of n taxa assigned to any of the w pulses is represented by f T , the proportion of n taxa belonging to each of the w pulses is described by the associated hyperparameter vector f s = {f 1 , . . ., f w }, and the proportions of n taxa across all Ψ events are indexed by the vector f = {f s , f i,1 , . . ., f i,r }.
Here, f T is a single proportion value that ranges from 0.0 to 1.0 and equals the total sum of f s (i.e., f T = P w j¼1 f j when w > 0), and both f s and f are hyperparameter vectors that index proportion values across events. Specifically, each of w elements within the vector f s ranges from 2/n to 1.0, and f comprises of f s as well as each f i element = 1/n. The proportion f T and proportions within the vector XUE AND HICKERSON | e213 f s may be converted to numbers of taxa S T = f T 9 n and S = f s 9 n, respectively. Synchronous pulse times are indexed in the vector s s = {s s,1 , . . ., s s,w }, whereas idiosyncratic times are indexed in the vector s i = {s i,1 , . . ., s i,r }, with both vectors arranged in ascending order from most recent to oldest. To clarify, synchronous pulses are indexed by the temporal order established by s s = {s s,1 , . . ., s s,w }, which thus determines the order of f s such that f 1 pertains to the most recent pulse and f w reflects the most ancient. In the case of Ψ = w and r = 0, accordingly f T = 1.0 such that all taxa are assigned to one of w synchronous pulses with no temporally idiosyncratic taxa. On the other extreme, when Ψ = r = n and w = 0, accordingly f T = 0.0 with zero elements in the associated f s vector such that there are no synchronous pulses with all taxa idiosyncratically experiencing population size change across r different times. Other taxon-specific demographic parameters include each taxon's ratio of size change from the ancestral effective population size to current effective population size is indexed by the vector e = {e 1 , . . ., e n } and each taxon's current effective population size indexed by the vector N = {N 1 , . . ., N n }. Additionally, population size change times may be indexed to coincide with the taxa arrangement of e and N such that s = {s 1 , . . ., s n } (Table 1).
When implementing this co-demographic model for comparative demographic inference, there exists flexibility in the hierarchical parameterization, with several options available in MULTI-DICE. One such option, similar to the approach described in Chan, Schanzenbach, and Hickerson (2014) and Xue and Hickerson (2015), is to constrain the hyperparameter w to the values within the set {0, 1} and condition Ψ and r on the hyperparameter f T , which freely varies according to the hyperprior distribution P(f T ). This allows scenarios of complete idiosyncrasy, absolute synchrony within a single pulse, and intermediate degrees of synchronicity belonging to one pulse with remaining taxa temporally idiosyncratic. Here, f 1 is the only element possible in f s whereby f T = f 1 when w = 1 and f T = 0.0 when w = 0, resulting in the joint posterior distribution P (f T , s, e, N | Data) / P(Data | f T , s, e, N) P(e, N) P(s | Ψ, r, f T ) P(Ψ, r | f T ) P(f T | w < 1). The values for Ψ and r are then determined by Ψ = 1 + n À S T (when w = 1) and r = n À S T , respectively.
An alternative scheme is to randomly assign the proportions of n taxa to Ψ times according to the hyperprior distribution for the vector f, which is conditional on the hyperprior distribution of Ψ, with w and r accordingly conditional on P(f | Ψ) and P(Ψ). This leads to the joint posterior distribution P(Ψ, f, s, e, N | Data) / P(Data | Ψ, f, s, e, N) P(e, N) P(s | Ψ, f, w, r) P(w, r | Ψ, f) P(f | Ψ) P(Ψ). The values for w and r are then determined by the number of Ψ draws for f that are above and equal to 1/n, respectively, yielding the so-called Chinese restaurant process (Aldous, 1985;Blei, Griffiths, Jordan, & Tenenbaum, 2003) that is similarly applied in msBayes (Hickerson et al., 2007;Huang, Takebayashi, Qi, & Hickerson, 2011). Similarly, a third scheme is to condition the hyperprior distribution for the vector f s , which must have a lower bound greater than 1/n, on the hyperprior distribution of w, with Ψ and r accordingly conditional on P(f s | w) and P(w), such that the joint posterior distribution is P(w, f s , s, e, N | Data) / P(Data | w, f s , s, e, N) P(e, N) P(s | w, f s , Ψ, r) P(Ψ, r | w, f s ) P (f s | w) P(w). The values for Ψ and r are then determined by Ψ = w + n À S T and r = n À S T , respectively. Optionally, for each possible value in the w hyperprior, the associated f s , Ψ and r values may be fixed to specified values rather than allowed to vary.

| Simulation experiments
We conducted a series of in silico experiments to quantify accuracy and bias for various inferential frameworks and hierarchical codemographic modelling variants. Data were simulated under known hyperparameter and parameter values with the coalescent simulator FASTSIMCOAL version 2.5 (Excoffier et al., 2013). To directly generate single-population folded SFS, the FREQ setting was enabled assuming a set number of independent genealogies per SFS, which was treated as an approximation for the number of SNPs sampled and differed between experiments. Each SFS contained 20 haploid samples, only polymorphic bins and proportional SNP frequencies rather than total SNP counts. Per individual simulation, a set of 10 SFS corresponding to n = 10 populations was converted to a single aSFS summary vector following Xue and Hickerson (2015). Simulation reference tables composed of hyperparameter and parameter values randomly drawn from their respective hyperprior and prior distributions and their corresponding aSFS summaries were separately ε 1 ε 9 ε 8 ε 7 ε 6 ε 5 ε 4 ε 3 ε 2 ε 10 ε 1 ε 9 ε 8 ε 7 ε 6 ε 5 ε 4 ε 3 ε 2 ε 10 under the simple rejection algorithm against pseudo-observed data sets (PODs). PODs were produced under one of two methods, either independently from the reference table or using the "leave-one-out" cross-validation procedure. In brief, the "leave-one-out" procedure involves iteratively treating a single randomly selected simulation from a reference

| Testing inferential frameworks
In addition to hRF and hABC, we coupled these frameworks with For PLS, the plsr() function in the R package PLS (Mevik & Wehrens, 2007) was applied to a random subset of 10,000 simulations against variation in w per reference Vector of proportions of taxa belonging to each event, thus including f s , ordered such: {f s , f i,1 , . . ., f i,r }, with each f i element = 1/n; hyperparameter that directly governs s f s Vector of proportions of taxa belonging to each pulse {f 1 , . . ., f w }, ordered from most recent to most ancient; hyperparameter that directly governs s s ; each element ranges from 2/n to 1.0 f i An element of f or f s as the index j iterates from 1 to Ψ or w, respectively S T Conversion of f T to numbers of taxa by f T * n; n = S T + r hCL, a custom pipeline that calls dadi to calculate the expected SFS (Gutenkunst et al., 2009) and incorporates the multinomially distributed CL equation utilized in FASTSIMCOAL2 (Excoffier et al., 2013) and the BFGS optimization algorithm (Liu & Nocedal, 1989) was implemented (Supporting Information).

| Pulse buffer on prior space
Estimation of w or Ψ can be problematic as it does not necessarily correlate well with true temporal variability in co-demographic events. For example, a large number of synchronous events closely clustered in time would signify a high w value yet have low temporal variability, whereas a history with two synchronous co-demographic events that are far apart in time would yield a lower w value (w = 2) but with higher variance in time. As is the case with previous implementations of hierarchical co-demographic models , this inconsistency can hinder the ability to capture meaningful signal of w contained within the aSFS. To improve w estimation, we deployed a user-defined temporal pulse buffer that defines a minimal threshold of time b surrounding each co-demographic event such that for each jth event, all other co-demographic events occur outside a s j AE b window. Mechanistically, this involves sequentially modifying the prior distribution with every subsequent s draw, with final assignment of {s s,1 , . . ., s s,w } in ascending order such that s s,1 is the most recent and s s,w is the most ancient. For example, given a simulation with values w = 2, s~U{10,000, 1,000,000} and b = 20,000, if the first s s draw is 100,000 generations, then the second s s draw would be from the set U{10,000, 79,999} ∪ U {120,001, 1,000,000}; and if the second s s draw is 15,000 generations, then {s s,1 , s s,2 } is assigned such that s s,1 = 15,000 and s s,2 = 100,000. Importantly, a limit on the allowable number of buffered co-demographic events is imposed by the total s prior distribution across these events and the magnitude of b.

| Testing pulse buffer on prior space
To gauge how b impacts hyperparameter estimation, two reference tables with b = 0 generations and b = 30,000 generations were generated. In the special case of w = 0 for the b = 30,000 reference table, b was reduced to 10,000 due to the constraint from the s prior range and to allow more flexibility in the temporal dispersion for the total idiosyncrasy scenario. Both reference tables contained 100,000 aSFS simulations of instantaneous co-expansion ( Figure 1a) per value of w~U{0, 5} for a total of 600,000 simulations each. For simplicity, idiosyncratic taxa were not permitted and f T = 1.0 was evenly distributed across the vector f s for each value of w > 0 (Table 2). Importantly, to accommodate the special case of w = 0, which is equivalent to Ψ = 10, whereas all other values of w result in Ψ = w, w values were converted to Ψ for estimation purposes.
Single-population SFS were generated from 5,000 independent genealogies and according to the prior distributions s~U{5,000, 250,000} (in units of number of generations), e~U(0.01, 0.10) and N~U{50,000, 250,000}.
The "leave-one-out" cross-validation procedure was performed on each reference

| Testing inferential frameworks
The inferential frameworks that demonstrated the highest accuracy and precision in estimating w were hRF (r = .600-.807, comparison with the "leave-one-out" cross-validation on the b = 0 reference table. Additionally, buffering appears to benefit hyperparameter estimation without substantially affecting hABC estimation of parameter summaries Ω s and E(s). Notably, hRF again outperformed hABC in Ψ estimation, although this was minimal.
Better performance in Ψ estimation is apparent when truncated hyperprior ranges were employed, with w~U{0, 3} possibly the best compromise here between a more flexible hyperprior and greater accuracy ( Lopes & Beaumont, 2010), especially when considering efficiency with respect to a finite sampling of parameter space .
In the case here of hierarchical co-demographic models, rather than using Ψ or w, as well as other parameters such as s, in an arbitrary manner to merely construct the model, it can instead be specified meaningfully to gain insight about the variability in demographic changes across taxa given the temporal scale of interest.   (Bertorelle et al., 2010) and assessing goodness of fit with techniques such as prior and posterior predictive checks (Gelman et al., 2003;Lemaire, Jay, Lee, Csill ery, & Blum, 2016 (Ewing & Hermisson, 2010;Kern & Schrider, 2016), or analytical calculations of the SFS (Kamm, Terhorst, & Song, 2017;Wakeley & Hey, 1997). Notably, although our focus here is on the aSFS and accordingly reduced representation data sets (e.g., SNPs, RAD-seq, GBS), we acknowledge the great value in utilizing widely available mitochondrial/barcode-type data (Burbrink et al., 2016)  In its simplest operation then, MULTI-DICE can construct a reference

| Workflow
The function build.dice() is deployed first to construct hyperpriors across discrete hyperparameter values (i.e., Ψ, w, f T , f and f s ), allowing the following distributions: (i) a discrete uniform hyperprior on Ψ or w, depending on how the associated f vector is specified, then for f T within each discrete Ψ or w value, and finally across all combinations of the vector f or f s , respectively, within each discrete f T value; (ii) a Dirichlet-process hyperprior (Oaks, 2014)  play.dice() use the sample() function for random draws, each value in a user-specified distribution is treated as unique even when values are repeated (e.g., w 2 {0, 0, 0, 1, 2}), thus any discrete distribution (e.g., ln, gamma, beta) may be deployed for hyperpriors and priors.
Together, build.dice(), roll.dice() and play.dice() specify the hierarchical co-demographic model, as well as administer hyperparameter, parameter summary and taxon-specific parameter draws given this model.
Notably, data partitioning may be performed here (Prates et al., 2016), which allows heterogeneous specification of demographic scenarios (e.g., expansion, contraction), prior distributions, and data content and format (e.g., sampling of individuals, sampling time, polarization) across taxa within a data set; for example, data partitioning can accommodate a co-demographic model of expanders mixed with contractors at a pre-determined ratio.
In succession is dice.sims(), where FASTSIMCOAL2 is called to simulate data independently per taxon. Here, heterogeneous generation times across taxa may be specified (Xue & Hickerson, 2015). Importantly, for genomic-scale data, either the FREQ setting may be activated to directly generate SFS, or the SNP setting may be employed, which allows the option of using a mutation rate prior and thus monomorphic sites; for single-locus data, the SNP setting is deployed. Simulated summary statistic vectors and associated hyperparameter draws, taxon-specific parameter values and optional parameter summaries are outputted to a user-specified directory as simple text files. The total number of outputted files equals the number of simulated taxa plus one file per hyperparameter, taxon-specific parameter vector and parameter summary chosen for output. As aforementioned, all the functions described thus far can be implemented together automatically within dice.sims(), although independently calling functions may afford enhanced customization.
Following dice.sims() is either dice.aSFS() or dice.sumstats(), depending on the data scale (i.e., genomic or single locus, respectively). For XUE AND HICKERSON | e219 dice.aSFS(), the independent taxon-specific SFS are rearranged into a single aSFS according to the procedure outlined in Xue and Hickerson (2015), and for dice.sumstats(), the first four moments (i.e., mean, variance, skewness and kurtosis) are calculated for each of the four summary statistics (i.e., number of haplotypes, haplotype diversity, nucleotide diversity and Tajima's D) of the single-locus sequence block across the multiple taxa, for a total of 16 multi-taxa summary statistics, following Chan et al. (2014). For both of these functions, the user specifies the directory containing the simulation files, with simple specification for multiple directories resulting from parallelized runs, and the subsequent conversion is outputted within R, enabling easy piping into an inferential package such as ABC and/or writing to a simple text file. Importantly, these two functions can be applied to convert empirical data as well. Additionally, neither function calls upon any other MULTI-DICE functions and thus must be used in conjunction with at least dice.sims().
Advantageously, data type is irrelevant in all functions until dice.sims(), for which the data type is easily specified in a single argument and there is no disparity in output format. Hence, hierarchical codemographic models can be specified with the same level of complexity and flexibility for single-locus data as genomic-scale data in MULTI-DICE. Furthermore, dice.aSFS() and dice.sumstats() operate analogously and have near identical arguments, resulting in equivalent procedures for both data types with negligible difference. This feature lends itself nicely to conveniently analysing both data types for the same system either consecutively or simultaneously.

| Data sampling and processing
Although not directly handled by the MULTI-DICE package, we discuss here our recommendations for the practice of obtaining and preparing data. We emphasize that our methodology assumes dice.sumstats() F I G U R E 2 Flowchart of MULTI-DICE usage. MULTI-DICE accomplishes multi-taxa co-demographic inference under a hierarchical model through three major steps: model specification, single-population simulation across multiple taxa and conversion of simulated data to multi-taxa summary statistics. Hierarchical co-demographic model specification is conducted across multiple functions in sequence, with preceding functions contained within successive functions. This sequential embedding of functions extends to dice.sims(), allowing the entire model specification process to be performed concurrently with data simulation. Simulated data can then be converted to multi-taxa summary statistics by either dice.aSFS() or dice.sumstats(), depending on the data type. Additionally, these functions can be applied to empirical data as well. To clarify, only two MULTI-DICE functions/command lines, dice.sims() and dice.aSFS()/dice.sumstats(), are needed for simplest usage to construct a reference table of multi-taxa summary statistics under a hierarchical co-demographic model. This reference table can then be exploited in a downstream software program for hRF or hABC purposes, where appropriate statistical practices should be used to examine robustness and fit. Importantly, exploratory analyses should be performed on the empirical data prior to deploying MULTI-DICE to better guide its usage, for example, to determine sensible prior distributions and evaluate differences among taxa population-level sampling of multiple independent taxa, thus necessitating a sufficient number of samples per panmictic population (Robinson, Coffman, Hickerson, & Gutenkunst, 2014), which would depend on the temporal scale under investigation (Keinan & Clark, 2012). Importantly though, there is greater statistical resolution gained with increasing numbers of taxa (Chan et al., 2014;Xue & Hickerson, 2015), such that more emphasis should be placed on producing data sets with greater taxa representation rather than population-level sampling. To achieve this, investigators can benefit from splitting species/complexes into multiple independent structured populations that are determined from a preliminary exploratory analysis (Frichot, Mathieu, Trouillon, Bouchard, & Franc ßois, 2014;Patterson, Price, & Reich, 2006). This is especially important as lumping samples from multiple subdivided populations can result in strong bias when estimating population size changes (Mazet et al., 2016). While splitting indeed neglects shared ancestry, this problem may be negligible if isolation times are older than the co-demographic events of interest. Relatedly, conducting a cross-validation analysis across various sampling schemes, including both number of samples per taxon and number of taxa, prior to data collection and sequencing can be particularly informative of the proper sampling required for a given study (Bertorelle et al., 2010).
Greater statistical strength is gained with increasing taxa membership, but a strategy of indiscriminately adding taxa without consideration of specific characteristics can restrict researchers to testing generic hypotheses about assemblage-level demographic responses to shared conditions (Papadopoulou & Knowles, 2016). In consideration of this, we encourage researchers to, whenever possible, delineate data sets based on guilds that share a trait of interest.
We highlight here that the aSFS is capturing information within multiple independent structured populations, particularly size change history, through an aggregation of independent single-population SFS vectors. This operates somewhat differently than a joint-SFS or multi-SFS across multiple related populations, which also contains information about divergence and migration from shared and fixed polymorphisms (Wakeley & Hey, 1997). By focusing on solely within-population polymorphisms and being exploited to test hypotheses about size change history across taxa that may have experienced shared responses to climatic and habitat change while ignoring inter-population relationships, the aSFS-based approach simplifies the modelling, eliminates certain assumptions (e.g., topology, nature and duration of migration) and allows the option to directly test hypotheses across co-distributed taxa. On a related note, if SNPs are pruned to one per locus to avoid linkage disequilibrium violations prior to constructing the observed SFS, and if SNP calls were conducted across populations, then fixed polymorphisms should be removed before pruning to maximize the total number of SNPs per population.
Although the focus here on the aSFS has been exclusively regarding SNPs, MULTI-DICE is capable of incorporating monomorphic sites and accordingly mutation rates. Importantly, considering how s scales with N E in a coalescent model, if prior distributions exceed one order of magnitude for both parameters, then nonidentifiable SFS at different parameter combinations may be produced by ignoring monomorphic sites, thus potentially inflating bias and inaccuracy.
Hence, models that cannot have priors informed at least to this level may need to incorporate monomorphic sites. Assuming SNPs are pruned to one per locus, the number of monomorphic sites may be re-scaled given its ratio to the total number of SNPs. A prior for mutation rates must then be applied as well, which may result in this same identifiability issue if it likewise exceeds one order of magnitude. For this reason, users are advised to calculate population genetic summary statistics beforehand to assess the risk of incorporating taxa that vary to such extreme degrees as to falsely signal synchrony (Figure 2), which may be exacerbated with extremely phylogenetically distant taxa. For example, if the range in ratio of monomorphic to polymorphic sites among a multi-taxa data set greatly exceeds one order of magnitude, then extra considerations may need to be taken.

| Informing hierarchical co-demographic model
When conducting a multi-taxa co-demographic analysis using MULTI-DICE, the user is expected to assume a priori the composition of the demographic scenarios within the data set with respect to number of expanders and contractors, as well as accompanying prior distributions ( Figure 2). Furthermore, the aSFS requires that all singlepopulation SFS are at the same sampling level of individuals. This can be easily accomplished with the program dadi (Gutenkunst et al., 2009), but considering that multi-taxa data sets usually do not consist of a uniform sampling level, an optimal sampling projection must be selected. This optimal sampling projection is typically not readily apparent as the number of SNPs varies at different projection levels, with more SNPs discarded at higher sampling projections due to missing data and decreased singleton resolution at lower sampling projections resulting in low-frequency SNPs being assigned as monomorphic. Hence, to determine the optimal sampling projection across all taxa given this interplay between sampling of individuals and SNPs, as well as infer demographic scenarios with reasonable priors, an initial model-based investigation can be performed for each single-population taxon separately.
While this may be performed with CL-based methods such as dadi (Gutenkunst et al., 2009) or FASTSIMCOAL2 (Excoffier et al., 2013), an exploratory analysis across many independent taxa can be more efficiently conducted with an ABC approach, which allows quick inference for multiple empirical data sets against a single reference table and provides posterior distributions simultaneously with point estimates. MULTI-DICE coupled with an ABC framework then is well suited for efficiently performing a high throughput of such singlepopulation analyses to test models of demographic scenarios, explore various prior distributions and employ several data sampling XUE AND HICKERSON | e221 levels/projections. Notably, such a preliminary analysis may also be informative for multi-population demographic models, as well as elucidating results of synchrony from a co-demographic analysis by identifying candidate taxa potentially involved with synchronous pulses.

| CONCLUSION
The MULTI-DICE software package is designed for comparative population genetics and phylogeography and offers flexibility in user specification of hierarchical co-demographic models within a command-line interface R environment, a popular scripting language for population genetics (Paradis et al., 2017). This includes operating at different hierarchical levels (i.e., Ψ/w and f T /f/f s ), applying various demographic trajectories (including co-expansion and cocontraction) and implementing buffering on parameter values in prior space (b), for either genomic-scale or single-locus sequence data. Furthermore, there are several other features not discussed here that are available in MULTI-DICE, such as partitioning taxa into different modelling and data specifications within a combined analysis (Prates et al., 2016). Additionally, there exist options that offer greater flexibility within the co-demographic modelling, including incorporating two-event/three-epoch size change models, employing exponential rather than instantaneous growth and detecting congruence in other demographic parameters. This flexibility extends to data content and format as well, as MULTI-DICE also allows exploiting ancient samples, incorporating generation time heterogeneity, using polarized data (i.e., unfolded SFS), removing/adding allele frequency classes (e.g., avoiding classes more prone to error such as singletons, or including monomorphic sites and thus mutation rates and whole-locus information), and operating simulations under FASTSIMCOAL2's SNP model instead of its FREQ setting. Moreover, prior distributions can be highly customized, for example assigning different prior distributions between taxa within a shared pulse and those that are idiosyncratic, allocating alternative prior distributions per shared pulse and conditional buffering through a customized user-written function that allows the b value to change depending on the prior draw rather than remain a static value across the parameter range.
In consideration of this wide range of potential applications, we emphasize that as in any modelling exercise, iterative exploration is likely necessary with MULTI-DICE and should be embraced when it is required. We anticipate that MULTI-DICE will be a valuable and convenient tool for comparative population geneticists and phylogeographers.