Substructured Population Growth in the Ashkenazi Jews Inferred with Approximate Bayesian Computation

Abstract

The Ashkenazi Jews (AJ) are a population isolate sharing ancestry with both European and Middle Eastern populations that has likely resided in Central Europe since at least the tenth century. Between the 11th and 16th centuries, the AJ population expanded eastward leading to two culturally distinct communities in Western/Central and Eastern Europe. Our aim was to determine whether the western and eastern groups are genetically distinct, and if so, what demographic processes contributed to population differentiation. We used Approximate Bayesian Computation to choose among models of AJ history and to infer demographic parameter values, including divergence times, effective population sizes, and levels of gene flow. For the ABC analysis, we used allele frequency spectrum and identical by descent-based statistics to capture information on a wide timescale. We also mitigated the effects of ascertainment bias when performing ABC on SNP array data by jointly modeling and inferring SNP discovery. We found that the most likely model was population differentiation between Eastern and Western AJ ∼400 years ago. The differentiation between the Eastern and Western AJ could be attributed to more extreme population growth in the Eastern AJ (0.250 per generation) than the Western AJ (0.069 per generation).

population genetics, Approximate Bayesian Computation, substructure, demography

Introduction

As we enter the age of genomic medicine, it is vital to fully understand the complex demographic factors that shape patterns of genetic variation. Without a better understanding of population dynamics, the promise of medical genetics cannot be fulfilled. Hidden population substructure characterized by genetic differences among closely related populations can affect results of association (Martin et al. 2017). Substructure resulting from different demographic histories can result in errors in estimates of variant frequencies. For example, rapid population growth causes an excess of rare variants, increasing the frequency of deleterious mutations in a population (Keinan and Clark 2012), whereas gene flow increases the frequency of common variants in a population. Knowledge of an individual’s ancestry and the demographic processes that contributed to their genomic architecture will help improve accuracy of medically relevant inferences.

The Ashkenazi Jewish (AJ) communities from Western/Central (WAJ) and Eastern Europe (EAJ) are an ideal population for studying the effects of complex demography because they have a well-documented history (Malamat and Ben-Sasson 1976) with census records (DellaPergola 2001), they have been relatively isolated, and they are known to have a large number founder mutations associated with disease (Einhorn et al. 2017). The AJ have experienced a complex demographic history, including bottlenecks, extreme growth, and gene flow, all of which can affect the incidence and prevalence of genetic diseases (see Gladstein and Hammer 2016 for a more in-depth review of the literature and Carmi et al. 2014; Xue et al. 2017 for recent demographic inference). Different frequencies of mutations associated with diseases have been found among AJ communities (Risch et al. 2003); however, most medical studies treat the AJ as a single isolated population (see for example Einhorn et al. 2017). The historical record suggests there may have been a substantial difference in growth rates between the communities in Western/Central and Eastern Europe (DellaPergola 2001). Some authors suggest that the EAJ population size is unlikely to have reached observed 20th century levels without massive influx from nonJewish Eastern Europeans (Straten 2007) or from the Khazars, a lost Turkic population that converted to Judaism (Brook 2018). However, no genetic work has examined structured growth within the AJ and whether the recorded population growth in Eastern Europe is plausible without major external contributions. If left unresolved, this gap in our understanding of AJ population history could lead to deficient conclusions from medical studies.

The eastward expansion of the AJ settlement started after the 11th century and continued into the 16th century due to religious and ethnic persecution in Western and Central Europe (Hundert 2004). This led to the establishment of two major culturally distinct communities in Western/Central and Eastern Europe (Meyer et al. 1996). This is reflected in the division between the Eastern and Eastern Yiddish dialects (Katz 2004; supplementary fig. S1, Supplementary Material online). Previous genetic studies have produced conflicting results on whether there is substructure within the AJ (Karlin et al. 1979; Livshits et al. 1991; Risch et al. 2003; Behar, Garrigan, et al. 2004; Behar, Hammer, et al. 2004; Klitz et al. 2010; Guha et al. 2012; Granot-Hershkovitz et al. 2018), and it is currently unclear whether AJ cultural differences correspond to genetic substructure.

