Population Genetic Simulations Research Articles

Gene flow biases population genetic inference of recombination rate.

Accurate estimates of the rate of recombination are key to understanding a host of evolutionary processes as well as the evolution of the recombination rate itself. Model-based population genetic methods that infer recombination rates from patterns of linkage disequilibrium in the genome have become a popular method to estimate rates of recombination. However, these linkage disequilibrium-based methods make a variety of simplifying assumptions about the populations of interest that are often not met in natural populations. One such assumption is the absence of gene flow from other populations. Here, we use forward-time population genetic simulations of isolation-with-migration scenarios to explore how gene flow affects the accuracy of linkage disequilibrium-based estimators of recombination rate. We find that moderate levels of gene flow can result in either the overestimation or underestimation of recombination rates by up to 20-50% depending on the timing of divergence. We also find that these biases can affect the detection of interpopulation differences in recombination rate, causing both false positives and false negatives depending on the scenario. We discuss future possibilities for mitigating these biases and recommend that investigators exercise caution and confirm that their study populations meet assumptions before deploying these methods.

Empirical estimates of the mutation rate for an alphabaculovirus.

Mutation rates are of key importance for understanding evolutionary processes and predicting their outcomes. Empirical mutation rate estimates are available for a number of RNA viruses, but few are available for DNA viruses, which tend to have larger genomes. Whilst some viruses have very high mutation rates, lower mutation rates are expected for viruses with large genomes to ensure genome integrity. Alphabaculoviruses are insect viruses with large genomes and often have high levels of polymorphism, suggesting high mutation rates despite evidence of proofreading activity by the replication machinery. Here, we report an empirical estimate of the mutation rate per base per strand copying (s/n/r) of Autographa californica multiple nucleopolyhedrovirus (AcMNPV). To avoid biases due to selection, we analyzed mutations that occurred in a stable, non-functional genomic insert after five serial passages in Spodoptera exigua larvae. Our results highlight that viral demography and the stringency of mutation calling affect mutation rate estimates, and that using a population genetic simulation model to make inferences can mitigate the impact of these processes on estimates of mutation rate. We estimated a mutation rate of μ = 1×10-7 s/n/r when applying the most stringent criteria for mutation calling, and estimates of up to μ = 5×10-7 s/n/r when relaxing these criteria. The rates at which different classes of mutations accumulate provide good evidence for neutrality of mutations occurring within the inserted region. We therefore present a robust approach for mutation rate estimation for viruses with stable genomes, and strong evidence of a much lower alphabaculovirus mutation rate than supposed based on the high levels of polymorphism observed.

RAD-Seq and Ecological Niche Reveal Genetic Diversity, Phylogeny, and Geographic Distribution of Kadsura interior and Its Closely Related Species.

Most plants of Kadsura have economic value and medicinal application. Among them, K. interior and its closely related species have been demonstrated to have definite efficacy. However, the taxonomy and phylogenetic relationship of Kadsura in terms of morphology and commonly used gene regions remain controversial, which adversely affects its rational application. In this study, a total of 107 individuals of K. interior, K. heteroclita, K. longipedunculata, K. oblongifolia, and K. coccinea were studied from the perspectives of genetic diversity, phylogeny, and ecology via single nucleotide polymorphisms (SNPs) developed through restriction site-associated DNA sequencing (RAD-seq). Based on these SNPs, the genetic diversity, phylogenetic reconstruction, and population genetic structure were analyzed. Subsequently, divergence time estimation and differentiation scenario simulation were performed. Meanwhile, according to the species distribution records and bioclimatic variables, the Last Glacial Maximum and current potential distributions of five species were constructed, and the main ecological factors affecting the distribution of different species were extracted. The FST calculated showed that there was a moderate degree of differentiation among K. heteroclita, K. longipedunculata, and K. oblongifolia, and there was a high degree of genetic differentiation between K. interior and the above species. The phylogenetic tree indicated that each of the species was monophyletic. The results of population genetic structure and divergence scenario simulation and D-statistics showed that there were admixture and gene flow among K. heteroclita, K. longipedunculata, and K. oblongifolia. The results of ecological niche modeling indicated that the distribution areas and the bioclimatic variables affecting the distribution of K. interior and its related species were different. This study explored the differences in the genetic divergence and geographical distribution patterns of K. interior and its related species, clarifying the uniqueness of K. interior compared to its relatives and providing a reference for their rational application in the future.

