Activity 1

Table of Contents

1. Species Tree Simulation

• 1.1. Birth-only model: Yule process
• 1.2. Birth-death model
• 1.3. Comparing species trees

2. Gene Tree Simulation

• 2.1. Simulating Gene Trees with ms
• 2.2. Generate gene trees under a species tree
• 2.3. Effects of Branch Lengths and Population Size on Gene Trees
  • 2.3.1. Effects of Internal Branch Length
  • 2.3.2. Effects of Population Size

---

Goal

In this activity, we will simulate and compare:

1. Simulate species trees under different diversification models
2. Simulate gene trees evolving within species trees under the coalescent process
3. Explore how incomplete lineage sorting generates gene tree-species tree discordance

---

Activity 1: Species Trees and Gene Trees

1. Species Tree Simulation

1.1. Birth-only model: Yule process

The Yule process is a simple stochastic model of species evolution in which lineages split through speciation events over time, but extinction is not allowed. In this model, each lineage independently gives rise to a new descendant lineage at a constant speciation rate. As time progresses, the number of lineages increases through repeated branching events.

Because only branching events occur and extinction is excluded, the Yule process is also called a pure-birth model. The Yule process is one of the most widely used models for simulating species trees in phylogenetics because it provides a simple mathematical description of diversification over evolutionary time. Although real evolutionary histories often include extinction and changing diversification rates in different parts of the tree, the Yule model serves as a useful starting point for understanding how phylogenetic trees are generated and how branching patterns emerge.

In the script below, we simulate a species tree under the Yule model using the R package phytools [1]. First, we load the required R packages for phylogenetic analysis and tree manipulation. We set a random seed because the tree simulation is stochastic. It ensures that the same simulated tree is generated each time the script is run with the same seed, making the analysis reproducible.

Next, we simulate a species tree with 10 species (n = 10) under a constant speciation rate (b = 1.0) and no extinction (d = 0.0). The simulated tree is stored as an R object. Finally, we visualize the tree, add a figure title, and save the tree to a Newick file (named species_tree_yule.tre and stored in the working directory) so it can be used in other phylogenetic software.

# Load packages for phylogenetic analysis and tree manipulation
library(phytools)

# Set random seed for reproducibility
set.seed(1)

# Simulate a species tree under the Yule (pure-birth) model
species_tree_yule <- pbtree(
  b = 1.0,   # speciation rate
  d = 0.0,   # extinction rate
  n = 10     # number of species (tips)
)

# Plot the simulated tree
plot(species_tree_yule)

# Add a title to the figure
title("Species tree under a Yule model")

# Write the tree to a Newick file
write.tree(
  species_tree_yule,
  file = "species_tree_yule.tre"
)

Example output

Important

Exercise

1. Fix b = 1.0 and change n = 10 to n = 30.
   • How does the simulated tree differ from the tree generated above?
   • Is the larger tree easier or harder to interpret visually?

Show Answer

• When n is increased from 10 to 30, the simulated tree contains more species (tips) and becomes larger and more complex.

• The larger tree has (if binary):

  1. more branching events
  2. more internal nodes
  3. a denser visual appearance

• As the number of species increases, the tree often becomes harder to interpret visually because more branches and labels must be displayed.

2. Changing the speciation rate b from 1.0 to 0.2 and then to 2.0.
   • How does changing the speciation rate affect the topology of the tree?
   • Which simulation produces more recent branching events?

Show Answer

• (Usually,) lower speciation rates (for example, b = 0.2) tend to produce trees with fewer recent branching events and longer terminal branches (i.e., branches that are linked to the tips).

• (Usually,) higher speciation rates (for example, b = 2.0) tend to produce trees with more frequent branching events and shorter intervals between speciation events.

• The simulation with the higher speciation rate (b = 2.0) usually produces more recent branching events.

1.2. Birth-death model

The birth-death model extends the Yule process by allowing both speciation and extinction events. In this model, lineages can split into new descendant lineages through speciation, but they can also disappear through extinction. Because both speciation and extinction occur, the birth-death model is often considered a more realistic model of macroevolution than the pure-birth Yule process as extinction can strongly influence the shape, balance, and branching patterns of phylogenetic trees.

In the script below, we simulate a species tree under a birth-death model. We again set a random seed to ensure reproducibility. Next, we simulate a species tree with 10 species (n = 10) under a constant speciation rate (b = 1.0) and a nonzero extinction rate (d = 0.2). The simulated tree is then visualized, labeled, and written to a Newick file for use in other phylogenetic software.

