A novel approach for choosing summary statistics in approximate Bayesian computation

Estimating mutation rates in a Markov branching process using approximate Bayesian computation

Estimating microbial mutation rates is an essential task in evolutionary biology, with wide range applications in related fields such as virology, epidemiology, clinic and public health, and antibiotic research. Significant progress has been made on this research since 1943 when Luria-Delbrück fluctuation analysis was first introduced. However, existing estimators of mutation rates are heavily reliant on model assumptions in fluctuation analysis, and become less applicable to real microbial experiments which deviate from the model assumptions. To overcome this difficulty, we propose to model fluctuation experimental data by a two-type Markov branching process (MBP) and use approximate Bayesian computation (ABC) to estimate the mutation probability parameters. Such an ABC-based mutation rate estimator is based on intensive simulations from the mutation process, thereby taking advantage of modern computing power. Most importantly, its likelihood-free feature allows more complex and realistic setups of the mutation process, especially when the distribution of the number of mutants cannot be easily derived. To further improve computation efficiency, we use a Gaussian process surrogate to substitute the simulator in the ABC algorithm, and call the resulting estimator GPS-ABC. Simulation studies show that, when used to estimate constant mutation rate in MBP, ABC-based estimators generally outperform traditional moment or likelihood-based estimators. When mutations occur in two stages, i.e., in MBP with a piece-wise constant mutation rate function, traditional mutation rate estimators become not applicable, yet GPS-ABC still achieves reasonable estimates. Finally, the proposed GPS-ABC estimator is used to analyze real fluctuation experimental datasets for studying drug resistance.

Population Genetics of Oncomelania hupensis Snails from New-Emerging Snail Habitats in a Currently Schistosoma japonicum Non-Endemic Area

Schistosomiasis is still one of the most significant neglected tropical diseases worldwide, and China is endemic for Schistosoma japonicum. With its great achievement in schistosomiasis control, the government of China has set the goal to eliminate the parasitic disease at the country level by 2030. However, one major challenge is the remaining huge areas of habitats for the intermediate host Oncomelania hupensis. This is further exacerbated by an increasing number of new emerging snail habitats reported each year. Therefore, population genetics on snails in such areas will be useful in evaluation of snail control effect and/or dispersal. We then sampled snails from new emerging habitats in Taicang of Jiangsu, China, a currently S. japonicum non-endemic area from 2014 to 2017, and performed population genetic analyses based on nine microsatellites. Results showed that all snail populations had low genetic diversity, and most genetic variations originated from within snail populations. The estimated effective population size for the 2015 population was infinitive. All snails could be separated into two clusters, and further DIYABC analysis revealed that both the 2016 and the 2017 populations may derive from the 2015, indicating that the 2017 population must have been missed in the field survey performed in 2016. These findings may have implications in development of more practical guidelines for snail monitoring and control.

Coalescent-based species delimitation meets deep learning: Insights from a highly fragmented cactus system

Delimiting species boundaries is a major goal in evolutionary biology. An increasing volume of literature has focused on the challenges of investigating cryptic diversity within complex evolutionary scenarios of speciation, including gene flow and demographic fluctuations. New methods based on model selection, such as approximate Bayesian computation, approximate likelihoods, and machine learning are promising tools arising in this field. Here, we introduce a framework for species delimitation using the multispecies coalescent model coupled with a deep learning algorithm based on convolutional neural networks (CNNs). We compared this strategy with a similar ABC approach. We applied both methods to test species boundary hypotheses based on current and previous taxonomic delimitations as well as genetic data (sequences from 41 loci) in Pilosocereus aurisetus, a cactus species complex with a sky-island distribution and taxonomic uncertainty. To validate our method, we also applied the same strategy on data from widely accepted species from the genus Drosophila. The results show that our CNN approach has a high capacity to distinguish among the simulated species delimitation scenarios, with higher accuracy than ABC. For the cactus data set, a splitter hypothesis without gene flow showed the highest probability in both CNN and ABC approaches, a result agreeing with previous taxonomic classifications and in line with the sky-island distribution and low dispersal of P. aurisetus. Our results highlight the cryptic diversity within the P. aurisetus complex and show that CNNs are a promising approach for distinguishing complex evolutionary histories, even outperforming the accuracy of other model-based approaches such as ABC.

ABCDP: Approximate Bayesian Computation with Differential Privacy

We developed a novel approximate Bayesian computation (ABC) framework, ABCDP, which produces differentially private (DP) and approximate posterior samples. Our framework takes advantage of the sparse vector technique (SVT), widely studied in the differential privacy literature. SVT incurs the privacy cost only when a condition (whether a quantity of interest is above/below a threshold) is met. If the condition is sparsely met during the repeated queries, SVT can drastically reduce the cumulative privacy loss, unlike the usual case where every query incurs the privacy loss. In ABC, the quantity of interest is the distance between observed and simulated data, and only when the distance is below a threshold can we take the corresponding prior sample as a posterior sample. Hence, applying SVT to ABC is an organic way to transform an ABC algorithm to a privacy-preserving variant with minimal modification, but yields the posterior samples with a high privacy level. We theoretically analyzed the interplay between the noise added for privacy and the accuracy of the posterior samples. We apply ABCDP to several data simulators and show the efficacy of the proposed framework.

The genomic history of the Aegean palatial civilizations