Dominance Can Increase Genetic Variance After a Population Bottleneck: A Synthesis of the Theoretical and Empirical Evidence.

Drastic reductions in population size, or population bottlenecks, can lead to a reduction in additive genetic variance and adaptive potential. Genetic variance for some quantitative genetic traits, however, can increase after a population reduction. Empirical evaluations of quantitative traits following experimental bottlenecks indicate that non-additive genetic effects, including both allelic dominance at a given locus and epistatic interactions among loci, may impact the additive variance contributed by alleles that ultimately influences phenotypic expression and fitness. The dramatic effects of bottlenecks on overall genetic diversity have been well studied, but relatively little is known about how dominance and demographic events like bottlenecks can impact additive genetic variance. Herein, we critically examine how the degree of dominance among alleles affects additive genetic variance after a bottleneck. We first review and synthesize studies that document the impact of empirical bottlenecks on dominance variance. We then extend earlier work by elaborating on 2 theoretical models that illustrate the relationship between dominance and the potential increase in additive genetic variance immediately following a bottleneck. Furthermore, we investigate the parameters that influence the maximum level of genetic variation (associated with adaptive potential) after a bottleneck, including the number of founding individuals. Finally, we validated our methods using forward-time population genetic simulations of loci with varying dominance and selection levels. The fate of non-additive genetic variation following bottlenecks could have important implications for conservation and management efforts in a wide variety of taxa, and our work should help contextualize future studies (e.g., epistatic variance) in population genomics.

Overcoming constraints on the detection of recessive selection in human genes from population frequency data

The identification of genes that evolve under recessive natural selection is a long-standing goal of population genetics research that has important applications to the discovery of genes associated with disease. We found that commonly used methods to evaluate selective constraint at the gene level are highly sensitive to genes under heterozygous selection but ubiquitously fail to detect recessively evolving genes. Additionally, more sophisticated likelihood-based methods designed to detect recessivity similarly lack power for a human gene of realistic length from current population sample sizes. However, extensive simulations suggested that recessive genes may be detectable in aggregate. Here, we offer a method informed by population genetics simulations designed to detect recessive purifying selection in gene sets. Applying this to empirical gene sets produced significant enrichments for strong recessive selection in genes previously inferred to be under recessive selection in a consanguineous cohort and in genes involved in autosomal recessive monogenic disorders.

Gene network simulations provide testable predictions for the molecular domestication syndrome.

The domestication of plant species leads to repeatable morphological evolution, often referred to as the phenotypic domestication syndrome. Domestication is also associated with important genomic changes, such as the loss of genetic diversity compared with adequately large wild populations, and modifications of gene expression patterns. Here, we explored theoretically the effect of a domestication-like scenario on the evolution of gene regulatory networks. We ran population genetics simulations in which individuals were featured by their genotype (an interaction matrix encoding a gene regulatory network) and their gene expressions, representing the phenotypic level. Our domestication scenario included a population bottleneck and a selection switch mimicking human-mediated directional and canalizing selection, i.e., change in the optimal gene expression level and selection toward more stable expression across environments. We showed that domestication profoundly alters genetic architectures. Based on four examples of plant domestication scenarios, our simulations predict (1) a drop in neutral allelic diversity; (2) a change in gene expression variance that depends upon the domestication scenario; (3) transient maladaptive plasticity; (4) a deep rewiring of the gene regulatory networks, with a trend toward gain of regulatory interactions; and (5) a global increase in the genetic correlations among gene expressions, with a loss of modularity in the resulting coexpression patterns and in the underlying networks. We provide empirically testable predictions on the differences of genetic architectures between wild and domesticated forms. The characterization of such systematic evolutionary changes in the genetic architecture of traits contributes to define a molecular domestication syndrome.

The maintenance of polygenic sex determination depends on the dominance of fitness effects which are predictive of the role of sexual antagonism.