# Set random seed for reproducibility
set.seed(2)

# Simulate a species tree under a birth-death model
species_tree_birth_death <- pbtree(
  b = 1.0,   # speciation rate
  d = 0.2,   # extinction rate
  n = 10     # number of species (tips)
)

# Plot the simulated tree
plot(species_tree_birth_death)

# Add a title to the figure
title("Species tree under a birth-death model")

# Write the tree to a Newick file
write.tree(
  species_tree_birth_death,
  file = "species_tree_birth_death.tre"
)

Example output

In the above birth-death simulation, n = 10 means that the process stops when there are 10 extant taxa, or 10 lineages alive at the present time. However, because extinction is allowed, some lineages may have existed in the past and then gone extinct before the present. Therefore, the full simulated tree can be non-ultrametric and contain more than 10 tips if extinct lineages are also shown. To focus only on the species alive at the present, we can remove extinct lineages and plot only the extant taxa adding extant.only=TRUE option (as shown below) when running pbtree() from phytools.

# Simulate a species tree under a birth-death model and only retain the extant taxa
species_tree_birth_death <- pbtree(
  b = 1.0,   # speciation rate
  d = 0.2,   # extinction rate
  n = 10,     # number of species (tips)
  extant.only=TRUE
)

# Plot the simulated tree
plot(species_tree_birth_death)

# Add a title to the figure
title("Species tree under a birth-death model")

# Write the tree to a Newick file
write.tree(
  species_tree_birth_death,
  file = "species_tree_birth_death.tre"
)

Example output

Important

Exercise

1. Fix b = 1.0 and increase the extinction rate d from 0.2 to 0.5.
   • Qualitatively, how does increasing extinction affect the shape of the tree?
   • Does the tree appear more balanced or imbalanced?

Note: In some simulations, no tree may be produced because all lineages go extinct before reaching the target number of extant taxa. In this case, simply re-run the simulation.

Show Answer

• Increasing the extinction rate generally produces trees with fewer surviving lineages and more uneven branching patterns.

• Higher extinction rates may create:

  1. longer isolated branches
  2. fewer recently diverged lineages
  3. a more irregular tree structure

• The tree often appears more imbalanced because some lineages survive while others disappear through extinction.

2. Fix d = 0.5 and increase the number of species from n = 10 to n = 30.
   • Qualitatively, how does the larger tree compare to the original tree?
   • Is the effect of extinction easier or harder to observe in larger trees?

Show Answer

• Increasing the number of species produces a larger and more complex tree with more branching events and internal nodes.

• In larger trees, the effects of extinction are often easier to observe because differences in branch lengths and tree balance become more visually apparent.

• However, larger trees can also become harder to interpret because of the increased number of branches and labels.

1.3. Comparing species trees

In this final species tree simulation, we compare several trees generated under different combinations of speciation rate b and extinction rate d. This allows us to see how changing diversification parameters can affect tree depth, branching patterns, and overall tree shape. Note because these simulations are stochastic, each tree is only one possible outcome under the chosen parameters. Even with the same values of b, d, and n, different runs can produce different tree shapes.

# Set random seed for reproducibility
set.seed(5)

# Simulate six species trees with different birth and death rates
tree1 <- pbtree(b = 0.5, d = 0.2, n = 30)
tree2 <- pbtree(b = 0.5, d = 0.5, n = 30)
tree3 <- pbtree(b = 1.0, d = 0.2, n = 30)
tree4 <- pbtree(b = 1.0, d = 0.5, n = 30)
tree5 <- pbtree(b = 2.0, d = 0.2, n = 30)
tree6 <- pbtree(b = 2.0, d = 0.5, n = 30)

# Plot the six trees together
par(mfrow = c(3, 2))
plot(tree1, main = "b = 0.5, d = 0.2")
plot(tree2, main = "b = 0.5, d = 0.5")
plot(tree3, main = "b = 1.0, d = 0.2")
plot(tree4, main = "b = 1.0, d = 0.5")
plot(tree5, main = "b = 2.0, d = 0.2")
plot(tree6, main = "b = 2.0, d = 0.5")

# Reset plotting layout
par(mfrow = c(1, 1))