The Cycladic, the Minoan, and the Helladic (Mycenaean) cultures define the Bronze Age (BA) of Greece. Urbanism, complex social structures, craft and agricultural specialization, and the earliest forms of writing characterize this iconic period. We sequenced six Early to Middle BA whole genomes, along with 11 mitochondrial genomes, sampled from the three BA cultures of the Aegean Sea. The Early BA (EBA) genomes are homogeneous and derive most of their ancestry from Neolithic Aegeans, contrary to earlier hypotheses that the Neolithic-EBA cultural transition was due to massive population turnover. EBA Aegeans were shaped by relatively small-scale migration from East of the Aegean, as evidenced by the Caucasus-related ancestry also detected in Anatolians. In contrast, Middle BA (MBA) individuals of northern Greece differ from EBA populations in showing ∼50% Pontic-Caspian Steppe-related ancestry, dated at ca. 2,600-2,000 BCE. Such gene flow events during the MBA contributed toward shaping present-day Greek genomes.

Waves Out of the Korean Peninsula and Inter- and Intra-Species Replacements in Freshwater Fishes in Japan

The Japanese archipelago is located at the periphery of the continent of Asia. Rivers in the Japanese archipelago, separated from the continent of Asia by about 17 Ma, have experienced an intermittent exchange of freshwater fish taxa through a narrow land bridge generated by lowered sea level. As the Korean Peninsula and Japanese archipelago were not covered by an ice sheet during glacial periods, phylogeographical analyses in this region can trace the history of biota that were, for a long time, beyond the last glacial maximum. In this study, we analyzed the phylogeography of four freshwater fish taxa, Hemibarbus longirostris, dark chub Nipponocypris temminckii, Tanakia ssp. and Carassius ssp., whose distributions include both the Korean Peninsula and Western Japan. We found for each taxon that a small component of diverse Korean clades of freshwater fishes migrated in waves into the Japanese archipelago to form the current phylogeographic structure of biota. The replacements of indigenous populations by succeeding migrants may have also influenced the phylogeography.

The Early Peopling of the Philippines based on mtDNA

Despite the efforts made to reconstruct the history of modern humans, there are still poorly explored regions that are key for understanding the phylogeography of our species. One of them is the Philippines, which is crucial to unravel the colonization of Southeast Asia and Oceania but where little is known about when and how the first humans arrived. In order to shed light into this settlement, we collected samples from 157 individuals of the Philippines with the four grandparents belonging to the same region and mitochondrial variants older than 20,000 years. Next, we analyzed the hypervariable I mtDNA region by approximate Bayesian computation based on extensive spatially explicit computer simulations to select among several migration routes towards the Philippines and to estimate population genetic parameters of this colonization. We found that the colonization of the Philippines occurred more than 60,000 years ago, with long-distance dispersal and from both north and south migration routes. Our results also suggest an environmental scenario especially optimal for humans, with large carrying capacity and population growth, in comparison to other regions of Asia. In all, our study suggests a rapid expansion of modern humans towards the Philippines that could be associated with the establisment of maritime technologies and favorable environmental conditions.

Evolutionary Dynamics in Structured Populations Under Strong Population Genetic Forces

In the long-term neutral equilibrium, high rates of migration between subpopulations result in little population differentiation. However, in the short-term, even very abundant migration may not be enough for subpopulations to equilibrate immediately. In this study, we investigate dynamical patterns of short-term population differentiation in adapting populations via stochastic and analytical modeling through time. We characterize a regime in which selection and migration interact to create non-monotonic patterns of population differentiation over time when migration is weaker than selection, but stronger than drift. We demonstrate how these patterns can be leveraged to estimate high migration rates using approximate Bayesian computation. We apply this approach to estimate fast migration in a rapidly adapting intra-host Simian-HIV population sampled from different anatomical locations. We find differences in estimated migration rates between different compartments, even though all are above [Formula: see text] = 1. This work demonstrates how studying demographic processes on the timescale of selective sweeps illuminates processes too fast to leave signatures on neutral timescales.

Low genetic diversity in a critically endangered primate: shallow evolutionary history or recent population bottleneck?

BACKGROUND

Current patterns of population genetic variation may have been shaped by long-term evolutionary history and contemporary demographic processes. Understanding the underlying mechanisms that yield those patterns is crucial for informed conservation of endangered species. The critically endangered white-headed langur, Trachypithecus leucocephalus, is endemic to a narrow range in southwest China. This species shows very low genetic diversity in its 2 main relict populations, Fusui and Chongzuo. Whether this has been caused by a short evolutionary history or recent population declines is unknown. Therefore, we investigated the contributions of historical and recent population demographic changes to population genetic diversity by using 15 nuclear microsatellite markers and mitochondrial DNA (mtDNA) control region sequences.

RESULTS

Using genetic data from 214 individuals we found a total of 9 mtDNA haplotypes in the Fusui population but only 1 haplotype in the Chongzuo population, and we found an overall low genetic diversity (haplotype and nucleotide diversities: h = 0.486 ± 0.036; π = 0.0028 ± 0.0003). The demographic history inferred from mtDNA and microsatellite markers revealed no evidence for historical population size fluctuations or recent population bottlenecks. Simulations of possible population divergence histories inferred by DIYABC analysis supported a recent divergence of the Chongzuo population from the Fusui population and no population bottlenecks.

CONCLUSIONS

Despite severe population declines caused by anthropogenic activities in the last century, the low genetic diversity of the extant white-headed langur populations is most likely primarily due to the species' shallow evolutionary history and to a recent, local population founder event.

Speciation and hybridization in invasive fire ants