Our aim was to determine whether the WAJ and EAJ are genetically distinct subpopulations, and if so, whether differential population growth and/or gene flow contributed to the differentiation (fig. 1). We performed coalescent-based Approximate Bayesian Computation (ABC) to infer the most likely model of AJ history based on genome-wide SNP data from 258 AJ with all four grandparents known to belong to either western (n = 19) or eastern (n = 239) cultural groups. With current computational capabilities we were able to perform millions of genomic simulations at an unprecedented chromosome-size scale, allowing us to use identical by descent (IBD) segments, which are informative for more recent demographic inference, and allele frequency based-statistics, which are informative for older inference (Carmi et al. 2013; Ralph and Coop 2013; Xue et al. 2017). Previous methods using ABC used small loci, making chromosome-wide haploblock identification impossible, which is necessary for recent demographic inference. In addition, we used an ABC approach that takes into account SNP array ascertainment bias (Quinto-Cortés et al. 2018), which must be taken into account when performing demographic inference on populations not represented in SNP discovery panels. Without correcting for ascertainment bias, allele-frequency spectra tend to be skewed toward a deficiency of rare alleles and an excess of common alleles (Quinto-Cortés et al. 2018). In this study, we pushed the limit of computing capabilities by simulating whole chromosomes of hundreds of individuals, paving the way for future model-based inference of recent demographic history.

Fig. 1.

Three modeled demographic histories of the AJ. The Null Model has no substructure within the AJ. The Substructure model has a population split between Eastern and Western AJ and one common admixture event from Europeans. The Substructure with Differential Admixture model has a population split between Eastern and Western AJ and separate admixture events from the Europeans. Populations are labeled as E, European; AJ, Ashkenazi Jews; EAJ, Eastern Ashkenazi Jews; WAJ, Western Ashkenazi Jews; J, Sephardic/Mizrahi Jews; ME, Middle Eastern.

Results

We first used Principal Component Analysis (PCA) and ADMIXTURE to look for qualitative signatures of population structure. Next we tested for a genetic difference between Eastern (EAJ) and Western (WAJ) AJ using a variety of population genetic measures. Finally, we used ABC to choose the best model out of the three hypothesized demographic histories and inferred the parameters of the best model.

The PCA showed a similar pattern as other studies, with the AJ forming a cluster between European and Middle Eastern populations on PC1, with the AJ separating from the other populations on PC2 (Behar et al. 2010; Bray et al. 2010; Granot-Hershkovitz et al. 2018). Additionally, there was separation of the EAJ and WAJ, with the WAJ further on PC2 from both the European and Middle Eastern clusters than the EAJ (fig. 2), as in Granot-Hershkovitz et al. (2018). Drop one in PCA and ADMIXTURE revealed similar structure as PCA, showing differentiation between EAJ and WAJ (supplementary figs. S6–S8, Supplementary Material online).

Fig. 2.

PC1 versus PC2 with Jewish, European, and Middle Eastern populations.

EAJ and WAJ were significantly different in pairwise F_ST (P < 0.0050), PCA mean coordinate values (ANOVA, $P=8.5e-22$⁠), and admixture proportions (Wilcoxon, $P=2.803e-08$⁠; supplementary figs. S10 and S11, Supplementary Material online). Additionally, the WAJ had significantly larger F_IS values than EAJ (Wilcoxon, $P=4.204e-13$⁠), more frequently had long ROH (>1.8 Mb) than the EAJ, and when the ROH length cutoff was 1 Mb, the WAJ had significantly larger median ROH than the EAJ (Kruskal–Wallis, P = 1.572e–08; supplementary figs. S12–S14, and table S8, Supplementary Material online). These results are consistent with a smaller effective population size in WAJ and/or greater gene flow in the EAJ from external populations.

We investigated whether the EAJ genomic archetype could result from Khazarian ancestry. We performed an ad hoc resampling simulation by using the WAJ data as a proxy for ancestral EAJ and randomly perturbing the alleles so that the allele frequencies matched various amounts of gene flow from the Turkic-speaking Chuvash host, a proxy for the Khazars. If the EAJ have substantial Khazarian genetic contribution, we expect them to appear more like the EAJ as we increase the amount of Chuvash ancestry in the WAJ. We found that adding Chuvash gene flow to the WAJ did not create an admixed population similar to EAJ (supplementary fig. S20, Supplementary Material online).

Approximate Bayesian Computation

In order to insure our ABC analysis would be sensitive to both recent and more ancient history we used a combination of AFS and IBD summary statistics. Recombination breaks up IBD segments over time. Therefore, large IBD segments are expected to be informative of recent history, including effective population size, population structure, and gene flow (Carmi et al. 2013; Ralph and Coop 2013; Xue et al. 2017).