# Compute tree depths
tree_depths <- c(
  tree1 = max(nodeHeights(tree1)),
  tree2 = max(nodeHeights(tree2)),
  tree3 = max(nodeHeights(tree3)),
  tree4 = max(nodeHeights(tree4)),
  tree5 = max(nodeHeights(tree5)),
  tree6 = max(nodeHeights(tree6))
)

# Plot tree depths
plot(
  1:6,
  tree_depths,
  xaxt = "n",
  xlab = "Simulation",
  ylab = "Tree depth",
  main = "Comparison of Tree Depths"
)

axis(1, at = 1:6, labels = 1:6)

Example output

Important

Exercise

1. Compare the six simulated trees.
   • Which tree is deepest?
   • How do different combinations of speciation rate b and extinction rate d affect tree shape?
   • Does using the same parameters always produce the same tree shape?

Show Answer

• The deepest tree is the one with the largest value in tree_depths. Likely to be one of the two cases where b = 0.5.

• Trees with higher speciation rates often have more frequent branching events because lineages split more rapidly.

• Trees with higher extinction rates may appear more irregular or imbalanced because some lineages disappear before the present.

• The overall tree shape depends on the combination of speciation rate b, extinction rate d, and stochastic variation in the simulation.

• Using the same parameters does not always produce the same tree shape unless the same random seed is used.

---

2. Gene Tree Simulation

Species trees describe the evolutionary history of populations or species. They represent the order and timing of speciation events through time. Gene trees, however, describe the ancestry of individual genes (or orthologous genes) sampled from those species. Because genes evolve within populations, and due to the stochastic nature of the coalescent process, different genes can produce different evolutionary trees. As a result, gene trees may disagree with the species tree, a phenomenon known as gene tree-species tree discordance. This discordance is especially common when internal branches in the species tree are short or ancestral populations are large, leading to so-called incomplete lineage sorting.

2.1. Simulating Gene Trees with ms

We will use the coalescent simulator ms to generate gene trees evolving inside the species tree. Note ms is a command-line software, so we will no longer use R for the rest of the tutorial but open your command line interface. The program ms, developed by Richard Hudson, is one of the most widely used coalescent simulators in population genetics and phylogenetics [2]. Although ms is a classic and highly influential simulator, many extensions and newer simulators have since been developed to model more complex evolutionary scenarios and improve computational efficiency. Examples include msms [3] for modeling selection and recombination, scrm [4] for scalable coalescent simulation, msprime [5] for efficient modern coalescent simulation using tree sequences, SLiM [6] for forward-time evolutionary simulation, and SimPhy [7] for species tree and gene tree simulation. In this activity, we use ms because it provides the foundational framework underlying many of these later coalescent simulation tools.

Unlike species tree simulations, which model speciation and extinction processes, ms models the ancestry of sampled gene copies backward in time under the coalescent process. Because coalescence is stochastic, different simulated genes may produce different gene trees even when evolving inside the same species tree. The basic structure of an ms command is:

ms nsam nreps [options]

This command generates random gene trees with n = nsam under the coalescent process where:

• nsam: number of sampled gene copies (i.e., number of tips)
• nreps: number of replicate simulations (i.e., number of gene trees)
• [options]: additional commands controlling output format, population structure, migration, divergence times, and demographic history

Simple example

ms 4 3 -T 

This command means:

• sample 4 gene copies (i.e., n = 4)
• simulate 3 independent gene trees
• print the simulated gene trees in Newick format using -T

Example output

skong@Sungsiks-iMac phylo-simulation % ./software/msdir/ms 4 3 -T
./software/msdir/ms 4 3 -T
58975 47308 65065

//
(4:0.383,(1:0.207,(2:0.115,3:0.115):0.092):0.175);

//
(1:0.700,(4:0.366,(2:0.232,3:0.232):0.134):0.334);

//
((2:0.024,4:0.024):0.209,(1:0.070,3:0.070):0.163);

The three integers on top represent the random number seeds used for the simulation. Because the coalescent process is random, different simulations can produce different tree topologies. As mentioned in the previous activity, these seeds make it possible to reproduce the exact same simulated gene trees later if needed. Then, each block beginning with // represents one independently simulated gene tree.

Alternatively, instead of printing the output directly to the terminal, the gene trees can be saved to the file for later analysis and visualization. This command redirects the simulated gene trees into a text file called gene_trees_basic.txt.:

ms 4 10 -T > gene_trees_basic.txt

2.2. Generate gene trees under a species tree