BACKGROUND

A major focus of evolutionary biology is the formation of reproductive barriers leading to divergence and ultimately, speciation. Often, it is not clear whether the separation of populations is complete or if there still is ongoing gene flow in the form of rare cases of admixture, known as isolation with migration. Here, we studied the speciation of two fire ant species, Solenopsis invicta and Solenopsis richteri, both native to South America, both inadvertently introduced to North America in the early twentieth century. While the two species are known to admix in the introduced range, in the native range no hybrids were found.

RESULTS

We conducted a population genomic survey of native and introduced populations of the two species using reduced representation genomic sequencing of 337 samples. Using maximum likelihood analysis over native range samples, we found no evidence of any gene flow between the species since they diverged. We estimated their time of divergence to 190,000 (100,000-350,000) generations ago. Modelling the demographic history of native and introduced S. invicta populations, we evaluated their divergence times and historic and contemporary population sizes, including the original founder population in North America, which was estimated at 26 (10-93) unrelated singly-mated queens.

CONCLUSIONS

We provide evidence for complete genetic isolation maintained between two invasive species in their natïve range, based, for the first time, on large scale genomic data analysis. The results lay the foundations for further studies into different stages in the formation of genetic barriers in dynamic, invasive populations.

Genomic evidence for shared common ancestry of East African hunting-gathering populations and insights into local adaptation

Anatomically modern humans arose in Africa ∼300,000 years ago, but the demographic and adaptive histories of African populations are not well-characterized. Here, we have generated a genome-wide dataset from 840 Africans, residing in western, eastern, southern, and northern Africa, belonging to 50 ethnicities, and speaking languages belonging to four language families. In addition to agriculturalists and pastoralists, our study includes 16 populations that practice, or until recently have practiced, a hunting-gathering (HG) lifestyle. We observe that genetic structure in Africa is broadly correlated not only with geography, but to a lesser extent, with linguistic affiliation and subsistence strategy. Four East African HG (EHG) populations that are geographically distant from each other show evidence of common ancestry: the Hadza and Sandawe in Tanzania, who speak languages with clicks classified as Khoisan; the Dahalo in Kenya, whose language has remnant clicks; and the Sabue in Ethiopia, who speak an unclassified language. Additionally, we observed common ancestry between central African rainforest HGs and southern African San, the latter of whom speak languages with clicks classified as Khoisan. With the exception of the EHG, central African rainforest HGs, and San, other HG groups in Africa appear genetically similar to neighboring agriculturalist or pastoralist populations. We additionally demonstrate that infectious disease, immune response, and diet have played important roles in the adaptive landscape of African history. However, while the broad biological processes involved in recent human adaptation in Africa are often consistent across populations, the specific loci affected by selective pressures more often vary across populations.

The impact of experimental design choices on parameter inference for models of growing cell colonies

To better understand development, repair and disease progression, it is useful to quantify the behaviour of proliferative and motile cell populations as they grow and expand to fill their local environment. Inferring parameters associated with mechanistic models of cell colony growth using quantitative data collected from carefully designed experiments provides a natural means to elucidate the relative contributions of various processes to the growth of the colony. In this work, we explore how experimental design impacts our ability to infer parameters for simple models of the growth of proliferative and motile cell populations. We adopt a Bayesian approach, which allows us to characterize the uncertainty associated with estimates of the model parameters. Our results suggest that experimental designs that incorporate initial spatial heterogeneities in cell positions facilitate parameter inference without the requirement of cell tracking, while designs that involve uniform initial placement of cells require cell tracking for accurate parameter inference. As cell tracking is an experimental bottleneck in many studies of this type, our recommendations for experimental design provide for significant potential time and cost savings in the analysis of cell colony growth.

The divergence history of European blue mussel species reconstructed from Approximate Bayesian Computation: the effects of sequencing techniques and sampling strategies

Genome-scale diversity data are increasingly available in a variety of biological systems, and can be used to reconstruct the past evolutionary history of species divergence. However, extracting the full demographic information from these data is not trivial, and requires inferential methods that account for the diversity of coalescent histories throughout the genome. Here, we evaluate the potential and limitations of one such approach. We reexamine a well-known system of mussel sister species, using the joint site frequency spectrum (jSFS) of synonymous mutations computed either from exome capture or RNA-seq, in an Approximate Bayesian Computation (ABC) framework. We first assess the best sampling strategy (number of: individuals, loci, and bins in the jSFS), and show that model selection is robust to variation in the number of individuals and loci. In contrast, different binning choices when summarizing the jSFS, strongly affect the results: including classes of low and high frequency shared polymorphisms can more effectively reveal recent migration events. We then take advantage of the flexibility of ABC to compare more realistic models of speciation, including variation in migration rates through time (i.e., periodic connectivity) and across genes (i.e., genome-wide heterogeneity in migration rates). We show that these models were consistently selected as the most probable, suggesting that mussels have experienced a complex history of gene flow during divergence and that the species boundary is semi-permeable. Our work provides a comprehensive evaluation of ABC demographic inference in mussels based on the coding jSFS, and supplies guidelines for employing different sequencing techniques and sampling strategies. We emphasize, perhaps surprisingly, that inferences are less limited by the volume of data, than by the way in which they are analyzed.

Two-scale dispersal estimation for biological invasions via synthetic likelihood