Model Choice

To test whether EAJ and WAJ diverged into two subpopulations, and if so, whether different histories of population growth or gene flow contributed to the differentiation, we chose the best out of three demographic models (fig. 1) with ABC. Our models incorporated population structure, growth, and gene flow. Our first model represents the null hypothesis of cultural differentiation without genetic division between EAJ and WAJ. This could be due to continuous gene flow between EAJ and WAJ, insufficient time for genetic drift to cause differing allele frequencies, and/or similar rates of gene flow with Europeans in Eastern and Central Europe. Next, the Substructure model includes a population split between EAJ and WAJ allowing for different population sizes, following one common gene flow event with Europeans after the initial founder event. This could reflect low rates of gene flow between WAJ and EAJ and/or different growth rates (with similar levels of gene flow from Europeans) after their split. Finally, the Substructure with Differential Admixture model is similar to the substructure model; however, it allows for gene flow with Europeans separately in EAJ and WAJ after the initial divergence of WAJ and EAJ.

We found that the Substructure model was the best fitting model with a Bayes factor of 2.01 and posterior probability of 0.67 (supplementary tables S9 and S10, Supplementary Material online). Our observed Bayes factor of the substructure model was >92% of the Substructure Model Bayes factors from the cross validations when the Null Model was the true model and >86% of the Substructure Model Bayes factors from the cross validations when the Substructure Model with Differential Admixture was the true model (supplementary fig. S16, Supplementary Material online). Therefore, given that the Substructure Model was found to be the best model, it is unlikely that the Null Model or the Substructure with Differential Admixture Model is the true best fit.

Parameter Estimates

The estimations for the Middle Eastern, Jewish, and AJ populations are of particular interest to this study (table 1, fig. 3). We used a generation time of 25 years to convert from generation to years. Since for the simulations we used a mutation rate of 2.5e–8, all parameter estimates were given in terms of that mutation rate. The parameter estimates are given in table 1. Parameter estimates from additional ABC analyses on the whole genome is given in supplementary tables S12 and S13, Supplementary Material online.

Table 1.

Prior Distributions and Parameter Estimates of AJ Demographic Parameters of the Substructure Model from 1,446,125 Simulations of Chromosome 1, with 10 PLS Components and 1,000 Retained Simulations.

Parameter	Description	Prior Distribution	Distribution	Mode	HPDI 95
log10(NWA)	Ne in Western AJ	(3.0,6.7)	Uniform	3.90	(3.10,6.48)
log10(NEA)	Ne in Eastern AJ	(4.0,6.7)	Uniform	6.18	(4.24,6.7)
log10(NAg)	Ne in ancestral AJ, before growth	[2.0, log10(NWA)]	Uniform	3.04	(2,4.33)
		or [2.0, log10(NEA)]
log10(NJ)	Ne in nonAJ Jews	(3.0,6.0)	Uniform	5.61	(4,6)
log10(NM)	Ne in Middle Easterners	(3.0,6.0)	Uniform	5.67	(4.04,6)
m	Geneflow from Europeans	(0,1)	Uniform	0.16	(0.015,0.94)
TEM	European—Middle Eastern split time	(400,1200)	Discrete uniform	480.81	(400,1007)
TMJ	Middle Eastern—Jewish split time	(TA,TEM)	Discrete uniform	211.38	(21,706.96)
TA	AJ—Jewish split time	(20,36)	Discrete uniform	28.73	(20,36)
TAEW	Eastern—Western AJ split time	(2,TA)	Discrete uniform	14.93	(2,29)

Note.—Divergence time is given in generations and the proportion of gene flow is given as the fraction of AJ that comprised European migrants each generation. See supplementary figure S3, Supplementary Material online for depiction of model parameters. See supplementary table S3, Supplementary Material online for priors of noninferred parameters.

Fig. 3.

Posterior, truncated prior, and prior distributions of parameters of the Substructure Model from 1,446,125 simulations of chromosome 1. Grey area is the simulated prior, dotted line is the truncated prior (1,000 simulations closest to the observed data), and the solid line is the approximated posterior. Black solid vertical line is the mode of the posterior density and the dotted vertical lines are the lower and upper limit of the 95 high posterior density interval. Units of time and migration rate are as given in table 1. See supplementary figure S3, Supplementary Material online for depiction of model parameters.

