An empirical examination of sample size effects on population demographic estimates in birds using single nucleotide polymorphism (SNP) data

Study system

We used eight datasets of ultraconserved elements (UCEs) from Beringian birds from McLaughlin et al. (2020; Table 1). Genomic data were generated using the methods outlined in Winker et al. (2019) using protocols from Dr. Travis Glenn’s lab at the University of Georgia (http://baddna.uga.edu/protocols.html). From these data, we generated repeatedly subsampled datasets at smaller sample sizes for analysis under a coalescent framework using δaδi (Gutenkunst et al., 2009). While we tried other programs, we were unable to get them to consistently run on our UCE datasets despite months of effort and over a hundred thousand hours of high-performance computing resources (e.g., jaatha, Mathew et al., 2013; IMa2p, Sethuraman & Hey, 2015). Thus, although we lack independent corroboration, we consider δaδi to be sufficient to answer the questions posed. Our empirical datasets represent taxonomically designated levels of population, subspecies and species pairs in three avian orders, contrasting pairs of Asian and North American populations of: Clangula hyemalis (long-tailed duck), Anas crecca crecca/A. c. carolinensis (green-winged teal), and Mareca penelope/M. americana (Eurasian and American wigeons) in Anseriformes; Numenius phaeopus variegatus/N. p. hudsonicus (whimbrel) and Tringa brevipes/T. incana (gray-tailed and wandering tattlers) in Charadriiformes; and Luscinia svecica (bluethroat), Pinicola enucleator kamschatkensis/P. e. flammula (pine grosbeak) and Pica pica/P. hudsonia (Eurasian and black-billed magpies) in Passeriformes. These datasets, which span divergence levels from populations with substantial levels of gene flow to effectively reproductively isolated species (albeit with low gene flow), enable us to explore how the effects of low sample sizes on demographic inference play out across these levels of divergence. Insofar as taxonomy is not a reliable indicator of genomic divergence levels (Humphries & Winker, 2011), we also include in our evaluations estimates of F_ST made from the full datasets (Table 1). Among the lineages in this study, pairwise comparisons fell out into two general groups, one with relatively low divergence and one with relatively high divergence (McLaughlin et al., 2020; Table 1).

These datasets consist of one single nucleotide polymorphism (SNP) per locus from over 1,500 UCE loci per lineage (each lineage is a pairwise, two-population sample of diverging populations, subspecies, or species). For bioinformatics methods, see Winker et al. (2019) and McLaughlin et al. (2020); a summary of our pipeline is given here: https://github.com/jfmclaughlin92/beringia_scripts. Each dataset consists of 100% coverage for all individuals (all individuals have phased, high-quality SNPs called at both alleles for all loci). Z-linked loci were removed because they have a different inheritance scalar from autosomal loci (Winker et al., 2019; McLaughlin et al., 2020). Original sequence data are deposited in the NCBI Sequence Read Archive (SRA; PRJNA393740). Complete data files analyzed for this study are available at https://doi.org/10.6084/m9.figshare.12622658.v1.

Subsampling datasets and analyses

To produce datasets of varying sample sizes, stepping down from the maximum number of individuals available for each population (6–8) to 1 individual per population, a custom Python script (https://github.com/jfmclaughlin92/beringia_scripts) was used. This script (ngapi_dadi.py) iteratively sampled individuals without replacement from the thinned .vcf files, created new .vcf files containing these individuals, converted these files to the proper δaδi input format (using a Perl script by Kun Wang, https://groups.google.com/forum/#!msg/dadi-user/p1WvTKRI9_0/1yQtcKqamPcJ), and ran δaδi models with predetermined, lineage-specific best-fit parameters for the split-migration (divergence-with-gene-flow) model that comes with the δaδi Demographics2D.py file (split-mig). For six of our eight lineages, split-migration models produced a best-fit model among multiple options, while for two of them a secondary contact model was a demonstrably better fit (Clangula hyemalis and Mareca penelope/americana; McLaughlin et al., 2020). Here we chose to include all eight datasets under a single model framework (split-migration, an isolation-with-migration model in δaδi, termed “split-mig” therein). We wished to focus here on changes due to sample size variation with multiple empirical datasets and not on more subtle variation due to differences among divergence-with-gene-flow models.

For each sample size, 25 subsampled datasets were created, which were each run five times. The best-fit run by highest maximum log composite likelihood value among those five runs was then selected for each dataset and used for subsequent analyses. Parameter estimates for effective population size (ν₁ and ν₂), migration (m), divergence time (T), and Θ (defined as 4N_refμ, with N_ref defined as ancestral population size and μ as mutation rate per generation), were then compared across different sample sizes. Raw parameter estimates are used throughout; we did not convert these values to individuals or years (except for three illustrative examples for individuals; see below) because that would introduce lineage-specific idiosyncrasies (e.g., through application of different mutation rate estimates) that would diminish the power of our focal among-lineage comparisons here. For the three exemplar cases in which we translated raw values into numbers of individuals, we used estimates of mutation rate and generation time given in Table S2.

The scaled root mean square error (SRMSE) was calculated, defined as $$SRMS{E}_{\theta }={\displaystyle \frac{\sqrt{{\displaystyle \frac{\mathrm{\Sigma }{\left(\hat{\theta }-\theta \right)}^2}{n}}}}{\underset{\_}{\theta }}}$$with $\theta$ in this context representing the estimate from the full dataset, $\hat{\theta }$ as the parameter estimate from the subsampled dataset, and n the number of datasets (25) considered, following Robinson et al. (2014). This was scaled by the mean of the parameter estimate at each sample size ($\underset{\_}{\theta })$ to enable inter-lineage comparisons of the changes in accuracy at lower sample sizes (SRMSE). This allowed us to quantify the changes in accuracy of estimates at different sample sizes relative to each species’ parameter estimates’ means.

Estimates of migration

We found the most variation in estimates of m occuring when samples were at 2:2 (e.g., Pinicola enucleator and Pica pica/hudsonia; Fig. 2). In most of the cases in which extreme estimates occurred at 2:2, pairings of individuals that caused geographic clustering of within-continent population samples were involved together with numerically very small estimates of m. The values of m were consistently small, but variation around the mean estimate was apparently magnified by more subtle within-continent variation than our study was designed to detect. Biologically, we reason that small values of m are the more informative takeaway, and that increased variation around those very small numbers at N of 2:2 is an artifact arising from a combination of relatively deep divergence and very low gene flow, probably coupled with some more subtle population structure within continental populations. In biological terms, although these variations can appear graphically substantial (Fig. 2, m), in Pinicola enucleator they represented estimates ranging (max–min) from 0.01 to 6.13 × 10⁻⁹ individuals per generation. In Pica, these max–min values were 0.03–2.29 × 10⁻⁹ individuals per generation.

Estimates of population size

The effective population sizes (ν parameters) were not as robust, with variation tending to begin to increase markedly below four diploid individuals per population and accuracy decreasing in all lineages (Tables 2 and 3; Fig. 1; Figs. S1 and S2). They were, however, still reasonably accurate in many lineages at relatively small samples sizes (Tables 2 and 3; Fig. 1; Figs. S1 and S2). The negative relationships between SRMSE values for each demographic parameter and sample size (N) should help users interpret how lineages and individual parameters are affected by smaller sample sizes (Table 4; Table S1).

The impact of divergence

Our results reinforce previous findings (Kalinowski, 2005; Morin, Martien & Taylor, 2009) that an important factor in determining the minimum sample size for a study is the level of divergence in the lineages under examination. Although this might be known at the start of a study, that might not always be true, potentially complicating sampling design. However, some general recommendations are possible, at least within a broader framework of higher-and lower-divergence groups. Lineages with considerable divergence (e.g., species-level, such as in Tringa) had accurate demographic parameters estimated at lower sample sizes (Fig. 2; Figs. S1–S5). Thus, it seems possible in such systems to reliably use fewer individuals. In shallowly diverged populations that might experience substantial gene flow, however, higher sample sizes may be required to overcome the impact of individuals with varying amounts of admixture, which appears to increase the variation in model performance at lower sample sizes among low-divergence lineages (Fig. 1; Figs. S1–S6; Table S1).

Our findings of the effects of divergence levels on the minimum sample sizes needed to accurately estimate population demographic parameters broadly agreed with previous findings in other genetic markers, with some exceptions. In lineages that are more shallowly split and have experienced more gene flow, greater sample sizes are required to reliably estimate multiple parameters, including not just the demographic parameters examined here, but also genetic distance (Kalinowski, 2005), F_ST (Morin, Martien & Taylor, 2009; Humphries & Winker, 2011), and recent demographic events (e.g., <100 generations; Beichman, Huerta-Sanchez & Lohmueller, 2018). The two population-level splits in our study, L. svecica and C. hyemalis, did not perform as well for most parameter estimates at sample sizes below 6 individuals per population, with accuracy (as measured by SRMSE; Table 3) decreasing rapidly; this fits our understanding that accurately estimating more recent demographic events requires the improved draw on more recent coalescent events that increased sample sizes bring (Beichman, Huerta-Sanchez & Lohmueller, 2018). The presence of a substantial amount of gene flow appears to increase variation in parameter estimates and decrease accuracy, as seen in L. svecica (Tables 2 and 3), and in practical terms would require increased sample sizes for accurate parameter estimation.

Model fit

Due to computational restrictions, we analyzed all subsampled datasets under the split-mig δaδi (Gutenkunst et al., 2009) model determined and optimized for the full dataset in each lineage, and we did not investigate the impact of sample size on model fit. Several subsample datasets (notably Clangula hyemalis, Mareca penelope/M. americana, and Luscinia svecica) showed signs in some parameters of beginning to consistently push the upper bounds of some model parameters. This means that both variation and over-estimation of the parameters were likely underestimated in these groups at smaller sample sizes. This situation has also been noted with simulated data, which have been found in some situations to have a better fit with a model type different than the one under which they were simulated (Robinson et al., 2014).

Implications for study design

Research efficiency requires attention not only to the minimum sample size required to meet an objective, but also to the point after which adding more samples begins to produce diminishing returns. In this context, this means the point above which the SRMSE becomes similar between sample sizes, but before the means of estimates start to change due to decreased sample size. This inflection point might represent the minimum reliable sample size, but not necessarily. In some lineages, SRMSE was very similar at larger sample sizes, began to slowly increase at intermediate sizes, and then at low sample sizes increased quickly (Table 3; Fig. 2). This again varied among lineages (Table 3; Figs. S1–S5). In some, such as the Pica and Tringa species lineages, this inflection point was reached at higher sample sizes than the minimum reliable sample sizes in some parameters (Table 3), whereas in others, such as in most estimates of m, these points were the same (e.g., Fig. S3). However, in some groups, particularly estimates of effective population size (ν₁) and migration (m) in L. svecica, this optimal point was not reached until the full dataset was analyzed, and might not have been reached at all in C. hyemalis in any of the parameter estimates (Figs. S1–S5). This is consistent with the findings of Robinson et al. (2014), in that although in some cases a small sample size could be used, larger sample sizes still led to more accurate parameter estimates. This was especially the case in our data for divergence times (T), Θ, and some effective population size (ν) estimates (Table 2; Fig. 1; Figs. S1–S5). Our linear regression models help generalize these relationships (Table 4; Table S1). In sum, two key sources of variation preclude our providing detailed suggestions for threshold sample sizes in future studies: the levels of divergence in a study’s focal lineage, and the demographic parameter of most interest for that study. We urge those designing their own studies to consider our results (Tables 2–4, Supplemental Information) at different divergence depths and for different demographic parameters, depending on objectives.