Biological invasions reshape environments and affect the ecological and economic welfare of states and communities. Such invasions advance on multiple spatial scales, complicating their control. When modeling stochastic dispersal processes, intractable likelihoods and autocorrelated data complicate parameter estimation. As with other approaches, the recent synthetic likelihood framework for stochastic models uses summary statistics to reduce this complexity; however, it additionally provides usable likelihoods, facilitating the use of existing likelihood-based machinery. Here, we extend this framework to parameterize multi-scale spatio-temporal dispersal models and compare existing and newly developed spatial summary statistics to characterize dispersal patterns. We provide general methods to evaluate potential summary statistics and present a fitting procedure that accurately estimates dispersal parameters on simulated data. Finally, we apply our methods to quantify the short and long range dispersal of Chagas disease vectors in urban Arequipa, Peru, and assess the feasibility of a purely reactive strategy to contain the invasion.

Supervised Machine Learning for Population Genetics: A New Paradigm

As population genomic datasets grow in size, researchers are faced with the daunting task of making sense of a flood of information. To keep pace with this explosion of data, computational methodologies for population genetic inference are rapidly being developed to best utilize genomic sequence data. In this review we discuss a new paradigm that has emerged in computational population genomics: that of supervised machine learning (ML). We review the fundamentals of ML, discuss recent applications of supervised ML to population genetics that outperform competing methods, and describe promising future directions in this area. Ultimately, we argue that supervised ML is an important and underutilized tool that has considerable potential for the world of evolutionary genomics.

Fundamentals and Recent Developments in Approximate Bayesian Computation

Bayesian inference plays an important role in phylogenetics, evolutionary biology, and in many other branches of science. It provides a principled framework for dealing with uncertainty and quantifying how it changes in the light of new evidence. For many complex models and inference problems, however, only approximate quantitative answers are obtainable. Approximate Bayesian computation (ABC) refers to a family of algorithms for approximate inference that makes a minimal set of assumptions by only requiring that sampling from a model is possible. We explain here the fundamentals of ABC, review the classical algorithms, and highlight recent developments. [ABC; approximate Bayesian computation; Bayesian inference; likelihood-free inference; phylogenetics; simulator-based models; stochastic simulation models; tree-based models.]

Inferring rooted species trees from unrooted gene trees using approximate Bayesian computation

Methods for inferring species trees from gene trees motivated by incomplete lineage sorting typically use either rooted gene trees to infer a rooted species tree, or use unrooted gene trees to infer an unrooted species tree, which is then typically rooted using one or more outgroups. Theoretically, however, it has been known since 2011 that it is possible to consistently infer the root of the species tree directly from unrooted gene trees without assuming an outgroup. Here, we use approximate Bayesian computation to infer the root of the species tree from unrooted gene trees assuming the multispecies coalescent model. It is hoped that this approach will be useful in cases where an appropriate outgroup is difficult to find and gene trees do not follow a molecular clock. We use approximate Bayesian computation to infer the root of the species tree from unrooted gene trees. This approach could also be useful when there is prior information that makes a small number of root locations plausible in an unrooted species tree.

Likelihood-free inference via classification

Increasingly complex generative models are being used across disciplines as they allow for realistic characterization of data, but a common difficulty with them is the prohibitively large computational cost to evaluate the likelihood function and thus to perform likelihood-based statistical inference. A likelihood-free inference framework has emerged where the parameters are identified by finding values that yield simulated data resembling the observed data. While widely applicable, a major difficulty in this framework is how to measure the discrepancy between the simulated and observed data. Transforming the original problem into a problem of classifying the data into simulated versus observed, we find that classification accuracy can be used to assess the discrepancy. The complete arsenal of classification methods becomes thereby available for inference of intractable generative models. We validate our approach using theory and simulations for both point estimation and Bayesian inference, and demonstrate its use on real data by inferring an individual-based epidemiological model for bacterial infections in child care centers.

Inferring epidemiological parameters from phylogenies using regression-ABC: A comparative study

Inferring epidemiological parameters such as the R₀ from time-scaled phylogenies is a timely challenge. Most current approaches rely on likelihood functions, which raise specific issues that range from computing these functions to finding their maxima numerically. Here, we present a new regression-based Approximate Bayesian Computation (ABC) approach, which we base on a large variety of summary statistics intended to capture the information contained in the phylogeny and its corresponding lineage-through-time plot. The regression step involves the Least Absolute Shrinkage and Selection Operator (LASSO) method, which is a robust machine learning technique. It allows us to readily deal with the large number of summary statistics, while avoiding resorting to Markov Chain Monte Carlo (MCMC) techniques. To compare our approach to existing ones, we simulated target trees under a variety of epidemiological models and settings, and inferred parameters of interest using the same priors. We found that, for large phylogenies, the accuracy of our regression-ABC is comparable to that of likelihood-based approaches involving birth-death processes implemented in BEAST2. Our approach even outperformed these when inferring the host population size with a Susceptible-Infected-Removed epidemiological model. It also clearly outperformed a recent kernel-ABC approach when assuming a Susceptible-Infected epidemiological model with two host types. Lastly, by re-analyzing data from the early stages of the recent Ebola epidemic in Sierra Leone, we showed that regression-ABC provides more realistic estimates for the duration parameters (latency and infectiousness) than the likelihood-based method. Overall, ABC based on a large variety of summary statistics and a regression method able to perform variable selection and avoid overfitting is a promising approach to analyze large phylogenies.