We found the probability that the EAJ effective population was greater than the WAJ effective population size from the joint posterior of the two parameters to be 0.69 (fig. 4). Given the estimated effective population sizes and divergence time from other Jewish populations, we estimated the exponential growth rate in EAJ and WAJ to be 0.25 and 0.069 per generation, respectively.

Fig. 4.

Joint posterior of effective population size in Eastern and Western AJ. The brighter colors indicate a higher posterior probability. The diagonal black line indicates where N_e in Western and Eastern AJ are equal. Sixty nine percentage of the joint density is greater in Eastern AJ than in Western AJ.

To compare our inferred growth estimates from ABC to census data, we found census data of Jews in Europe from the 12th to 20th centuries (DellaPergola 2001; supplementary fig. S18 and table S11, Supplementary Material online). However, the source for the census data were not specific to AJ, and grouped Jews in Western and Southern Europe together, likely including Sephardi Jews from Southern Europe and the Balkans. Thus, we expect that the estimates of the population sizes for the Jews of Western and Southern Europe is an overestimate of the population sizes of WAJ. We fit an exponential model to the census data from the years 1170 to 1900. We found the same order of magnitude in EAJ and high concordance in WAJ for our estimates of growth rates based on genomic (0.25 per generation in EAJ and 0.069 per generation in WAJ) and census data (0.18 ± 0.017 per generation in Eastern European Jews and 0.081 ± 0.026 per generation in Western/Southern European Jews). Coventry et al. (2010) estimated similar growth rate in Europeans from census data (0.115 per generation since 1600; Livi-Bacci 2001), and a 10-fold higher growth rate in Europeans from genetic data (1.094 per generation).

Discussion

AJ traditionally trace their origin to ancient Hebrews who lived as semiherders in the Levant over 3000 years ago and have experienced a long subsequent history of migrations. The first known written accounts of “Israel,” in the central hill country of the southern Levant, is from the Merneptah Stele in 1207 BCE (Hjelm 2016). The destruction of the First (587 BCE) and Second Temple (70 CE) contributed to the establishment of major Jewish centers outside of Israel (Gruen 2002). While it is not clear what Jewish movement took place in Europe between the first and fourth centuries, migrations northward from Italy led to an established AJ community in the Rhine Valley by the tenth century (DellaPergola 2001). In the late Middle Ages and early modern period, Jews were expelled from much of Central Europe, and the Polish kings and nobility invited Jews to Poland (Hundert 2004).

From the end of the 16th century two major culturally distinct communities developed in Central and Eastern Europe (Meyer et al. 1996), each speaking a distinct Yiddish dialect (Katz 2004) and following different folk customs (e.g., dress, food, liturgical melodies, superstitions; Meyer et al. 1996). Throughout the Middle Ages in Central Europe, Jews were often expelled from their settlements and there were strict regulations on where Jews could live and what they could do to earn a living (Meyer et al. 1996), whereas in Eastern Europe Jews could generally move freely and were protected by nobles, which fostered a sense of belonging (Hundert 2004). In the 19th century, Jews in Central Europe became more integrated in general life, with less focus on traditional and religious institutions (Meyer et al. 1996), whereas Jews in Eastern Europe maintained Jewish learning and institutions and for the most part their lives revolved around Jewish religious tradition (Hundert 2004).

Our aim was to determine whether AJ with ancestry tracing to eastern and western cultural groups were genetically subdivided, and if so, what contributed to rapid population differentiation (fig. 1). Our ABC analysis based on AFS and IBD statistics from extensive simulations of chromosome 1 rejected the null hypothesis of no substructure in the AJ, and favored the Substructure Model over the Substructure Model with Differential Admixture (supplementary table S9, Supplementary Material online). While several factors could contribute to subdivision, we found that differential growth rates in the two subpopulations were a major factor. Our results support a model in which the AJ started from a small founder population upon initial arrival in Europe, followed by moderate population growth in Western/Central Europe and massive population growth in Eastern Europe. Our conclusions from ABC are consistent with the observation of longer ROH and IBD and higher F_IS in the WAJ, which could reflect a smaller N_e than in EAJ. The ADMIXTURE plot in supplementary figure S8, Supplementary Material online is also consistent with more genetic drift in the WAJ, which appears to have a higher proportion of an “AJ” component. Greater genetic drift may also be the explanation for the more extreme position of the WAJ in the PCA (fig. 2).