In species with polygenic sex determination (PSD), multiple male- and female-determining loci on different proto-sex chromosomes segregate as polymorphisms within populations. The extent to which these polymorphisms are at stable equilibria is not yet resolved. Previous work demonstrated that PSD is most likely to be maintained as a stable polymorphism when the proto-sex chromosomes have opposite (sexually antagonistic) fitness effects in males and females. However, these models usually consider PSD systems with only two proto-sex chromosomes, or they do not broadly consider the dominance of the alleles under selection. To address these shortcomings, I used forward population genetic simulations to identify selection pressures that can maintain PSD under different dominance scenarios in a system with more than two proto-sex chromosomes (modeled after the house fly). I found that overdominant fitness effects of male-determining proto-Y chromosomes are more likely to maintain PSD than dominant, recessive, or additive fitness effects. The overdominant fitness effects that maintain PSD tend to have proto-Y chromosomes with sexually antagonistic effects (male-beneficial and female-detrimental). In contrast, dominant fitness effects that maintain PSD tend to have sexually antagonistic multi-chromosomal genotypes, but the individual proto-sex chromosomes do not have sexually antagonistic effects. These results demonstrate that sexual antagonism can be an emergent property of the multi-chromosome genotype without individual sexually antagonistic chromosomes. My results further illustrate how the dominance of fitness effects has consequences for both the likelihood that PSD will be maintained as well as the role sexually antagonistic selection is expected to play in maintaining the polymorphism.

Strongly deleterious mutations are a primary determinant of extinction risk due to inbreeding depression.

Human‐driven habitat fragmentation and loss have led to a proliferation of small and isolated plant and animal populations with high risk of extinction. One of the main threats to extinction in these populations is inbreeding depression, which is primarily caused by recessive deleterious mutations becoming homozygous due to inbreeding. The typical approach for managing these populations is to maintain high genetic diversity, increasingly by translocating individuals from large populations to initiate a “genetic rescue.” However, the limitations of this approach have recently been highlighted by the demise of the gray wolf population on Isle Royale, which declined to the brink of extinction soon after the arrival of a migrant from the large mainland wolf population. Here, we use a novel population genetic simulation framework to investigate the role of genetic diversity, deleterious variation, and demographic history in mediating extinction risk due to inbreeding depression in small populations. We show that, under realistic models of dominance, large populations harbor high levels of recessive strongly deleterious variation due to these mutations being hidden from selection in the heterozygous state. As a result, when large populations contract, they experience a substantially elevated risk of extinction after these strongly deleterious mutations are exposed by inbreeding. Moreover, we demonstrate that, although genetic rescue is broadly effective as a means to reduce extinction risk, its effectiveness can be greatly increased by drawing migrants from small or moderate‐sized source populations rather than large source populations due to smaller populations harboring lower levels of recessive strongly deleterious variation. Our findings challenge the traditional conservation paradigm that focuses on maximizing genetic diversity in small populations in favor of a view that emphasizes minimizing strongly deleterious variation. These insights have important implications for managing small and isolated populations in the increasingly fragmented landscape of the Anthropocene.

Global adaptation complicates the interpretation of genome scans for local adaptation.

Spatially varying selection promotes variance in allele frequencies, increasing genetic differentiation between the demes of a metapopulation. For that reason, outliers in the genome‐wide distribution of summary statistics measuring genetic differentiation, such as FST, are often interpreted as evidence for alleles that contribute to local adaptation. However, theoretical studies have shown that in spatially structured populations the spread of beneficial mutations with spatially uniform fitness effects can also induce transient genetic differentiation. In recent years, numerous empirical studies have suggested that such species‐wide, or global, adaptation makes a substantial contribution to molecular evolution. In this perspective, we discuss how commonly such global adaptation may influence the genome‐wide distribution of FST and generate genetic differentiation patterns, which could be mistaken for local adaptation. To illustrate this, we use forward‐in‐time population genetic simulations assuming parameters for the rate and strength of beneficial mutations consistent with estimates from natural populations. We demonstrate that the spread of globally beneficial mutations in parapatric populations may frequently generate FST outliers, which could be misinterpreted as evidence for local adaptation. The spread of beneficial mutations causes selective sweeps at flanking sites, so in some cases, the effects of global versus local adaptation may be distinguished by examining patterns of nucleotide diversity within and between populations in addition to FST. However, when local adaptation has been only recently established, it may be much more difficult to distinguish from global adaptation, due to less accumulation of linkage disequilibrium at flanking sites. Through our discussion, we conclude that a large fraction of FST outliers that are presumed to arise from local adaptation may instead be due to global adaptation.

Adaptive gene misregulation.