Given the rapid evolution of many pathogens, analysing their genomes by means of phylogenies can inform us about how they spread. This is the focus of the field known as “phylodynamics”. Most existing methods inferring epidemiological parameters from virus phylogenies are limited by the difficulty of handling complex likelihood functions, which commonly incorporate latent variables. Here, we use an alternative method known as regression-based Approximate Bayesian Computation (ABC), which circumvents this problem by using simulations and dataset comparisons. Since phylogenies are difficult to compare to one another, we introduce many summary statistics to describe them and take advantage of current machine learning techniques able to perform variable selection. We show that the accuracy we reach is comparable to that of existing methods. This accuracy increases with phylogeny size and can even be higher than that of existing methods for some parameters. Overall, regression-based ABC opens new perspectives to infer epidemiological parameters from large phylogenies.

Sequence diversity patterns suggesting balancing selection in partially sex-linked genes of the plant Silene latifolia are not generated by demographic history or gene flow

DNA sequence diversity in genes in the partially sex-linked pseudoautosomal region (PAR) of the sex chromosomes of the plant Silene latifolia is higher than expected from within-species diversity of other genes. This could be the footprint of sexually antagonistic (SA) alleles that are maintained by balancing selection in a PAR gene (or genes) and affect polymorphism in linked genome regions. SA selection is predicted to occur during sex chromosome evolution, but it is important to test whether the unexpectedly high sequence polymorphism could be explained without it, purely by the combined effects of partial linkage with the sex-determining region and the population's demographic history, including possible introgression from Silene dioica. To test this, we applied approximate Bayesian computation-based model choice to autosomal sequence diversity data, to find the most plausible scenario for the recent history of S. latifolia and then to estimate the posterior density of the most relevant parameters. We then used these densities to simulate variation to be expected at PAR genes. We conclude that an excess of variants at high frequencies at PAR genes should arise in S. latifolia populations only for genes with strong associations with fully sex-linked genes, which requires closer linkage with the fully sex-linked region than that estimated for the PAR genes where apparent deviations from neutrality were observed. These results support the need to invoke selection to explain the S. latifolia PAR gene diversity, and encourage further work to test the possibility of balancing selection due to sexual antagonism.

Importance of incomplete lineage sorting and introgression in the origin of shared genetic variation between two closely related pines with overlapping distributions

Genetic variation shared between closely related species may be due to retention of ancestral polymorphisms because of incomplete lineage sorting (ILS) and/or introgression following secondary contact. It is challenging to distinguish ILS and introgression because they generate similar patterns of shared genetic diversity, but this is nonetheless essential for inferring accurately the history of species with overlapping distributions. To address this issue, we sequenced 33 independent intron loci across the genome of two closely related pine species (Pinus massoniana Lamb. and Pinus hwangshanensis Hisa) from Southeast China. Population structure analyses revealed that the species showed slightly more admixture in parapatric populations than in allopatric populations. Levels of interspecific differentiation were lower in parapatry than in allopatry. Approximate Bayesian computation suggested that the most likely speciation scenario explaining this pattern was a long period of isolation followed by a secondary contact. Ecological niche modeling suggested that a gradual range expansion of P. hwangshanensis during the Pleistocene climatic oscillations could have been the cause of the overlap. Our study therefore suggests that secondary introgression, rather than ILS, explains most of the shared nuclear genomic variation between these two species and demonstrates the complementarity of population genetics and ecological niche modeling in understanding gene flow history. Finally, we discuss the importance of contrasting results from markers with different dynamics of migration, namely nuclear, chloroplast and mitochondrial DNA.

Likelihood-Free Inference in High-Dimensional Models

Methods that bypass analytical evaluations of the likelihood function have become an indispensable tool for statistical inference in many fields of science. These so-called likelihood-free methods rely on accepting and rejecting simulations based on summary statistics, which limits them to low-dimensional models for which the value of the likelihood is large enough to result in manageable acceptance rates. To get around these issues, we introduce a novel, likelihood-free Markov chain Monte Carlo (MCMC) method combining two key innovations: updating only one parameter per iteration and accepting or rejecting this update based on subsets of statistics approximately sufficient for this parameter. This increases acceptance rates dramatically, rendering this approach suitable even for models of very high dimensionality. We further derive that for linear models, a one-dimensional combination of statistics per parameter is sufficient and can be found empirically with simulations. Finally, we demonstrate that our method readily scales to models of very high dimensionality, using toy models as well as by jointly inferring the effective population size, the distribution of fitness effects (DFE) of segregating mutations, and selection coefficients for each locus from data of a recent experiment on the evolution of drug resistance in influenza.

Deep Learning for Population Genetic Inference

Given genomic variation data from multiple individuals, computing the likelihood of complex population genetic models is often infeasible. To circumvent this problem, we introduce a novel likelihood-free inference framework by applying deep learning, a powerful modern technique in machine learning. Deep learning makes use of multilayer neural networks to learn a feature-based function from the input (e.g., hundreds of correlated summary statistics of data) to the output (e.g., population genetic parameters of interest). We demonstrate that deep learning can be effectively employed for population genetic inference and learning informative features of data. As a concrete application, we focus on the challenging problem of jointly inferring natural selection and demography (in the form of a population size change history). Our method is able to separate the global nature of demography from the local nature of selection, without sequential steps for these two factors. Studying demography and selection jointly is motivated by Drosophila, where pervasive selection confounds demographic analysis. We apply our method to 197 African Drosophila melanogaster genomes from Zambia to infer both their overall demography, and regions of their genome under selection. We find many regions of the genome that have experienced hard sweeps, and fewer under selection on standing variation (soft sweep) or balancing selection. Interestingly, we find that soft sweeps and balancing selection occur more frequently closer to the centromere of each chromosome. In addition, our demographic inference suggests that previously estimated bottlenecks for African Drosophila melanogaster are too extreme.