So far, we have not provided any information about the species tree when generating gene trees. In this subsection, we will provide species tree information to ms and simulate gene trees evolving within a given species tree. In particular, we will specify the divergence times and the taxa (or populations) involved in each divergence event. Consider the following species tree:

(1:4,(2:2,(3:1,4:1):1):2);

Because larger branch lengths represent longer evolutionary times, this species tree indicates that 3 and 4 are sister taxa, 2 diverged earlier from their common ancestor, and 1 diverged earliest from the other taxa. In this example, branch lengths are measured in coalescent units, a time scale commonly used in species tree models that incorporates both elapsed generations and effective population size. Consequently, these branch lengths do not correspond directly to calendar time. Nevertheless, longer branches provide more time for gene lineages to coalesce into common ancestors, whereas shorter branches increase the probability of incomplete lineage sorting (ILS) and gene tree-species tree discordance.

Baseline simulation

The following command simulates 1,000 independent gene trees under the species tree shown above, with one sampled lineage from each species:

ms 4 1000 -T \
-I 4 1 1 1 1 \
-n 1 2 -n 2 2 -n 3 2 -n 4 2 \
-ej 1 4 3 \
-ej 2 3 2 \
-ej 4 2 1

This command means:

• sample 4 gene copies (i.e., there are four taxa in the species tree)
• simulate 1000 independent gene trees
• -I: create 4 populations (taxa) with one sampled gene copy each
• -n sets the effective population size of each population. We double the population size here to match with the conventional multispecies coalescent model.
• The -ej option represents a population joining event backward in time under the coalescent process. Because coalescent simulations work backward from the present into the past, -ej merges one population into another ancestral population at a specified time. The general structure is:

-ej time source destination

which means at a given time in the past, merge the source population into the destination population.

For example:

-ej 1 4 3

means that at time 1, population 4 merges into population 3. This corresponds to taxa 3 and 4 sharing a common ancestor at that time.

Similarly:

-ej 2 3 2

means that at time 2, the ancestral lineage containing (3,4) merges into population 2. This corresponds to taxon 2 sharing a deeper common ancestor with (3,4).

Finally:

-ej 4 2 1

means that at time 4, the lineage containing (2,(3,4)) merges into population 1, corresponding to the root of the species tree. Together, these commands reconstruct the species tree backward in time through successive population mergers.

Visualizing Gene Trees Using DensiTree

The ms command prints gene trees in Newick format together with additional simulation information. To create a tree file readable by DensiTree, extract only the Newick trees. This can be done by adding | grep "^(" > filename in the command as shown in the example below:

ms 4 1000 -T \
-I 4 1 1 1 1 \
-n 1 2 -n 2 2 -n 3 2 -n 4 2 \
-ej 1 4 3 \
-ej 2 3 2 \
-ej 4 2 1 | grep "^(" > gene_trees.trees

Open gene_trees.trees in DensiTree.

DensiTree overlays many gene trees in a single figure. If most gene trees have the same topology, the plot will appear relatively clean and concentrated. If there is substantial gene tree discordance, the plot will appear more tangled because different loci support different evolutionary relationships. This allows us to visually compare how changing the internal branch length affects gene tree discordance. An example DensiTree output using gene_trees.trees should look something like:

Here, each line represents a gene tree, and darker regions indicate topologies that occur more frequently across simulations. The variation among trees reflects gene tree-species tree discordance caused by ILS.

Important

Exercise

Now try a different species tree topology:

((1:1,2:1):3,(3:2,4:2):2);

This species tree says that 1 and 2 are sister taxa, 3 and 4 are sister taxa, and the two pairs diverged from each other deeper in the past. Try to construct a command.

Example Command

ms 4 1000 -T -I 4 1 1 1 1 -n 1 2 -n 2 2 -n 3 2 -n 4 2 -ej 1 2 1 -ej 2 4 3 -ej 4 3 1 | grep "^(" > gene_trees_balanced.trees

Run the command and compare the simulated gene trees using DensiTree.

• Do all gene trees match the species tree topology?
• Which taxa are sisters in the species tree?
• Do 1 and 2 always appear as sister taxa in the simulated gene trees?
• Why might some gene trees differ from the species tree?

Show Answer

The species tree topology is ((1,2),(3,4)).

In the species tree, 1 and 2 are sister taxa, and 3 and 4 are sister taxa.

Simulated gene trees should generally match the species tree topology, but not all of them.

Even though 1 and 2 are sisters in the species tree, they may not always appear as sisters in the gene trees.