Because gene expression is important for evolutionary adaptation, its misregulation is an important cause of maladaptation. A misregulated gene can be incorrectly silent ("off") when a transcription factor (TF) that is required for its activation does not binds its regulatory region. Conversely, a misregulated gene can be incorrectly active ("on") when a TF not normally involved in its activation binds its regulatory region, a phenomenon also known as regulatory crosstalk. DNA mutations that destroy or create TF binding sites on DNA are an important source of misregulation and crosstalk. Although misregulation reduces fitness in an environment to which an organism is well-adapted, it may become adaptive in a new environment. Here, I derive simple yet general mathematical expressions that delimit the conditions under which misregulation can be adaptive. These expressions depend on the strength of selection against misregulation, on the fraction of DNA sequence space filled with TF binding sites, and on the fraction of genes that must be expressed for optimal adaptation. I then use empirical data from RNA sequencing, protein-binding microarrays, and genome evolution, together with population genetic simulations to ask when these conditions are likely to be met. I show that they can be met under realistic circumstances, but these circumstances may vary among organisms and environments. My analysis provides a framework in which improved theory and data collection can help us demonstrate the role of misregulation in adaptation. It also shows that misregulation, like DNA mutation, is one of life's many imperfections that can help propel Darwinian evolution.

The effect of gene flow from unsampled demes in landscape genetic analysis.

An assumption of correlative landscape genetic methods is that genetic differentiation at neutral markers arises solely from the degree to which the intervening landscape between individuals or populations resists gene flow. However, this assumption is violated when gene flow occurs into the sampled population from an unsampled, differentiated deme. This may happen when sampling within only a portion of a population's extent or when closely related species hybridize with the sampled population. In both cases, violation of the modelling assumptions has the potential to reduce landscape genetic model selection accuracy and result in poor inferences. We used individual-based population genetic simulations in complex landscapes within a model selection framework to explore the potential confounding effect of gene flow from unsampled demes. We hypothesized that as gene flow from outside the sampling extent increased, model selection accuracy would decrease due to the formation of a hybrid zone where allele frequencies were perturbed in a way that was not correlated with effective distances between sampled individuals. Surprisingly, we found this expectation was unfounded, because the reduced accuracy due to admixture was counteracted by an increase in allelic diversity as alleles spread from the unsampled deme into the sampled population. These new alleles increased the power to detect landscape genetic relationships and even slightly improving model selection accuracy overall. This is a reassuring result, suggesting that sampling the full extent of a population or related species that may hybridize may be unnecessary, as long as other well-established sampling requirements are met.

Standardizing population genetics simulations.

Stdpopsim facilitates realistic population genetics simulations with a catalog of widely used models.

A community-maintained standard library of population genetic models.

The explosion in population genomic data demands ever more complex modes of analysis, and increasingly, these analyses depend on sophisticated simulations. Recent advances in population genetic simulation have made it possible to simulate large and complex models, but specifying such models for a particular simulation engine remains a difficult and error-prone task. Computational genetics researchers currently re-implement simulation models independently, leading to inconsistency and duplication of effort. This situation presents a major barrier to empirical researchers seeking to use simulations for power analyses of upcoming studies or sanity checks on existing genomic data. Population genetics, as a field, also lacks standard benchmarks by which new tools for inference might be measured. Here, we describe a new resource, stdpopsim, that attempts to rectify this situation. Stdpopsim is a community-driven open source project, which provides easy access to a growing catalog of published simulation models from a range of organisms and supports multiple simulation engine backends. This resource is available as a well-documented python library with a simple command-line interface. We share some examples demonstrating how stdpopsim can be used to systematically compare demographic inference methods, and we encourage a broader community of developers to contribute to this growing resource.

Antagonistic pleiotropy conceals molecular adaptations in changing environments.

The importance of positive selection in molecular evolution is debated. Evolution experiments under invariant laboratory conditions typically show a higher rate of nonsynonymous nucleotide change than that of synonymous change, demonstrating prevalent molecular adaptations. Natural evolution inferred from genomic comparisons, however, almost always exhibits the opposite pattern even among closely related conspecifics, indicative of a paucity of positive selection. We hypothesize that this apparent contradiction is at least in part attributable to ubiquitous and frequent environmental changes in nature, causing nonsynonymous mutations beneficial at one time deleterious soon after because of antagonistic pleiotropy and hindering their fixations relative to synonymous mutations despite continued population adaptations. To test this hypothesis, we performed yeast evolution experiments in changing and corresponding constant environments, followed by genome sequencing of the evolving populations. We observed a lower nonsynonymous to synonymous rate ratio in antagonistic changing environments than in the corresponding constant environments, and the population dynamics of mutations supports our hypothesis. These findings and the accompanying population genetic simulations suggest that molecular adaptation is consistently underestimated in nature due to the antagonistic fitness effects of mutations in changing environments.