Genomic Signatures of Selective Pressures and Introgression from Archaic Hominins at Human Innate Immunity Genes

Human genes governing innate immunity provide a valuable tool for the study of the selective pressure imposed by microorganisms on host genomes. A comprehensive, genome-wide study of how selective constraints and adaptations have driven the evolution of innate immunity genes is missing. Using full-genome sequence variation from the 1000 Genomes Project, we first show that innate immunity genes have globally evolved under stronger purifying selection than the remainder of protein-coding genes. We identify a gene set under the strongest selective constraints, mutations in which are likely to predispose individuals to life-threatening disease, as illustrated by STAT1 and TRAF3. We then evaluate the occurrence of local adaptation and detect 57 high-scoring signals of positive selection at innate immunity genes, variation in which has been associated with susceptibility to common infectious or autoimmune diseases. Furthermore, we show that most adaptations targeting coding variation have occurred in the last 6,000-13,000 years, the period at which populations shifted from hunting and gathering to farming. Finally, we show that innate immunity genes present higher Neandertal introgression than the remainder of the coding genome. Notably, among the genes presenting the highest Neandertal ancestry, we find the TLR6-TLR1-TLR10 cluster, which also contains functional adaptive variation in Europeans. This study identifies highly constrained genes that fulfill essential, non-redundant functions in host survival and reveals others that are more permissive to change-containing variation acquired from archaic hominins or adaptive variants in specific populations-improving our understanding of the relative biological importance of innate immunity pathways in natural conditions.

Reliable ABC model choice via random forests

MOTIVATION

Approximate Bayesian computation (ABC) methods provide an elaborate approach to Bayesian inference on complex models, including model choice. Both theoretical arguments and simulation experiments indicate, however, that model posterior probabilities may be poorly evaluated by standard ABC techniques.

RESULTS

We propose a novel approach based on a machine learning tool named random forests (RF) to conduct selection among the highly complex models covered by ABC algorithms. We thus modify the way Bayesian model selection is both understood and operated, in that we rephrase the inferential goal as a classification problem, first predicting the model that best fits the data with RF and postponing the approximation of the posterior probability of the selected model for a second stage also relying on RF. Compared with earlier implementations of ABC model choice, the ABC RF approach offers several potential improvements: (i) it often has a larger discriminative power among the competing models, (ii) it is more robust against the number and choice of statistics summarizing the data, (iii) the computing effort is drastically reduced (with a gain in computation efficiency of at least 50) and (iv) it includes an approximation of the posterior probability of the selected model. The call to RF will undoubtedly extend the range of size of datasets and complexity of models that ABC can handle. We illustrate the power of this novel methodology by analyzing controlled experiments as well as genuine population genetics datasets.

AVAILABILITY AND IMPLEMENTATION

The proposed methodology is implemented in the R package abcrf available on the CRAN.

CONTACT

jean-michel.marin@umontpellier.fr

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Evolutionary processes driving spatial patterns of intraspecific genetic diversity in river ecosystems

Describing, understanding and predicting the spatial distribution of genetic diversity is a central issue in biological sciences. In river landscapes, it is generally predicted that neutral genetic diversity should increase downstream, but there have been few attempts to test and validate this assumption across taxonomic groups. Moreover, it is still unclear what are the evolutionary processes that may generate this apparent spatial pattern of diversity. Here, we quantitatively synthesized published results from diverse taxa living in river ecosystems, and we performed a meta-analysis to show that a downstream increase in intraspecific genetic diversity (DIGD) actually constitutes a general spatial pattern of biodiversity that is repeatable across taxa. We further demonstrated that DIGD was stronger for strictly waterborne dispersing than for overland dispersing species. However, for a restricted data set focusing on fishes, there was no evidence that DIGD was related to particular species traits. We then searched for general processes underlying DIGD by simulating genetic data in dendritic-like river systems. Simulations revealed that the three processes we considered (downstream-biased dispersal, increase in habitat availability downstream and upstream-directed colonization) might generate DIGD. Using random forest models, we identified from simulations a set of highly informative summary statistics allowing discriminating among the processes causing DIGD. Finally, combining these discriminant statistics and approximate Bayesian computations on a set of twelve empirical case studies, we hypothesized that DIGD were most likely due to the interaction of two of these three processes and that contrary to expectation, they were not solely caused by downstream-biased dispersal.

Using ABC and microsatellite data to detect multiple introductions of invasive species from a single source

The introduction of invasive species to new locations (that is, biological invasions) can have major impact on biodiversity, agriculture and public health. As such, determining the routes and modality of introductions with genetic data has become a fundamental goal in molecular ecology. To assist with this goal, new statistical methods and frameworks have been developed, such as approximate Bayesian computation (ABC) for inferring invasion history. Here, we present a model of invasion accounting for multiple introductions from a single source (MISS), a heretofore largely unexplored model. We simulate microsatellite data to evaluate the power of ABC to distinguish between single and multiple introductions from the same source, under a range of demographic parameters. We also apply ABC to microsatellite data from three invasions of bumblebee in New Zealand. In addition, we assess the performance of several methods of summary statistics selection. Our simulated results suggested good ability to distinguish between one- and two-wave models over much but not all of the parameter space tested, independent of summary statistics used. Globally, parameter estimation was good except for bottleneck timing. For one of the bumblebee species, we clearly rejected the MISS model, while for the other two we found inconclusive results. Since a second wave may provide genetic reinforcement to initial colonists, help relieve inbreeding among founders, or increase the hazard of the invasion, its detection may be crucial for managing invasions; we suggest that the MISS model could be considered as a potential model in future theoretical and empirical studies of invasions.