Gene trees can differ from the species tree because of incomplete lineage sorting. This happens when lineages fail to coalesce within the branches of the species tree before deeper speciation events.

Short internal branches in coalescent units usually increase the chance that gene trees differ from the species tree.

2.3. Effects of Branch Lengths and Population Size on Gene Trees

As you observed, gene trees do not always match the species tree from which they evolved. One of the primary causes of this discordance is ILS, which occurs when gene lineages fail to coalesce within the branches of the species tree before deeper speciation events. Two important factors that influence the amount of ILS are the length of internal branches in the species tree and the effective population size of ancestral populations. To explore these effects, we will first examine how changing internal branch lengths affects gene tree variation and then investigate how population size influences gene tree-species tree discordance.

2.3.1. Effects of Internal Branch Length

One of the most important factors affecting gene tree-species tree discordance is the length of internal branches in the species tree. Short internal branches provide less time for ancestral lineages to coalesce before the next speciation event, increasing the probability of incomplete lineage sorting. In contrast, long internal branches provide more time for coalescence, increasing agreement between gene trees and the species tree.

Consider the following species tree:

(1:4,(2:2,(3:1,4:1):1):2);

Gene trees generated under this species tree is saves as gene_trees.trees in the previous section. In DensiTree, the plot that visualizes these gene trees appears concentrated and well-defined. If substantial discordance exists, the plot appears more tangled because different loci support different evolutionary histories.

Important

Exercise

1. Short Internal Branch and Gene Tree Discordance

Generate a set of gene trees from the following species tree that contains a short internal branch:

(1:4,(2:1.2,(3:1,4:1):0.2):2.8);

• What happens to the length of the internal branch between (3,4) and 2?
• Do you expect more or less gene tree discordance?

Show Answer and Code

This can be done by modifying the baseline command.

ms 4 1000 -T \
-I 4 1 1 1 1 \
-n 1 2 -n 2 2 -n 3 2 -n 4 2 \
-ej 1 4 3 \
-ej 1.2 3 2 \
-ej 4 2 1

This shortens the internal branch separating (3,4) from 2. A shorter internal branch leaves less time for ancestral lineages to coalesce before the next speciation event. As a result, we expect more incomplete lineage sorting (ILS) and therefore greater gene tree discordance.

---

2. Long Internal Branch and Gene Tree Concordance

Create a long internal branch by generating a set of gene trees from the following species tree:

(1:4,(2:3,(3:1,4:1):2):1);

• What happens to the length of the internal branch between (3,4) and 2?
• Do you expect more or less agreement with the species tree?

Show Answer and Code

This can be done by modifying the baseline command.

ms 4 1000 -T \
-I 4 1 1 1 1 \
-n 1 2 -n 2 2 -n 3 2 -n 4 2 \
-ej 1 4 3 \
-ej 3 3 2 \
-ej 4 2 1

This lengthens the internal branch separating (3,4) from 2. A longer internal branch gives ancestral lineages more time to coalesce before the next speciation event. As a result, we expect less incomplete lineage sorting (ILS) and therefore greater agreement between gene trees and the species tree.

---

3. Suppose the species tree is a star tree (a completely unresolved tree with no internal branches).
   • What would you expect the relative frequencies of the possible gene tree topologies to be?

Show Answer

In a star tree, all lineages diverge simultaneously and there are no internal branches in which gene lineages can coalesce before deeper ancestral events. As a result, no particular gene tree topology is favored by the species tree.

Under the multispecies coalescent, all possible gene tree topologies are therefore expected to occur with equal probability. For four taxa, the three possible rooted gene tree topologies each occur approximately one-third of the time, giving an expected ratio of:

1 : 1 : 1

This illustrates that when the species tree contains no phylogenetic information, the gene trees do not preferentially support any particular topology.

---

4. Imagine analyzing thousands of genes from these species.
   • Would all genes necessarily support the same phylogenetic tree?
   • Why is it important to analyze many loci instead of only one gene?

Show Answer

Different genes may support different phylogenetic trees because gene trees evolve stochastically under the coalescent process. As a result, some genes may match the species tree, while others may show discordant relationships.

By analyzing many loci, we can estimate the underlying species history more accurately by averaging across gene tree variation.

2.3.2. Effects of Population Size