Molecular assessment of natural disturbance regime in a temperate swamp forest

Advanced knowledge of forest disturbance regimes (the temporal and spatial patterns of periodic disturbances) has key relevance to implementation of effective managements for forest communities. For many forest ecosystems such knowledge is inadequate or is lacking. Spatially explicit simulations, testing effects of alternative disturbance regimes on fine-scale spatial genetic structure (FSGS) may provide information on key components of forest disturbance regimes. Here, we compare FSGSs assessed using 18 microsatellite markers with spatially explicit simulations, to characterise a natural disturbance regime in a European temperate swamp forest dominated by Alnus glutinosa. Individual-based population genetic simulations were carried out for three previously introduced models of disturbance regime, to determine which best fits the FSGS assessed by the kinship coefficient in a natural old-growth forest. We hypothesised a low proportion of trees removed (the gap dynamic model) accentuates high-relatedness among individuals over short distances (i.e. it accentuates FSGS), while a high proportion of trees eliminated by a disturbance (the cyclic succession model and the stand-replacing dynamic model) leads either to random FSGS (i.e. a point pattern of genotypes that cannot be differentiated from a complete spatial randomness) or a high spatial clustering of genotypes. Spatial genetic analyses found significant positive autocorrelations of genotypes at < 30 m in the natural old-growth forest, which indicates A. glutinosa individuals were more highly clustered than for a null model with complete spatial randomness. However, the strength of FSGS was low as indicated by a kinship coefficient < 0.02. Comparing the FSGSs of simulated disturbance regimes with the FSGS recorded in the old-growth forest, the gap dynamic model was the most probable scenario, with small canopy gaps occurring continually through time. The cyclic succession model and the stand-replacing dynamic model also led to non-random FSGSs but in both these models, the FSGS was stronger and more coarse-grained than that observed. The approach of comparing the FSGS recorded in a natural forest with simulations of disturbance regimes allows support to be found for (or against) previously hypothesised disturbance patterns. This approach should be considered alongside the well-established dendrochronological and palaeoecological methods.

Adaptive introgression during environmental change can weaken reproductive isolation

Anthropogenic climate change is an urgent threat to species diversity1,2. One aspect of this threat is the merging of species through increased hybridization3. The primary mechanism for this collapse is thought to be the weakening of ecologically mediated reproductive barriers, as demonstrated in cases of ‘reverse speciation’4,5. Here, we expand on this idea and show that adaptive introgression between species adapting to a shared, moving climatic optimum can readily weaken any reproductive barrier, including those that are completely independent of climate. Using genetically explicit forward-time simulations, we show that genetic linkage between alleles conferring adaptation to a changing climate and alleles conferring reproductive isolation (intrinsic and/or non-climatic extrinsic) can lead to adaptive introgression facilitating the homogenization of reproductive isolation alleles. This effect causes the decay of species boundaries across a broad and biologically realistic parameter space. We explore how the magnitude of this effect depends on the rate of climate change, the genetic architecture of adaptation, the initial degree of reproductive isolation, the degree to which reproductive isolation is intrinsic versus extrinsic and the mutation rate. These results highlight a previously unexplored effect of rapid climate change on species diversity. Increased hybridization has the potential to threaten species diversity. Here population genetic computer simulations show that climate-induced adaptive introgression could readily lead to hybridization even when reproductive isolation is independent from climate.

Limited fine-scale larval dispersal of the threatened brooding corals Heliopora spp. as evidenced by population genetics and numerical simulation