The estimated census population growth in the EAJ has been referred to as the “Demographic Miracle” (Malamat and Ben-Sasson 1976), and has been the center of much debate with some arguing that it is unfeasible without massive conversion or intermarriage (Straten 2007). However, because the Substructure Model was favored over the Substructure with Differential Admixture Model, we did not find that gene ow from Europeans was necessary to increase the effective population size in the EAJ relative to WAJ. Our genetically based results corroborate DellaPergola (2001)’s conclusion based on demographic analysis of census data that the rapid growth of EAJ is feasible without mass immigration or large-scale conversions to Judaism (at least from nonJewish Europeans). However, unmodeled admixture from a ghost population could cause an increase in effective population size.

The Khazars have been hypothesized to have made a large genetic contribution to EAJ through a process of admixture and conversion (supplementary fig. S19, Supplementary Material online). Historical documents suggest that the Khazars were a Turkic people whose empire spanned parts of modern day Ukraine, southern Russia, and the Caucasus from the seventh to 11th century (Brook 2018). Although some documents suggest that Khazarian royalty converted to Judaism for political reasons, it is unknown how much of the Khazarian population converted (Brook 2018). Elhaik (2013) and Behar et al. (2013) have both addressed the possibility of substantial Khazarian contribution into an unstructured AJ gene-pool. Elhaik (2013) favored of the Khazarian admixture/conversion hypothesis whereas Behar et al. (2013) found no evidence to support this hypothesis based on a large survey of SNPs in Jewish and nonJewish population samples. While we did not explicitly test a model of Turkic gene flow to the EAJ in our ABC analysis, our results indicate that the EAJ could have undergone the observed growth without contribution from outside populations, such as the Khazars.

In addition, we performed an ad hoc resampling simulation in which we assumed that our WAJ data represented a reasonable proxy for the ancestral EAJ, and then randomly mixed WAJ genes with those of a sample from the contemporary Chuvash population. The Chuvash are a Turkic-speaking population that currently live near the Volga River in a region that was part of the ancient Khazar empire. Our simulation results did not support the Khazarian admixture model because the simulated population did not appear more like the EAJ as we increased the amount of Chuvash ancestry in the WAJ (supplementary fig. S23, Supplementary Material online). A major is that our interpretation of the simulation depends on the assumption that the level of genetic drift in the WAJ was low enough that they are a reasonable proxy for ancestral EAJ. Interestingly, census data show that the EAJ experienced high growth rates relative to Western/Southern Jews in the 16th and 20th centuries, well after the fall of Khazariaa result that is not consistent with a large Khazarian contribution to the EAJ (DellaPergola 2001).

An alternative explanation to account for the higher EAJ growth rate involves more conducive political, economic, and social conditions in Eastern Europe (Meyer et al. 1996). In Central Europe, discriminatory legislation of the late Middle Ages and early modern period meant that very few Jews lived there (Meyer et al. 1996). There were often legal limitations on the number of Jewish families—regardless of where they lived (Meyer et al. 1996). In the 19th century cramped ghettos and integration into nonJewish society contributed to low birthrates (Bachi and DellaPergola 1984). Unlike in other regions of Europe, Jewish communities in Eastern Europe neither had limitations imposed by the government on the number of Jewish marriages, nor strict residence restrictions (Bachi and DellaPergola 1984). Adherence to religious and traditional norms and the economic structures encouraged early marriage and high fertility in Eastern Europe (Bachi and DellaPergola 1984). Additionally, generally low mortality rates among the AJ, particularly low infant mortality, contributed to the higher growth rate.

Implication for Genetic Diseases

The AJ have an unusually high prevalence of >50 known disease-associated mutations (Consortium 2011), primarily because of founder effects resulting from population bottlenecks (Risch et al. 2003; Slatkin 2004). Both Risch et al. (2003) and Slatkin (2004) found evidence of founder events ∼11 and ∼5 centuries ago. Slatkin (2004) tested for neutrality and founder effect in several disease-associated alleles found predominantly in AJ. He used LD with a linked marker allele and a historical demographic model, assuming no subdivision. Additionally, Risch et al. (2003) examined allele-frequency distributions and estimated coalescence times of several disease-causing mutations at increased frequency in the AJ. Although the older founder mutations were present across all AJ subpopulations, the more recent founder mutations were restricted to Lithuanians (Risch et al. 2003). As far as we are aware no previous work has modeled population structure and population size changes in the AJ to examine the effects these processes have had on the distribution and frequency of disease alleles.

