Activity 1b

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, the evolutionary history of a gene is not always identical to the history of the species themselves.

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. Again, this discordance is especially common when internal branches in the species tree are short, ancestral populations are large, or incomplete lineage sorting occurs.

2.1. Simulating Gene Trees with ms

We will use the coalescent simulator ms to generate gene trees evolving inside the species tree. The program ms, developed by Richard Hudson, is one of the most widely used coalescent simulators in population genetics and phylogenetics. 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 for modeling selection and recombination, scrm for scalable coalescent simulation, msprime for efficient modern coalescent simulation using tree sequences, SLiM for forward-time evolutionary simulation, and SimPhy 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
• nreps: number of replicate simulations
• [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. 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. Effects of Internal Branch Length on Gene Trees

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.

To explore this effect, we will use the following species tree:

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

This species tree shows that 3 and 4 are sister taxa, 2 diverged earlier, and 1 diverged earliest from the other taxa. The branch lengths represent divergence times measured in coalescent units.

In coalescent theory, branch lengths are typically scaled by effective population size and correspond to the expected amount of time available for gene lineages to coalesce. Longer branches provide more time for gene lineages to merge into common ancestors, while shorter branches increase the probability of incomplete lineage sorting and gene tree-species tree discordance.

---

Baseline simulation

The following command simulates gene trees under the species tree above:

ms 4 1000 -T \
-I 4 1 1 1 1 \
-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 with one sampled gene copy each
• The -ej option in ms 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:

ms 4 1000 -T \
-I 4 1 1 1 1 \
-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.

Exercise

1. Short Internal Branch and Gene Tree Discordance

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

(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 \
-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 \
-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 internal branch becomes extremely short.
   • Why might many gene trees disagree with the species tree even if the species tree is correct?
   • What evolutionary process causes this pattern?

Show Answer

When the internal branch is extremely short, ancestral gene lineages may not have enough time to coalesce before the next speciation event occurs. As a result, different genes can trace back through different ancestral histories, producing different gene tree topologies even when the species tree itself is correct.

This pattern is caused by an evolutionary process called incomplete lineage sorting (ILS).

---

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. Effects of Population Size on Gene Trees

In this activity, we will keep the species tree and divergence times fixed, but change the population size. This allows us to examine how population size affects gene tree-species tree discordance under the coalescent process. In ms, population size can be changed using the -t option. The -t parameter controls the population-scaled mutation rate, but in this activity we use it as a simple way to explore the effect of changing the scale of coalescent time. Larger population sizes make coalescence slower, which increases the chance that gene lineages fail to coalesce within short internal branches. This can lead to more incomplete lineage sorting.

We will use the same 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.

---

Small population size

ms 4 20 -T -t 1 \
-I 4 1 1 1 1 \
-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 \
-ej 1 4 3 \
-ej 2 3 2 \
-ej 4 2 1

---

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.