In the previous activity, we changed the length of internal branches while keeping population size fixed. We now keep the species tree fixed and instead vary population size. Larger populations tend to coalesce more slowly, increasing the chance that gene lineages fail to coalesce within species-tree branches before deeper speciation events. Consequently, larger populations often produce more incomplete lineage sorting and greater gene tree-species tree discordance.

We will use the following species tree:

(1:4,(2:2,(3:1,4:1):1):2);

The species tree and branch lengths remain fixed. Only the population-size-related parameter is changed.

In ms, population size can be modified using the -t parameter. Although -t technically controls the population-scaled mutation rate, we use it here as a simple way to explore how changes in the scale of coalescent time affect gene tree variation.

Small population size

ms 4 20 -T -t 1 \
-I 4 1 1 1 1 \
-n 1 2 -n 2 2 -n 3 2 -n 4 2 \
-ej 1 4 3 \
-ej 2 3 2 \
-ej 4 2 1

Large population size

ms 4 20 -T -t 10 \
-I 4 1 1 1 1 \
-n 1 2 -n 2 2 -n 3 2 -n 4 2 \
-ej 1 4 3 \
-ej 2 3 2 \
-ej 4 2 1

Important

Exercise

1. Compare the gene trees produced using -t 1 and -t 10.
   • Which setting produces more gene tree variation?
   • Which setting produces more agreement with the species tree?

Show Answer

• The larger value, -t 10, usually produces more gene tree variation.

• Larger population sizes mean that gene lineages tend to coalesce more slowly.

• Because lineages remain uncoalesced for longer, they may pass through multiple ancestral populations before coalescing.

• This increases incomplete lineage sorting and gene tree-species tree discordance.

• The smaller value, -t 1, usually produces more agreement with the species tree because lineages coalesce more quickly.

---

2. Why does increasing population size increase incomplete lineage sorting?

Show Answer

• In larger populations, two gene lineages are less likely to share a recent common ancestor quickly.

• This means lineages need more time to coalesce.

• If speciation events occur before gene lineages coalesce, the gene tree may differ from the species tree.

• Therefore, larger population sizes increase the probability of incomplete lineage sorting.

---

3. Suppose the population size becomes extremely large.
   • Would you expect gene trees to become more similar or more different from one another?
   • Why?

Show Answer

• Gene trees would generally become more different from one another.

• In very large populations, lineages coalesce slowly and may persist through several speciation events before merging.

• This increases gene tree variation and gene tree-species tree discordance.

---

4. Why is it important to analyze many loci when population sizes are large?

Show Answer

• In large populations, different genes may have very different evolutionary histories because of incomplete lineage sorting.

• A single gene tree may therefore provide a misleading estimate of the species history.

• Analyzing many loci helps estimate the underlying species tree more accurately by averaging across gene tree variation.

References

[1] Revell LJ. 2024. phytools 2.0: an updated R ecosystem for phylogenetic comparative methods (and other things). PeerJ 12:e16505. DOI: https://doi.org/10.7717/peerj.16505

[2] Hudson RR. 2002. Generating samples under a Wright–Fisher neutral model of genetic variation. Bioinformatics 18(2):337–338. DOI: https://doi.org/10.1093/bioinformatics/18.2.337

[3] Ewing G, Hermisson J. 2010. MSMS: a coalescent simulation program including recombination, demographic structure and selection at a single locus. Bioinformatics 26(16):2064–2065. DOI: https://doi.org/10.1093/bioinformatics/btq322

[4] Staab PR, Zhu S, Metzler D, Lunter G. 2015. scrm: efficiently simulating long sequences using the approximated coalescent with recombination. Bioinformatics 31(10):1680–1682. DOI: https://doi.org/10.1093/bioinformatics/btu861

[5] Baumdicker F, Bisschop G, Goldstein D, Gower G, Ragsdale AP, Tsambos G, Zhu S, Eldon B, Ellerman EC, Galloway JG, et al. 2022. Efficient ancestry and mutation simulation with msprime 1.0. Genetics 220(3):iyab229. DOI: https://doi.org/10.1093/genetics/iyab229

[6] Haller BC, Messer PW. 2023. SLiM 4: multispecies eco-evolutionary modeling. The American Naturalist 201(5):E127–E139. DOI: https://doi.org/10.1086/723601

[7] Mallo D, De Oliveira Martins L, Posada D. 2016. SimPhy: Phylogenomic simulation of gene, locus, and species trees. Systematic Biology 65(2):334–344. DOI: https://doi.org/10.1093/sysbio/syv082