Different growth rates in WAJ and EAJ may have varying effects on the frequency of deleterious mutations in these subpopulations. Extreme recent population growth from a relatively small ancestral population increases the likelihood of new mutations on a background of high homozygosity. It has been suggested that the high incidence and frequency of deleterious founder mutations in the EAJ resulted from strong growth of the Lithuanian Jewish community, as seen in the census data between the 14th century and the 1800 s (Meiner et al. 1991; Risch et al. 2003). However, it is unclear the extent to which continuous small population size (e.g., as in the WAJ), versus founder effect with extreme population growth (e.g., as in the EAJ), has influenced increased frequency of disease alleles in the AJ.

Because we used a combination of AFS and IBD based statistics, we expect our inferences of older demographic events (TEM and TMJ) to be influenced more by the mutation rate, whereas recent events should be influenced more by recombination rate. Therefore the mutation rate used in the simulation should not affect the majority of the inferred parameters. However, the parameter estimates can be adjusted to the 1.2e–8 (Kong et al. 2012) or 1.44e–8 (Gravel et al. 2013) mutation rate by doubling the estimated divergence times and effective population sizes.

Our estimate of the Middle Eastern and European divergence time corresponds with other published estimates (Haber et al. 2013; Carmi et al. 2014).Our estimate of the time of divergence of the AJ and nonAshkenazi Jews (AJ) overlaps with estimates of a population size reduction of the AJ (Carmi et al. 2014). Our estimate of current N_e in the EAJ is in agreement with that of Carmi et al. (2014) while our estimate of the founder effect size is slightly larger with overlap of the 95 high posterior density interval.

Although we have used a generation time of 25 years (as in Haber et al. 2013; Carmi et al. 2014), longer generation times may be more reasonable given the recent timescale being investigated (as in Ralph and Coop 2013). If we assume a generation time of 30 years, then the divergence times in years will be pushed back, giving better concordance with historical events in the recent past.

Previous genomic scale demographic studies did not take into account population substructure in their sampling strategies (Palamara et al. 2012; Carmi et al. 2013, 2014; Xue et al. 2017). While our conclusion of substructure in the AJ is consistent with results from studies using haploid loci and classical markers (Karlin et al. 1979; Livshits et al. 1991; Risch et al. 2003; Behar, Garrigan, et al. 2004; Behar, Hammer, et al. 2004; Klitz et al. 2010), it is contrary to Guha et al. (2012) conclusion that the AJ are not genetically subdivided. However, Guha et al. (2012)’s PCA results of AJ with ancestry from different countries does not provide support for or against our conclusion of subdivision between the EAJ and WAJ, because they do not use historically/culturally motivated AJ groupings.

Xue et al. (2017) inferred an admixture event with Southern Europeans ∼24–49 generations ago, and at least one admixture event ∼10–20 generations ago; however, they were not able to infer whether this more recent event involved western or eastern nonJewish Europeans. Although the uncertainty could be due to limitations in their methods to distinguish between Western or Eastern European sources, it is plausible that their uncertainty was due to their use of the AJ as a single population. Potentially, if they had performed their analyses separately with EAJ and WAJ, they would have been able to more clearly identify the sources of the later admixture events.

Caveats and Future Directions

We attempted to keep the model as simple as possible, while still providing insight into substructure in the AJ. Because we did not include European gene flow to Sephardic Jews, we may have overestimated the divergence time between the Sephardic Jews and Middle Eastern population (e.g., as a result of European lineages in the Sephardic Jews). This may also have contributed to a more recent estimate of divergence time between the AJ and Sephardic populations and to an overestimate of the Sephardic N_e. Additionally, we ignored possible gene flow between EAJ and WAJ to avoid loss of power. If significant gene flow actually occurred, we may have overestimated the divergence time between the EAJ and WAJ. Additionally, while the Substructure Model includes growth after EAJ diverged from WAJ, we did not explicitly test whether growth started before or after divergence.

While we a priori attempted to choose summary statistics that would capture older and more recent history, it is possible that our set of summary statistics were not informative for all aspects of the model. In particular, for the Substructure Model we inferred a surprisingly low proportion of European gene ow into the AJ. In comparison, recent studies have estimated ∼50% admixture from Europeans (Carmi et al. 2014; Xue et al. 2017). We believe that our combination of summary statistics may not have been informative for European gene flow, and we are not confident in the admixture proportion estimates returned in this study.