Coral reefs are declining worldwide due to regional and global anthropogenic stress. Sekisei Lagoon and areas in its vicinity are impacted by the same stresses. Sekisei Lagoon is the largest coral reef area in Japan covering an area of almost 15 km × 20 km. Being situated upstream of Kuroshio Current, coral reefs in this area could supply larvae to downstream such as the Ryukyu Archipelago, facilitating high biodiversity. To conserve corals in such an area, estimation of larval dispersal among populations is essential because connectivity via larval dispersal is associated with the time for replenishment after disturbance. Isolated coral populations have a higher risk of local extinction. The purpose of this study is to reveal genetically isolated populations of brooding corals with short larval durations in and around Sekisei Lagoon. For this purpose, we utilized both population genetic analysis and larval dispersal simulations. Based on genotyping of 929 colonies of blue coral, Heliopora spp., collected from 22 sites around Sekisei Lagoon, genetic analysis suggested self-recruitment with connectivity in the southeastern sites and northwestern sites. Particle-tracking larval dispersal simulations also revealed similar patterns, corroborating the genetic results. These results suggest that lagoonal scale connectivity among populations of Heliopora spp. in Sekisei Lagoon is limited, as expected based on its short larval duration. Our results highlight the fact that most of the Heliopora spp. populations are genetically isolated, even inside Sekisei Lagoon, implying that the likelihood of being replenished by larval supply from surrounding populations after a disturbance is low and have a higher risk of local extinction than other coral species with high larval dispersal potential.

Population genetic simulation study of power in association testing across genetic architectures and study designs.

While it is well established that genetics can be a major contributor to population variation of complex traits, the relative contributions of rare and common variants to phenotypic variation remains a matter of considerable debate. Here, we simulate genetic and phenotypic data across different case/control panel sampling strategies, sequencing methods, and genetic architecture models based on evolutionary forces to determine the statistical performance of rare variant association tests (RVATs) widely in use. We find that the highest statistical power of RVATs is achieved by sampling case/control individuals from the extremes of an underlying quantitative trait distribution. We also demonstrate that the use of genotyping arrays, in conjunction with imputation from a whole‐genome sequenced (WGS) reference panel, recovers the vast majority (90%) of the power that could be achieved by sequencing the case/control panel using current tools. Finally, we show that for dichotomous traits, the statistical performance of RVATs decreases as rare variants become more important in the trait architecture. Our results extend previous work to show that RVATs are insufficiently powered to make generalizable conclusions about the role of rare variants in dichotomous complex traits.

The global male-bias in sex ratio at birth is sustained by the sex ratio genotypes of replacement offspring

Throughout the world, there is a male-bias in the sex ratio at birth (SRB). It is not known whether this phenomenon has a genetic basis, though there is tentative evidence from genealogical and genomic studies that it may have. It has been proposed that the higher rate of male childhood mortality in humans is linked to the male-bias in SRB through parental investment, but this may only apply to facultative not genetic sex ratio adjustment. In previous population genetic models, elevated mortality in one sex (prior to breeding) has been shown not to affect the SRB, but these models did not consider the role of replacement births (i.e. births that only occur because a sibling died prematurely). In a set of population genetic modelling simulations, in which sex ratio is controlled by an autosomal gene expressed in the male line, this study shows that when there is replacement of dead offspring, this leads to a sustained bias in the SRB in the direction of the sex suffering the highest mortality. In the example of higher male mortality, this occurs, because replacement offspring are disproportionately drawn from fathers who were genetically predisposed to have initially had sons (because sons were more likely to die prematurely), and more likely to pass on male-biasing alleles to replacement offspring. To test the empirical basis for replacement births, an analysis of birth data from the Demographic and Health Survey program was conducted, which shows that parents do indeed tend to replace children who die prematurely.

LearnPopGen: An R package for population genetic simulation and numerical analysis.

Here, I briefly present a new R package called learnPopGen that has been designed primarily for the purposes of teaching evolutionary biology, population genetics, and evolutionary theory. Functions of the package can be used to conduct simulations and numerical analyses of a wide range of evolutionary phenomena that would typically be covered in advanced undergraduate through graduate‐level curricula in population genetics or evolution. For instance, learnPopGen functions can be used to visualize gene frequency changes through time under multiple deterministic and stochastic processes, to compute and animate the changes in phenotypic trait values or distributions under natural selection, to numerically analyze and graph the outcome of simple game theory models, and to plot coalescence within a population experiencing genetic drift, along with a number of other things. Functions have been designed to be maximally didactic and frequently employ compelling animated visualizations. Furthermore, it is straightforward to export plots and animations from R in the form of flat or animated graphics, or as videos. For maximum flexibility, students working with the package can run functions directly in R; however, instructors may choose to guide students less adept in the R environment to one of various web interfaces that I have built for a number of the functions of the package and that are already available online.