Modeling tuberculosis dynamics, detection and control in cattle herds

Epidemiological models are key tools for designing and evaluating detection and control strategies against animal infectious diseases. In France, after decades of decrease of bovine tuberculosis (bTB) incidence, the disease keeps circulating. Increasing prevalence levels are observed in several areas, where the detection and control strategy could be adapted. The objective of this work was to design and calibrate a model of the within-herd transmission of bTB. The proposed model is a stochastic model operating in discrete-time. Three health states were distinguished: susceptible, latent and infected. Dairy and beef herd dynamics and bTB detection and control programs were explicitly represented. Approximate Bayesian computation was used to estimate three model parameters from field data: the transmission parameter when animals are inside (β_inside) and outside (β_outside) buildings, and the duration of the latent phase. An independent dataset was used for model validation. The estimated median was 0.43 [0.16–0.84] month⁻¹ for β_inside and 0.08 [0.01–0.32] month⁻¹ for β_outside. The median duration of the latent period was estimated 3.5 [2]–[8] months. The sensitivity analysis showed only minor influences of fixed parameter values on these posterior estimates. Validation based on an independent dataset showed that in more than 80% of herds, the observed proportion of animals with detected lesions was between the 2.5% and 97.5% percentiles of the simulated distribution. In the absence of control program and once bTB has become enzootic within a herd, the median effective reproductive ratio was estimated to be 2.2 in beef herds and 1.7 in dairy herds. These low estimates are consistent with field observations of a low prevalence level in French bTB-infected herds.

ABC inference of multi-population divergence with admixture from unphased population genomic data

Rapidly developing sequencing technologies and declining costs have made it possible to collect genome-scale data from population-level samples in nonmodel systems. Inferential tools for historical demography given these data sets are, at present, underdeveloped. In particular, approximate Bayesian computation (ABC) has yet to be widely embraced by researchers generating these data. Here, we demonstrate the promise of ABC for analysis of the large data sets that are now attainable from nonmodel taxa through current genomic sequencing technologies. We develop and test an ABC framework for model selection and parameter estimation, given histories of three-population divergence with admixture. We then explore different sampling regimes to illustrate how sampling more loci, longer loci or more individuals affects the quality of model selection and parameter estimation in this ABC framework. Our results show that inferences improved substantially with increases in the number and/or length of sequenced loci, while less benefit was gained by sampling large numbers of individuals. Optimal sampling strategies given our inferential models included at least 2000 loci, each approximately 2 kb in length, sampled from five diploid individuals per population, although specific strategies are model and question dependent. We tested our ABC approach through simulation-based cross-validations and illustrate its application using previously analysed data from the oak gall wasp, Biorhiza pallida.

Matrix inversions for chromosomal inversions: a method to construct summary statistics in complex coalescent models

Chromosomal inversions allow genetic divergence of locally adapted populations by reducing recombination between chromosomes with different arrangements. While patterns of genetic variation within inverted regions are increasingly documented, inferential methods are largely missing to analyze such data. Previous work has provided expectations for coalescence patterns of neutral sites linked to an inversion polymorphism in two locally adapted populations. Here, we define a method to construct summary statistics in such complex population structure models. Under a scenario of selection on the inversion breakpoints, we first construct estimators of the migration rate between the two habitats, and of the recombination rate of a nucleotide site between the two inversion backgrounds. Next, we analyze the disequilibrium between two sites within an inversion and provide an estimator of the distinct recombination rate between these two sites in homokaryotypes and heterokaryotypes. These estimators should be suitable summary statistics for simulation-based methods that can handle the complex dependences in the data.

Approximate Bayesian inference for complex ecosystems

Mathematical models have been central to ecology for nearly a century. Simple models of population dynamics have allowed us to understand fundamental aspects underlying the dynamics and stability of ecological systems. What has remained a challenge, however, is to meaningfully interpret experimental or observational data in light of mathematical models. Here, we review recent developments, notably in the growing field of approximate Bayesian computation (ABC), that allow us to calibrate mathematical models against available data. Estimating the population demographic parameters from data remains a formidable statistical challenge. Here, we attempt to give a flavor and overview of ABC and its applications in population biology and ecology and eschew a detailed technical discussion in favor of a general discussion of the advantages and potential pitfalls this framework offers to population biologists.

Amount of information needed for model choice in Approximate Bayesian Computation

Approximate Bayesian Computation (ABC) has become a popular technique in evolutionary genetics for elucidating population structure and history due to its flexibility. The statistical inference framework has benefited from significant progress in recent years. In population genetics, however, its outcome depends heavily on the amount of information in the dataset, whether that be the level of genetic variation or the number of samples and loci. Here we look at the power to reject a simple constant population size coalescent model in favor of a bottleneck model in datasets of varying quality. Not only is this power dependent on the number of samples and loci, but it also depends strongly on the level of nucleotide diversity in the observed dataset. Whilst overall model choice in an ABC setting is fairly powerful and quite conservative with regard to false positives, detecting weaker bottlenecks is problematic in smaller or less genetically diverse datasets and limits the inferences possible in non-model organism where the amount of information regarding the two models is often limited. Our results show it is important to consider these limitations when performing an ABC analysis and that studies should perform simulations based on the size and nature of the dataset in order to fully assess the power of the study.