Regarding AJ demographic history we would like future studies to address three topics—more detail regarding population growth, the source of small effective population size, and more complex admixture history. Future studies should compare models of instantaneous, exponential, and logistic growth, and attempt to distinguish between genetic drift caused by small N_e in a randomly mating population versus preferential mating of more closely related individuals. Finally, future studies could use more sensitive local ancestry estimation with multiple European sources and multiple admixture events, with distinction between multiple events and continuous admixture.

Materials and Methods

Data Set

High density SNP array in Jewish, European, and Middle Eastern populations, and whole genome sequence data in world-wide populations were used (supplementary tables S1 and S2, Supplementary Material online). An additional 231 AJ samples were genotyped on Illumina Omni Express chip for 730 K markers. AJ samples were collected in the United States and Israel and information on the grandparental country of origin was provided, with a total of 13 European countries. We defined the AJ as Eastern or Western based on their grandparental country of origin, according to Yiddish dialect borders (supplementary fig. S1, Supplementary Material online). All included samples had four AJ grandparents from Eastern or Central Europe (we did not include AJ individuals of mixed Eastern and Central European ancestry).

Approximate Bayesian Computation

We used ABC to choose the best supported demographic model and to obtain the posterior distribution of model parameters values under the best-supported model. We used the ABC-GLM (General Linear Model) approach introduced by Leuenberger and Wegmann (2010) and implemented in the program ABCtoolbox 2.0 (Wegmann et al. 2009). An overview of the ABC pipeline is shown in figure 5. The following sections provide more details on each of the ABC steps (reproducible step-by-step instructions are available on the GitHub wiki, https://github.com/agladstein/AJ_ABC/wiki; last accessed March 18, 2019).

Fig. 5.

Overview of ABC pipeline. The eight main steps of the pipeline with the tools used to perform them (Chen et al. 2009; Wegmann et al. 2009; Mevik et al. 2016; Gladstein 2018).

Demographic Models

We simulated three models of AJ demographic history, with an incorporated SNP array ascertainment scheme based on the three HapMap populations, CEU, CHB, and YRI (Quinto-Cortés et al. 2018; supplementary figs. S2–S4, Supplementary Material online). We randomly generated all the model parameter values for the simulations according to prior distributions based on previously inferred parameter values (Gutenkunst et al. 2009; Carmi et al. 2014; Quinto-Cortés et al. 2018) and documented history (Malamat and Ben-Sasson 1976; supplementary table S3, Supplementary Material online). On the basis of the sample sizes from the observed data, we simulated 9 YRI, 9 CEU, 4 CHB, 38 AJ (or 19 Eastern and 19 Western for models 2 and 3), 14 Jewish, and 14 Middle Eastern. We reduced the simulation sample size of the AJ (specifically Eastern), Jewish, and Middle Eastern populations to reduce the run-time of the simulations. The three models are the same until the AJ diverge from the other Jewish populations. The Null Model has no population structure in the AJ, whereas in the Substructure Model and the Substructure with Differential Admixture Model there is a population split in the AJ. In the Substructure Model, gene flow from Europeans to the AJ occurs before the AJ split. In the Substructure with Differential Admixture Model, gene flow from Europeans occurs after the split in the AJ, allowing for different proportions and times of gene flow in EAJ and WAJ. Both the Substructure Model and the Substructure with Differential Admixture Model allow for independent population growth in the EAJ and WAJ. The three models are described in detail in supplementary tables S4–S6, Supplementary Material online, with their appropriate MaCS (Chen et al. 2009) commands.

Simulations

We performed ∼1 millions simulations of chromosome 1 and calculated summary statistics with SimPrily-alpha (Gladstein 2018), which uses a modified version of Markovian Coalescent Simulator (MaCS; Chen et al. 2009). We ran the simulations in parallel on high throughput clusters (Pordes et al. 2007; Sfiligoi et al. 2009; Moore et al. 2014; Towns et al. 2014; Nystrom et al. 2015; Stewart et al. 2015; supplementary fig. S5, Supplementary Material online). We used a mutation rate of μ = 2.5e–8 (Nachman and Crowell 2000) and a HapMap recombination map (Fraser et al. 2007), with the recombination rate, ρ = 1e–8.

We made pseudo SNP arrays from the simulated whole chromosome 1 based on real SNP arrays and a discovery sample. We followed the method of Quinto-Cortés et al. (2018) and used random samples from YRI, CEU, and CHB as the SNP discovery set (uniform prior of 2–20 chromosomes for each population). We used a random derived allele threshold (uniform prior of 0.05–0.1) to find SNPs to use for the pseudo array. To create the pseudo array we found the closest simulated sites to the real array that had a minor allele frequency greater than the threshold in the discovery sample. By taking into account ascertainment bias in our simulated model were able to recover rare alleles not present in the SNP array data (supplementary fig. S17, Supplementary Material online).

Summary Statistics Choice

We included AFS based summary statistics to capture information on the older parts of the demography and included IBD statistics to summarize information on recent history (supplementary table S7, Supplementary Material online). We used a total of 181 summary statistics, consisting of the number of segregating sites, number of singletons, number of doubletons, Tajima’s D, F_ST, mean and median length, variance of length, and number of shared IBD >3 Mb and 30 Mb. We used the same methods to calculate all summary statistics, including the IBD calling with GERMLINE, from the simulated and real data (see Supplementary Methods and Results section 4.2).

To choose the best model, we found the best sets of summary statistics for model choice from a pruned subset of summary statistics, using the greedy search algorithm in ABCtoolbox (see supplementary methods, Supplementary Material online).

In order to avoid correlations among statistics and to preserve the informativeness of the data, we transformed our statistics with Partial Least Squares (PLS) and used the first 10 PLS components for the parameter estimation.

Model Choice and Parameter Estimation

For model choice and parameter estimation, to maximize our power we used a reduced set of 9, 11, and 13 parameters for the three models respectively models, leaving the remaining parameters as noninferred free parameters. We compared the three models using 1,275,807 simulations with ABCtoolbox, which chooses the best model by calculating Bayes factors for each compared model. The Bayes factors are calculated from the marginal densities of the models (Wegmann et al. 2009). We then inferred the parameters of the best model with 1,446,124 simulations using ABCtoolbox, which applies a General Linear Model regression adjustment to estimate the posterior densities of parameters (Wegmann et al. 2009). We used the mode of the posterior densities as the parameter estimates.

Validation of Model Choice and Parameter Estimation

We used cross-validation with ABCtoolbox to assess the accuracy of model choice and parameter estimates of the best model. Cross-validation executes the model choice and parameter estimation process 1,000 times by choosing one of the simulated data sets as the observed data. From this, we obtained a confusion matrix with estimates of false positives and false negatives. Because of the way the cross validation works, it is necessary for all models to have the same summary statistics. Therefore, we forced the Null Model to have the same summary statistics as the Substructure models by randomly splitting the simulated AJ into two groups and calling one Eastern and one Western.

Data Availability

Data used in the study are available from the figshare public repository under the project “Substructured population growth in the Ashkenazi Jews inferred with Approximate Bayesian Computation,” at https://figshare.com/projects/Substructured_population_growth_in_the_Ashkenazi_Jews_inferred_with_Approximate_Bayesian_Computation/37856, last accessed March 18, 2019. Additionally, the code for the ABC analysis is accessible at https://github.com/agladstein/AJ_ABC, last accessed March 18, 2019.

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Acknowledgments

An allocation of computer time from the UA Research Computing High Performance Computing (HPC) at the University of Arizona is gratefully acknowledged.

We would like to thank Mats Rynge for his invaluable help setting up the Pegasus workflow and running it on the Open Science Grid. This research used resources provided by the Open Science Grid, which is supported by the National Science Foundation award 1148698, and the U.S. Department of Energy’s Office of Science.

This material is based upon work supported by the National Science Foundation under Award Numbers DBI-0735191 and DBI-1265383 and we thank CyVerse for their External Collaborative Support program.

This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. The computations were conducted on the Comet supercomputer, which is supported by National Science Foundation award number ACI-1341698, at the San Diego Supercomputing Center (SDSC), the Jetstream cloud environment, which is supported by National Science Foundation award number NSF-1445604, at Indiana University and Texas Advanced Computing Center, and the Bridges supercomputer, which is supported by National Science Foundation award number ACI-1445606, at the Pittsburgh Supercomputing Center (PSC).

This research was performed using the compute resources and assistance of the UW-Madison Center For High Throughput Computing (CHTC) in the Department of Computer Sciences.

Pegasus is funded by the National Science Foundation under the Offce of Advanced Cyberinfrastructure S12-SS1 program, grant #1664162.