Human-facilitated metapopulation dynamics in an emerging pest species, Cimex lectularius

The number and demographic history of colonists can have dramatic consequences for the way in which genetic diversity is distributed and maintained in a metapopulation. The bed bug (Cimex lectularius) is a re-emerging pest species whose close association with humans has led to frequent local extinction and colonization, that is, to metapopulation dynamics. Pest control limits the lifespan of subpopulations, causing frequent local extinctions, and human-facilitated dispersal allows the colonization of empty patches. Founder events often result in drastic reductions in diversity and an increased influence of genetic drift. Coupled with restricted migration, this can lead to rapid population differentiation. We therefore predicted strong population structuring. Here, using 21 newly characterized microsatellite markers and approximate Bayesian computation (ABC), we investigate simplified versions of two classical models of metapopulation dynamics, in a coalescent framework, to estimate the number and genetic composition of founders in the common bed bug. We found very limited diversity within infestations but high degrees of structuring across the city of London, with extreme levels of genetic differentiation between infestations (F_ST = 0.59). ABC results suggest a common origin of all founders of a given subpopulation and that the numbers of colonists were low, implying that even a single mated female is enough to found a new infestation successfully. These patterns of colonization are close to the predictions of the propagule pool model, where all founders originate from the same parental infestation. These results show that aspects of metapopulation dynamics can be captured in simple models and provide insights that are valuable for the future targeted control of bed bug infestations.

Origin and demographic history of the endemic Taiwan spruce (Picea morrisonicola)

Taiwan spruce (Picea morrisonicola) is a vulnerable conifer species endemic to the island of Taiwan. A warming climate and competition from subtropical tree species has limited the range of Taiwan spruce to the higher altitudes of the island. Using seeds sampled from an area in the central mountain range of Taiwan, 15 nuclear loci were sequenced in order to measure genetic variation and to assess the long-term genetic stability of the species. Genetic diversity is low and comparable to other spruce species with limited ranges such as Picea breweriana, Picea chihuahuana, and Picea schrenkiana. Importantly, analysis using approximate Bayesian computation (ABC) provides evidence for a drastic decline in the effective population size approximately 0.3–0.5 million years ago (mya). We used simulations to show that this is unlikely to be a false-positive result due to the limited sample used here. To investigate the phylogenetic origin of Taiwan spruce, additional sequencing was performed in the Chinese spruce Picea wilsonii and combined with previously published data for three other mainland China species, Picea purpurea, Picea likiangensis, and P. schrenkiana. Analysis of population structure revealed that P. morrisonicola clusters most closely with P. wilsonii, and coalescent analyses using the program MIMAR dated the split to 4–8 mya, coincidental to the formation of Taiwan. Considering the population decrease that occurred after the split, however, led to a much more recent origin.

Robust demographic inference from genomic and SNP data

We introduce a flexible and robust simulation-based framework to infer demographic parameters from the site frequency spectrum (SFS) computed on large genomic datasets. We show that our composite-likelihood approach allows one to study evolutionary models of arbitrary complexity, which cannot be tackled by other current likelihood-based methods. For simple scenarios, our approach compares favorably in terms of accuracy and speed with ∂a∂i, the current reference in the field, while showing better convergence properties for complex models. We first apply our methodology to non-coding genomic SNP data from four human populations. To infer their demographic history, we compare neutral evolutionary models of increasing complexity, including unsampled populations. We further show the versatility of our framework by extending it to the inference of demographic parameters from SNP chips with known ascertainment, such as that recently released by Affymetrix to study human origins. Whereas previous ways of handling ascertained SNPs were either restricted to a single population or only allowed the inference of divergence time between a pair of populations, our framework can correctly infer parameters of more complex models including the divergence of several populations, bottlenecks and migration. We apply this approach to the reconstruction of African demography using two distinct ascertained human SNP panels studied under two evolutionary models. The two SNP panels lead to globally very similar estimates and confidence intervals, and suggest an ancient divergence (>110 Ky) between Yoruba and San populations. Our methodology appears well suited to the study of complex scenarios from large genomic data sets.

Inferring the geographic mode of speciation by contrasting autosomal and sex-linked genetic diversity

When geographic isolation drives speciation, concurrent termination of gene flow among genomic regions will occur immediately after the formation of the barrier between diverging populations. Alternatively, if speciation is driven by ecologically divergent selection, gene flow of selectively neutral genomic regions may go on between diverging populations until the completion of reproductive isolation. It may also lead to an unsynchronized termination of gene flow between genomic regions with different roles in the speciation process. Here, we developed a novel Approximate Bayesian Computation pipeline to infer the geographic mode of speciation by testing for a lack of postdivergence gene flow and a concurrent termination of gene flow in autosomal and sex-linked markers jointly. We applied this approach to infer the geographic mode of speciation for two allopatric highland rosefinches, the vinaceous rosefinch Carpodacus vinaceus and the Taiwan rosefinch C. formosanus from DNA polymorphisms of both autosomal and Z-linked loci. Our results suggest that the two rosefinch species diverged allopatrically approximately 0.5 Ma. Our approach allowed us further to infer that female effective population sizes are about five times larger than those of males, an estimate potentially useful when comparing the intensity of sexual selection across species.