Species delimitation and introgression in Penstemon

Speciation Genomics · PB 495/595 · Dasanthera exercise

0. Technical details

This notebook uses data from whole genomes of 18 plants across seven species of Penstemon from Stone & Wessinger (2024). The system involves a recent radiation of bee-pollinated and hummingbird-pollinated species.

Create a new directory for your work and navigate to it in R with setwd("your_directory")
All of the inputs are in data/, which you can download from drive and place in your working directory.
A few R packages will be used (some new some old): ape, ggplot2, dplyr, tidyr, patchwork, ggrepel, and maps.
The code is a bit dense in places. Go for it yourself if you want to understand the details, but I will lead the discussion today.
Because some of the analyses are dependant on software outside of R that may not easily work on your computers, those outputs have been provided for us to assess. I will run IQ-TREE and ASTRAL to show interactively to provide some discussion.
Some analyses are done with separate software, but we recreate the D-stats ourselves in R.
There are more questions than usual, but we will just attempt what we can in the time.

A little bit of setup to help us along the way. Ensure your packages are installed and loaded.

library(ape)         # trees and DNA sequences
library(ggplot2)     # plots
library(dplyr)       # filtering and summarising tables
library(tidyr)       # reshaping tables
library(patchwork)   # putting several plots side by side
library(ggrepel)     # labels that don't overlap
theme_set(theme_bw(base_size = 11))

# one colour per species, used throughout
species_cols <- c(`P. cardwellii` = "#1b9e77", `P. davidsonii` = "#d95f02",
                  `P. fruticosus` = "#7570b3", `P. newberryi` = "#e7298a",
                  `P. rupicola` = "#66a61e", `P. lyallii` = "grey55", `P. montanus` = "grey20")

# a colour legend for the tree plots (which are drawn with base R)
species_legend <- function(where = "bottomleft") {
  legend(where, legend = names(species_cols), text.col = species_cols,
         text.font = 3, bty = "n", cex = 0.8, inset = c(0, 0.08))
}

1. The plants and the data

Dasanthera penstemons are woody shrubs of rocky slopes from the Sierra Nevada to the northern Rockies. Most are bee-pollinated with blue-violet flowers. Two, P. newberryi and P. rupicola, have switched to hummingbird pollination and magenta flowers, and each forms a long-standing hybrid zone with the bee-pollinated P. davidsonii.

Stone & Wessinger sequenced 18 plants from different localities: several each of P. cardwellii, P. davidsonii, P. fruticosus, P. newberryi and P. rupicola, plus one each of P. lyallii and P. montanus as outgroups.

# read the table of plants: one row per plant
samples <- read.csv("data/samples.csv")

# a lookup from sample ID to species, e.g. sp["car_28"] gives "P. cardwellii"
sp <- setNames(samples$species, samples$sample_id)

# show the columns we need
samples[, c("sample_id", "species", "locality", "hybrid_zone", "role")]

   sample_id       species       locality hybrid_zone     role
1     car_28 P. cardwellii Mt. St. Helens          no  ingroup
2     car_91 P. cardwellii        Siuslaw          no  ingroup
3     car_92 P. cardwellii      Tillamook          no  ingroup
4     dav_87 P. davidsonii    Mt. Ashland          no  ingroup
5    dav_116 P. davidsonii    Crater Lake          no  ingroup
6    dav_117 P. davidsonii   Paulina Peak          no  ingroup
7    dav_118 P. fruticosus                         no  ingroup
8    fru_106 P. fruticosus    Mt. Rainier          no  ingroup
9     new_75  P. newberryi         Dorado          no  ingroup
10    new_80  P. newberryi         Plumas          no  ingroup
11    new_84  P. newberryi                         no  ingroup
12    new_85  P. newberryi    Castle Lake         yes  ingroup
13    rup_86   P. rupicola    Castle Lake         yes  ingroup
14    rup_98   P. rupicola    Rattlesnake          no  ingroup
15   rup_101   P. rupicola       Cle Elum          no  ingroup
16   rup_105   P. rupicola    Mt. Rainier          no  ingroup
17    lya_44    P. lyallii        Kaniksu          no outgroup
18    mon_61   P. montanus        Challis          no outgroup

The species column is the taxonomic assignment from the paper (its Table S1), and it is what we use to colour and group plants from here on. The three-letter prefix on each sample ID is just a collector’s code, and it doesn’t always match: dav_118 is assigned to P. fruticosus. Two plants, new_85 and rup_86, come from the same place, the Castle Lake hybrid zone. Keep an eye on them.

# keep only the 16 ingroup plants; the outgroups are left off the map
ingroup_samples <- filter(samples, role == "ingroup")

# outlines of the US states, used as a grey background (needs the maps package)
us_states <- map_data("state")

# nudge points by up to 0.25 degrees so plants from the same site don't hide each other
jit <- position_jitter(width = 0.25, height = 0.25, seed = 1)

# the map's limits: the range of the plants plus 1.5 degrees of margin
lon_limits <- range(ingroup_samples$longitude) + c(-1.5, 1.5)
lat_limits <- range(ingroup_samples$latitude) + c(-1.5, 1.5)

ggplot() +
  # layer 1: the state outlines
  geom_polygon(data = us_states, aes(long, lat, group = group),
               fill = "grey90", colour = "grey65", linewidth = 0.3) +
  # layer 2: one point per plant, coloured by species (triangles = Castle Lake)
  geom_point(data = ingroup_samples,
             aes(longitude, latitude, colour = species, shape = hybrid_zone),
             size = 3, position = jit) +
  # layer 3: the sample ID next to each point
  geom_text_repel(data = ingroup_samples,
                  aes(longitude, latitude, label = sample_id),
                  size = 3, show.legend = FALSE, seed = 1, position = jit) +
  scale_colour_manual(values = species_cols) +
  scale_shape_manual(values = c(no = 16, yes = 17)) +
  coord_quickmap(xlim = lon_limits, ylim = lat_limits) +
  labs(x = "longitude", y = "latitude", colour = "species",
       shape = "Castle Lake hybrid zone")

A map of the western United States in grey with state outlines, showing 16 collection sites from northern California to Washington and east to southeast Oregon, coloured by species. Two triangles, one magenta and one green, sit next to each other at Castle Lake in northern California. — Collection sites of the 16 ingroup plants on a grey map of the western US (outgroups not shown; points are jittered slightly so plants from the same site don’t hide each other). Colour is the species; triangles are the two plants from the Castle Lake hybrid zone.

data/genotypes.csv.gz has one row per plant and one column per SNP, thinned and pruned for linkage. Each value is the number of non-reference allele copies (0, 1 or 2), or NA if the genotype wasn’t called.

# read the genotype table: first column is the plant, the other columns are SNPs
geno <- read.csv("data/genotypes.csv.gz", check.names = FALSE)

# turn it into a matrix G with plants as rows and SNPs as columns
G <- as.matrix(geno[, -1])
rownames(G) <- geno$sample_id

dim(G)                      # number of plants and number of SNPs

[1]    18 50000

# Have a look at the first 5 snps from the 18 individuals
G[1:18,1:5]

        scaffold_1086:1805 scaffold_1086:2997 scaffold_1086:15143
car_28                   1                  2                   2
car_91                   1                  2                   2
car_92                   1                  2                   2
dav_116                  1                  2                  NA
dav_117                  1                  2                   2
dav_118                  1                  2                   2
dav_87                   1                  1                   2
fru_106                  1                  2                   2
lya_44                  NA                  1                   2
mon_61                  NA                  2                  NA
new_75                   2                  2                   2
new_80                   2                  2                   2
new_84                   1                  1                   2
new_85                   1                  2                   2
rup_101                  1                  2                   2
rup_105                  1                  2                   2
rup_86                   1                  2                   2
rup_98                   1                  2                   0
        scaffold_1086:15926 scaffold_1086:16005
car_28                    2                   1
car_91                    2                   2
car_92                    2                   2
dav_116                   1                   1
dav_117                   2                   2
dav_118                   2                   0
dav_87                    2                   0
fru_106                   2                   2
lya_44                    2                   0
mon_61                    2                  NA
new_75                    2                   2
new_80                    1                   1
new_84                    2                   0
new_85                    2                   2
rup_101                   2                   0
rup_105                   2                  NA
rup_86                    2                   0
rup_98                    2                   0

2. PCA

PCA again! It is an easy and first check to look at genetic similarities among closely related species. Since we are reading in the data as a table rather than from the VCF, we need to handle missing genotypes ourselves. We impute the missing genotypes by filling them with the SNP’s mean, since PCA can’t use NA.

# the ingroup plants only (the outgroups are removed)
ingroup <- samples$sample_id[samples$role == "ingroup"]
Gi <- G[ingroup, ]

# drop SNPs that don't vary among these plants, since they carry no information
snp_variance <- apply(Gi, 2, var, na.rm = TRUE)
Gi <- Gi[, snp_variance > 0]

# fill each missing genotype (NA) with the mean genotype of that SNP
for (j in 1:ncol(Gi)) {
  is_missing <- is.na(Gi[, j])
  Gi[is_missing, j] <- mean(Gi[, j], na.rm = TRUE)
}

# run the PCA: each individual gets a score on each principal component (PC)
pca <- prcomp(Gi)

# the percent of the variance each PC explains
pve <- round(100 * pca$sdev^2 / sum(pca$sdev^2), 1)

# put the scores for the first three PCs in a table, next to the sample information
pc_scores <- data.frame(sample_id = rownames(pca$x), pca$x[, 1:3])
pcs <- left_join(pc_scores, samples, by = "sample_id")

pve[1:3]

[1] 15.9 13.1 11.4

# a function that draws one PCA panel: PC "x" against PC "y"
pc_panel <- function(x, y) {
  # axis labels that include the percent variance, e.g. "PC1 (15.9%)"
  x_number <- as.integer(sub("PC", "", x))   # "PC1" -> 1
  y_number <- as.integer(sub("PC", "", y))
  x_label <- sprintf("%s (%.1f%%)", x, pve[x_number])
  y_label <- sprintf("%s (%.1f%%)", y, pve[y_number])

  ggplot(pcs, aes(.data[[x]], .data[[y]], colour = species)) +
    geom_point(aes(shape = hybrid_zone), size = 3) +
    geom_text_repel(aes(label = sample_id), size = 2.8, show.legend = FALSE, seed = 1) +
    scale_colour_manual(values = species_cols) +
    scale_shape_manual(values = c(no = 16, yes = 17)) +
    labs(x = x_label, y = y_label, colour = "species", shape = "Castle Lake")
}

# two panels side by side, sharing one legend
pc_panel("PC1", "PC2") + pc_panel("PC2", "PC3") + plot_layout(guides = "collect")

Two PCA scatter plots of 16 plants coloured by species. Left, PC1 against PC2: cardwellii and newberryi on the left, rupicola top right, davidsonii and fruticosus bottom right with dav_118 among the davidsonii plants. Right, PC2 against PC3: the two fruticosus plants, dav_118 and fru_106, are both low on PC3 while the three davidsonii plants are high. — PCA of the 16 ingroup plants. Colour is the species.

Q1: Look at both PCA panels. Are there a couple of species that like to stick close together?

3. Trees

A tree is consistent with how we think about evolution - individuals share common ancestors. A distance based tree does not utilize an explicit evolutionary model, but it can still get the tree right under some conditions. They are computationally easy and feasible for very large datasets.

Here, we calculate pairwise distances (p-distances) between all individuals. The distance between two individuals is the average difference over the SNPs both were called at (missing data does not count as a difference), divided by 2. The Neighbour-joining (NJ) algorithm can build a tree from these distances, and we root it on P. montanus.

# a function that computes the genotype distance between every pair of plants
gdist <- function(X) {
  n <- nrow(X)                                    # number of plants
  # start with a matrix of zeros: one row and one column per plant
  D <- matrix(0, n, n, dimnames = list(rownames(X), rownames(X)))

  # go through every pair of plants (i, j)
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      # difference in genotype (0, 1 or 2) at each SNP, ignoring SNPs where either plant is missing
      differences <- abs(X[i, ] - X[j, ])
      # average over SNPs, then divide by 2: roughly the proportion of allele copies that differ
      d <- mean(differences, na.rm = TRUE) / 2
      D[i, j] <- d
      D[j, i] <- d
    }
  }
  as.dist(D)
}

# build the neighbour-joining tree, then root it on the outgroup P. montanus
distance_matrix <- gdist(G)
nj_tree <- nj(distance_matrix)
nj_tree <- root(nj_tree, "mon_61", resolve.root = TRUE)

# ladderize() just tidies the order of the branches for drawing
nj_l <- ladderize(nj_tree)

# colour each tip label by the species of that plant
tip_colours <- species_cols[sp[nj_l$tip.label]]

plot(nj_l, tip.color = tip_colours, cex = 0.9, underscore = TRUE, no.margin = TRUE)
species_legend()
add.scale.bar()

A neighbour-joining tree of 18 plants with tips coloured by species. The tree has two main groups: davidsonii plus fruticosus, and rupicola plus cardwellii plus newberryi. Every species forms a clade, and the two fruticosus plants, dav_118 and fru_106, are sisters. — Neighbour-joining tree from genotype distances, rooted on *P. montanus*. Tip colour is the species.

Q2: Are there notable clades or relationships in the tree? Does each accepted species form its own clade?

A maximum-likelihood tree with IQ-TREE

A maximum-likelihood (ML) tree uses a model of sequence evolution. We could do this for NJ too, but it is common to assess uncertainty in an ML framework using bootstrap replicates too. data/snps.phy is an alignment of 14,488 variable SNPs, with heterozygous sites written as IUPAC ambiguity codes. Because it only has variable sites, the model needs an ascertainment-bias correction (+ASC), otherwise branch lengths would be overestimated. We will not worry about the ascertainment-bias correction, but we can estimate the ML tree with a model and look at some of the components we discussed in class. Below is the command I used:

iqtree3 -s data/snps.phy -st DNA -m MFP+ASC -B 1000 -T 4 --prefix penstemon_tree

This output already exists for you in the /data/trees directory.

# provide the path to the tree file
treefile <- "data/trees/iqtree_reference.treefile"

# read the ML tree and root it on P. montanus
# (edgelabel = TRUE keeps the bootstrap values on the right branches when re-rooting)
ml <- read.tree(treefile)
ml <- root(ml, "mon_61", resolve.root = TRUE, edgelabel = TRUE)

# We can check the newly estimated tree against the author's tree
authors <- read.tree("data/trees/concatML_authors.tre")
authors <- root(authors, "mon_61", resolve.root = TRUE)

# The Robinson-Foulds distance can be used to assess if trees are different. A 0 means the topology is the same.
dist.topo(unroot(ml), unroot(authors))

      tree1
tree2     0

dist.topo(unroot(ml), unroot(nj_tree))

      tree1
tree2     0

#Let's go ahead and plot the tree, and we can look at the branch length and node support
ml_l <- ladderize(ml)

# the node labels are bootstrap values stored as text; turn them into numbers
# (the root has no value, so it becomes NA and is left blank)
support <- suppressWarnings(as.numeric(ml_l$node.label))
support_labels <- ifelse(is.na(support), "", support)

tip_colours <- species_cols[sp[ml_l$tip.label]]
plot(ml_l, tip.color = tip_colours, cex = 0.9, underscore = TRUE, no.margin = TRUE)
nodelabels(support_labels, frame = "none", adj = c(1.2, -0.4), cex = 0.7)
species_legend()
add.scale.bar()

A maximum-likelihood tree of 18 plants with bootstrap values at the nodes, identical in topology to the neighbour-joining tree. The two fruticosus plants are sisters with 100 percent support; the node joining them to the davidsonii clade has 70 percent. — Maximum-likelihood tree from the SNP alignment, rooted on *P. montanus*. Tip colour is the species; node labels are ultrafast bootstrap support.

Every species clade has 100% bootstrap support. The weakest node is the clade of P. fruticosus and P. davidsonii. For these types of bootstrap values, we consider anything 95% and above strong support. Thus, there is some uncertainty in the data for this relationship, but the species hypotheses are well-supported.

Do we need a model? p-distances against ML distances

If sequences had diverged a lot more, some sites will have changed more than once and p-distances would underestimate the true amount of change. We call this scenario saturation. An ML tree’s branch lengths come from a model that corrects for those hidden changes, so the ML distance between two individuals should climb above the p-distance as divergence grows. You can plot the p-distance against the ML distance to assess if your data is saturated.

# read the alignment that IQ-TREE used
aln <- read.dna("data/snps.phy", format = "sequential")

# p-distance: the proportion of sites that differ between each pair of plants
# (pairwise.deletion = TRUE ignores sites with missing data in either plant)
p_dist <- dist.dna(aln, model = "raw", pairwise.deletion = TRUE)
p_dist <- as.matrix(p_dist)

# ML distance: add up the branch lengths along the path between two plants in the ML tree
ml_dist <- cophenetic(ml)
ml_dist <- ml_dist[rownames(p_dist), rownames(p_dist)]   # same plant order as p_dist

# list every pair of plants once, using the lower triangle of the matrix
pair <- which(lower.tri(p_dist), arr.ind = TRUE)         # row and column number of each pair
dd <- data.frame(a  = rownames(p_dist)[pair[, 1]],
                 b  = colnames(p_dist)[pair[, 2]],
                 p  = p_dist[pair],
                 ml = ml_dist[pair])

# look up the species of each plant in the pair
dd$species_a <- unname(sp[dd$a])
dd$species_b <- unname(sp[dd$b])

# label each pair: the later lines overwrite the earlier ones where they apply
outgroups <- c("P. lyallii", "P. montanus")
dd$type <- "different species"
dd$type[dd$species_a == dd$species_b] <- "same species"
dd$type[dd$species_a %in% outgroups | dd$species_b %in% outgroups] <- "involves an outgroup"

# summarise each type of pair; the ratio ml / p shows how much the model adds
dd |>
  group_by(type) |>
  summarise(pairs = n(),
            median_p = round(median(p), 3),
            median_ml = round(median(ml), 3),
            median_ratio = round(median(ml / p), 2))

# A tibble: 3 × 5
  type                 pairs median_p median_ml median_ratio
  <chr>                <int>    <dbl>     <dbl>        <dbl>
1 different species      101    0.263     0.308         1.14
2 involves an outgroup    33    0.347     0.356         1.08
3 same species            19    0.097     0.088         0.89

# correlation between the two distances
cor(dd$p, dd$ml)

[1] 0.9378177

# each point is one pair of individuals
# the one-to-one dashed line shows the expectation where distances are equal
ggplot(dd, aes(p, ml, colour = type)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_point(alpha = 0.7, size = 2) +
  coord_equal(xlim = c(0, 0.45), ylim = c(0, 0.45)) +
  labs(x = "p-distance (proportion of differing sites)",
       y = "ML distance (patristic, from the IQ-TREE tree)", colour = NULL)

A scatter plot of ML distance against p-distance for all pairs of plants, with a dashed 1:1 line. Points rise steadily with p-distance and mostly sit slightly above the line for pairs from different species and below it for pairs within a species, with no sign of the ML distance shooting upwards at large p-distances. — Pairwise p-distance against ML (patristic) distance for all 153 pairs of plants. The dashed line is the one-to-one line.

The relationship is close to linear (r = 0.94). Thus, no saturation. The NJ tree is often fine for a few closely-related species, but estimating trees from multiple approaches is usually necessary for publication.

Q3: When would you expect a model-based distance to matter much more than it does here?

The coalescent (ASTRAL) tree

The ML tree above is a concatenation tree. It treats all 14,488 SNPs as if they shared one history. But, different parts of the genome can have different histories, so the authors also built a coalescent species tree. They inferred a tree for each of 28,255 10 kb windows and used ASTRAL to find the species tree that agrees with the most gene-tree quartets. This builds off of our discussion about variation in relationships for the three-taxon trees.

The author’s provided the ASTRAL result as a NEXUS file that has some additional annotations. This includes the propotion of gene trees that recover a given node in the species tree. If the genome carried one history, these values (called q1 here) would be near 1. If gene trees carried no signal about the branch, q1 would be 1/3, as the three possible arrangements would be equally common.

# read the NEXUS file as lines of text and pick out the line that holds the tree
nex <- readLines("data/trees/astral_10kb_authors.tre")
tree_line <- grep("TREE 1 =", nex, value = TRUE)
s <- sub("^\\s*TREE 1 = ", "", tree_line)      # keep only the tree itself

# step 1: rename the plants, e.g. 'car 28-5 S435' becomes car_28
s <- gsub("'([a-z]{3}) ([0-9]+)-[0-9]+ S[0-9]+'", "\\1_\\2", s)

# step 2: make q1 the node label, e.g. "):0.5[&q1=0.6,...]" becomes ")0.6:0.5"
s <- gsub("\\)(:[0-9.eE-]+)\\[&q1=([0-9.eE-]+),[^]]*\\]", ")\\2\\1", s)

# step 3: delete the other annotations, which are the text inside square brackets
s <- gsub("\\[[^]]*\\]", "", s)

# read the cleaned text as a tree, and root it on P. montanus
astral <- read.tree(text = s)
astral <- root(astral, "mon_61", resolve.root = TRUE, edgelabel = TRUE)

# the node labels are q1 values stored as text; turn them into numbers
q1 <- suppressWarnings(as.numeric(astral$node.label))

# do the trees have the same topology? 0 means yes
dist.topo(unroot(astral), unroot(ml))

      tree1
tree2     0

dist.topo(unroot(astral), unroot(nj_tree))

      tree1
tree2     0

summary(q1)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.3823  0.4075  0.4997  0.5162  0.5774  0.8365       2

astral_l <- ladderize(astral)

# node labels: q1 to two decimals (the root has no value, so it is left blank)
q1_l <- suppressWarnings(as.numeric(astral_l$node.label))
q1_labels <- ifelse(is.na(q1_l), "", sprintf("%.2f", q1_l))

tip_colours <- species_cols[sp[astral_l$tip.label]]

# use.edge.length = FALSE draws the tree without branch lengths (a cladogram)
plot(astral_l, tip.color = tip_colours, cex = 0.9, use.edge.length = TRUE,
     underscore = TRUE, no.margin = TRUE)
nodelabels(q1_labels, frame = "none", adj = c(1.2, -0.4), cex = 0.7)
species_legend()

A species tree of 18 plants drawn as a cladogram, with q1 quartet support values at the nodes. The topology is the same as the ML tree. Support for most nodes is between about 0.4 and 0.65, with a higher value of 0.84 for the node uniting everything except the outgroups. — ASTRAL species tree from 10 kb window trees (authors’ result), drawn without branch lengths. Tip colour is the species; node labels are q1, the fraction of gene-tree quartets supporting the branch.

Q5: The ML tree gave 100% bootstrap support to every species clade. All of the species tree estimation approaches match. The q1 values across nodes imply variation across the genome though. How can both be true?

A note on bootstrap support - it asks a very specific question: how much uncertainty is there due to sampling error? When alignments get long, like this one, it becomes less likely to get disagreements among bootstrap replicates. The coalescent tree is based on gene tree variation though. Most nodes have q1 between 0.4 and 0.65, against 1/3 for no signal. The fruticosus node (dav_118, fru_106) has q1 = 0.52, and the node joining fruticosus to davidsonii is among the weakest at 0.39.

4. Tests of gene flow:Genome-wide D across the species

The ABBA-BABA test asks whether two non-sister species share more derived alleles than incomplete lineage sorting alone would give. Take three species related as ((P1, P2), P3) and an outgroup O carrying the ancestral allele A. At an ABBA site, P2 and P3 share the derived allele B. At a BABA site, P1 and P3 share it. Incomplete lineage sorting produces both equally often, and gene flow between P2 and P3 produces extra ABBA. Significance comes from a block jackknife along the genome, reported as a Z-score, with |Z| > 3 the usual threshold.

An external program, Dsuite, ran this test for every trio of the five species, with P. montanus as outgroup. It orders each trio so that D is positive. Groups in the Dsuite files are labelled with short species codes: car = P. cardwellii, dav = P. davidsonii (three plants), fru = P. fruticosus (two plants, including dav_118), new = P. newberryi and rup = P. rupicola.

# Dsuite's results: one row per trio of species
dtab <- read.csv("data/dsuite_trios_species.csv", check.names = FALSE)

# keep the columns we need, round them, and sort by Z-score (largest first)
dtab_short <- data.frame(P1 = dtab$P1, P2 = dtab$P2, P3 = dtab$P3,
                         D = round(dtab$Dstatistic, 3),
                         Z = round(dtab$`Z-score`, 1),
                         f4_ratio = round(dtab$`f4-ratio`, 3))
dtab_short[order(-dtab_short$Z), ]

    P1  P2  P3     D    Z f4_ratio
2  car new dav 0.088 18.8    0.088
9  dav new rup 0.076 12.3    0.082
1  dav fru car 0.048  6.8    0.043
7  fru dav new 0.042  5.6    0.044
8  dav fru rup 0.040  5.5    0.042
5  rup car fru 0.032  4.9    0.031
3  rup car dav 0.026  4.3    0.028
10 rup new fru 0.023  3.9    0.022
4  new car fru 0.010  2.2    0.008
6  car new rup 0.003  0.9    0.003

Q6: How many of the ten trios show significant gene flow? Which pair of species stands out most? Does the result for P. newberryi and P. rupicola surprise you?

Eight of the ten trios have Z > 3. The strongest is newberryi sharing more derived alleles with davidsonii than its sister cardwellii does (D = 0.088, Z = 18.8), the davidsonii–newberryi hybrid zone. The last row, car–new–rup, shows nothing (D = 0.003, Z = 0.9), which we take apart next.

5. ABBA-BABA and fd by hand

The equations

With allele frequencies instead of single genomes, each site contributes a fractional ABBA and BABA. Let $p_1, p_2, p_3$ be the derived-allele frequencies in P1, P2 and P3. When the outgroup is ancestral ($p_O = 0$):

\[\text{ABBA} = (1 - p_1)\,p_2\,p_3 \qquad \text{BABA} = p_1\,(1 - p_2)\,p_3\]

\[D = \frac{\sum \text{ABBA} - \sum \text{BABA}}{\sum \text{ABBA} + \sum \text{BABA}}\]

D is approximately 0 without gene flow and positive with P2–P3 gene flow. But in small windows it behaves erratically due to sampling error. fd (Martin et al. 2015) keeps the same numerator and divides by the ABBA − BABA you would get if P2 and P3 were randomly mating mixed. At each site, the donor frequency $p_D$ stands in for both P2 and P3 and is whichever has the higher derived-allele frequency:

\[f_d = \frac{\sum \left[(1-p_1)\,p_2\,p_3 \;-\; p_1\,(1-p_2)\,p_3\right]} {\sum \left[(1-p_1)\,p_D\,p_D \;-\; p_1\,(1-p_D)\,p_D\right]}, \qquad p_D = \max(p_2,\,p_3)\]

fd only makes sense when D > 0. fdM (Malinsky et al. 2015) fixes that by using the same numerator but switching the denominator with the sign of D, so it is symmetric around zero: positive means P2–P3 gene flow, negative means P1–P3.

\[f_{dM} = \frac{\sum \left[(1-p_1)\,p_2\,p_3 \;-\; p_1\,(1-p_2)\,p_3\right]}{\text{denominator}}, \qquad \text{denominator} = \begin{cases} \sum \left[(1-p_1)\,p_D\,p_D - p_1\,(1-p_D)\,p_D\right], \; p_D=\max(p_2,p_3) & D \ge 0 \\[6pt] -\sum \left[(1-p_D)\,p_2\,p_D - p_D\,(1-p_2)\,p_D\right], \; p_D=\max(p_1,p_3) & D < 0 \end{cases}\]

Where the outgroup is polymorphic ($p_O > 0$), Dsuite weights each term by $(1-p_O)$ and adds the same site counted with the roles of A and B swapped, weighted by $p_O$. That’s what the code below does, and the extra terms vanish when $p_O = 0$.

D by hand

data/trio_af_scaffold_1086.csv.gz has one row per SNP on the largest scaffold (2.1 million). Each column is the derived-allele frequency in one group:

column	group
`p_car`	P. cardwellii (P1)
`p_new`	P. newberryi (P2)
`p_rup`	P. rupicola (P3)
`p_new_noCL`, `p_rup_noCL`	the same two species without the Castle Lake plants (`NA` where no other plant was called)
`p_out`	the outgroup, P. montanus

First, a worked example for a single SNP. Suppose P1 carries only the ancestral allele (frequency of the derived allele p1 = 0), and P2 and P3 both carry only the derived allele (p2 = 1, p3 = 1). The outgroup is ancestral. This site is a perfect ABBA site:

ex_p1 <- 0
ex_p2 <- 1
ex_p3 <- 1

(1 - ex_p1) * ex_p2 * ex_p3     # ABBA: P1 ancestral, P2 and P3 derived. 1

[1] 1

ex_p1 * (1 - ex_p2) * ex_p3     # BABA: P1 and P3 derived, P2 ancestral. 0

[1] 0

Real SNPs have in-between frequencies (for example p2 = 0.5, because one of two plants carries the derived allele), so each SNP contributes a fraction of an ABBA and a fraction of a BABA. We add these up over all SNPs on the scaffold. The function below does this for every SNP at once.

# read the table: one row per SNP, with the derived-allele frequency in each group
af <- read.csv("data/trio_af_scaffold_1086.csv.gz")
nrow(af)     # number of SNPs

[1] 2143343

head(af)     # the first few rows

   pos p_car p_new  p_rup p_new_noCL p_rup_noCL p_out
1 1841   0.5 0.250 0.5000     0.1667     0.5000   0.0
2 1844   0.0 0.000 0.1667     0.0000     0.0000   0.0
3 2772   0.0 0.000 0.1250     0.0000     0.1667   0.0
4 2783   1.0 1.000 0.8750     1.0000     0.8333   0.5
5 2789   0.0 0.250 0.0000     0.1667     0.0000   0.0
6 2792   0.0 0.375 0.0000     0.3333     0.0000   0.0

# a function that gives the ABBA and BABA contribution of every SNP.
# p1, p2, p3 and pO are columns of derived-allele frequencies in P1, P2, P3 and the outgroup.
abba_baba <- function(p1, p2, p3, pO) {
  # the usual case, outgroup ancestral: P2 and P3 share the derived allele (ABBA),
  # or P1 and P3 share it (BABA). (1 - pO) is the chance the outgroup is ancestral.
  abba_main <- (1 - p1) * p2 * p3 * (1 - pO)
  baba_main <- p1 * (1 - p2) * p3 * (1 - pO)

  # if the outgroup is itself polymorphic (pO > 0), Dsuite also counts the site with the
  # roles of the two alleles swapped, weighted by pO. These terms are 0 when pO = 0.
  abba_swapped <- p1 * (1 - p2) * (1 - p3) * pO
  baba_swapped <- (1 - p1) * p2 * (1 - p3) * pO

  # return a table with one row per SNP
  data.frame(abba = abba_main + abba_swapped,
             baba = baba_main + baba_swapped)
}

# a function that turns the per-SNP table into D: add up ABBA and BABA over all SNPs
# (na.rm = TRUE skips SNPs with missing values)
calc_D <- function(sites) {
  abba_total <- sum(sites$abba, na.rm = TRUE)
  baba_total <- sum(sites$baba, na.rm = TRUE)
  (abba_total - baba_total) / (abba_total + baba_total)
}

# trio ((cardwellii, newberryi), rupicola), with all plants...
all_samples <- abba_baba(af$p_car, af$p_new, af$p_rup, af$p_out)

# ...and again without the two Castle Lake plants
no_castle <- abba_baba(af$p_car, af$p_new_noCL, af$p_rup_noCL, af$p_out)

c(all_samples = calc_D(all_samples), without_Castle_Lake = calc_D(no_castle))

        all_samples without_Castle_Lake 
        -0.00664271         -0.03844378

We can check our by-hand numbers against Dsuite on the same scaffold. Dsuite lists the trio as new car rup, because it always puts first whichever sister makes D positive, so its D has the opposite sign to yours.

# Dsuite's own result for the same scaffold
dsuite_check <- read.csv("data/dsuite_car_new_rup_scaffold_1086.csv", check.names = FALSE)
dsuite_check[, c("samples", "P1", "P2", "P3", "Dstatistic", "Z-score", "ABBA", "BABA")]

         samples  P1  P2  P3 Dstatistic  Z-score    ABBA    BABA
1            all new car rup 0.00664306 0.552154 51845.6 51161.3
2 no Castle Lake new car rup 0.03844450 3.279800 51992.1 48142.4

Q8: Removing the two hybrid zone plants flipped our D stat to be significant. Any thoughts on why including these could mask a signal of historical gene flow with this approach?

Some weaknesses about D in window-based tests

To find where introgression happened, you might compute D in windows. In a low-diversity window ABBA + BABA is tiny, so D can become very noisy. Here we compute D and fd in 100 kb windows, without the Castle Lake plants. Since fd only works when D > 0, we swap P1 and P2 so the trio becomes ((newberryi, cardwellii), rupicola). We also keep only sites where the outgroup is homozygous.

# step 1: keep the SNPs where the outgroup is ancestral (p_out is 0), and where both
# newberryi and rupicola still have a called plant once Castle Lake is removed
sites <- filter(af, p_out == 0, !is.na(p_new_noCL), !is.na(p_rup_noCL))

# step 2: name the three groups. P1 and P2 are swapped, so P1 = newberryi, P2 = cardwellii
p1 <- sites$p_new_noCL
p2 <- sites$p_car
p3 <- sites$p_rup_noCL

# step 3: the donor frequency pD. At each SNP it is the higher of p2 and p3.
# pmax() takes the larger of two values at each position
pD <- pmax(p2, p3)

# step 4: three numbers for every SNP
abba   <- (1 - p1) * p2 * p3                   # fraction of an ABBA site
baba   <- p1 * (1 - p2) * p3                   # fraction of a BABA site
fd_max <- (1 - p1) * pD * pD - p1 * (1 - pD) * pD   # ABBA - BABA if P2 and P3 were completely mixed

# step 5: a table of these, with the 100 kb window each SNP falls in
# (a SNP at position 250,000 is in window 2, and so on)
site_table <- data.frame(window = floor(sites$pos / 1e5),
                         abba = abba, baba = baba, fd_max = fd_max)

# step 6: add up the SNPs within each window
win <- site_table |>
  group_by(window) |>
  summarise(n_sites = n(),
            abba_total = sum(abba),
            baba_total = sum(baba),
            fd_max_total = sum(fd_max))

# step 7: D and fd for each window. They share the same numerator (ABBA - BABA);
# D divides by ABBA + BABA, fd divides by the maximum possible ABBA - BABA
win$numerator <- win$abba_total - win$baba_total
win$denom <- win$abba_total + win$baba_total
win$D  <- win$numerator / win$denom
win$fd <- win$numerator / win$fd_max_total

summary(win[, c("n_sites", "D", "fd")])

    n_sites           D                  fd          
 Min.   : 613   Min.   :-0.63984   Min.   :-0.41916  
 1st Qu.:2790   1st Qu.:-0.07722   1st Qu.:-0.01954  
 Median :3566   Median : 0.04465   Median : 0.01279  
 Mean   :3502   Mean   : 0.03852   Mean   : 0.01822  
 3rd Qu.:4220   3rd Qu.: 0.14780   3rd Qu.: 0.05205  
 Max.   :6450   Max.   : 0.78542   Max.   : 0.60411

# put D and fd in a single column called "value", with a column "statistic" saying which
# is which, so ggplot can draw one panel for each
win_long <- pivot_longer(win, cols = c(D, fd), names_to = "statistic", values_to = "value")

# each dot is one 100 kb window; the x axis is the amount of ABBA + BABA in the window
p <- ggplot(win_long, aes(x=statistic, y=value, fill=statistic)) + 
  geom_violin() + 
  geom_jitter(height = 0, width = 0.2, aes(color = n_sites)) + 
  scale_color_viridis_c(option = "plasma")
  print(p)

Violin plots showing the distribution of D and fd in 100 kb windows, with points colored by the number of sites. — Distribution of D and fd in 100 kb windows on scaffold 1086 for ((*newberryi*, *cardwellii*), *rupicola*), without the Castle Lake plants.

The main takeaway is that people are always playing with this type of statistic to make it more interpretable and robust. In this case, fd is about small genomic windows. But if you grasp the concept of D, it should translate to other scenarios.

8. Flower colour: one origin or two?

Both hummingbird-pollinated species, P. newberryi and P. rupicola, have magenta flowers, while their bee-pollinated relatives are blue-violet. The blue colour comes from delphinidin-type anthocyanin pigments, which need the enzyme F3′5′H. Magenta flowers contain only the simpler pink cyanidin pigments, so they are consistent with a loss of F3′5′H function. The paper will show some excepetions. In our trees P. newberryi is sister to P. cardwellii. P. rupicola sits outside them, so the two magenta species are not sisters. The paper considered three hhypotheses for the loss of F3′5′H:

Two independent losses, one in each species.
One loss in the common ancestor, followed by reversals in the species that make blue pigments again. A broken gene is hard to repair, so this seems unlikely.
One loss shared between newberryi and rupicola, by incomplete lineage sorting or by introgression (adaptive introgression of the magenta allele).

Results rejecting hypothesis three

Could magenta flowers have a single origin shared between P. newberryi and P. rupicola due to gene flow? We can test this by looking at fdM around the F3′5′H gene. If the magenta allele was shared by introgression, we would expect a peak at the F3′5′H gene.

The values are Dsuite fdM for (cardwellii, newberryi, rupicola) in 500-SNP windows sliding by 250 SNPs, for all samples and without the two Castle Lake plants. The dashed line is the 95th percentile of windows on the scaffold, which we call “unusually high”. Following the paper, we look at the 400 kb around the gene, and also rank that 400 kb against every other 400 kb block on the scaffold.

# where the gene is: its scaffold, start and end, and its midpoint
locus <- read.csv("data/f3p5h_locus.csv")
gene_mid <- (locus$start + locus$end) / 2

# fdM in windows for the trio (car, new, rup), for both sets of samples
nr <- read.csv("data/fdm_car_new_rup_500SNP.csv.gz")

# keep the windows on the scaffold that holds the gene, and add each window's midpoint
scaf <- filter(nr, scaffold == locus$scaffold)
scaf$mid <- (scaf$start + scaf$end) / 2

# the paper looked at 400 kb around the gene
block <- 4e5

# the start of every 400 kb block along the scaffold, sliding along 50 kb at a time
block_starts <- seq(min(scaf$mid), max(scaf$mid) - block, by = 5e4)

# the two sets of samples: "all" and "no Castle Lake"
sample_sets <- unique(scaf$samples)

# an empty table that we fill in, one row per sample set
results <- data.frame()

for (set in sample_sets) {
  # the windows for this sample set
  d <- scaf[scaf$samples == set, ]

  # fdM in the windows within 200 kb either side of the gene
  near_gene <- d$fdM[abs(d$mid - gene_mid) <= block / 2]

  # the mean fdM of every 400 kb block on the scaffold
  block_means <- numeric(length(block_starts))
  for (b in seq_along(block_starts)) {
    in_block <- d$mid >= block_starts[b] & d$mid < block_starts[b] + block
    block_means[b] <- mean(d$fdM[in_block])
  }

  # how does the gene's neighbourhood compare with the rest of the scaffold?
  one_row <- data.frame(
    samples = set,
    highest_window_near_gene = max(near_gene),
    scaffold_95th_percentile_window = unname(quantile(d$fdM, 0.95)),
    mean_near_gene = mean(near_gene),
    percentile_among_400kb_blocks = mean(block_means <= mean(near_gene)))
  results <- rbind(results, one_row)
}

# round the numbers (everything except the first column) for display
results[, -1] <- signif(results[, -1], 2)
results

         samples highest_window_near_gene scaffold_95th_percentile_window
1            all                    0.037                           0.120
2 no Castle Lake                    0.025                           0.089
  mean_near_gene percentile_among_400kb_blocks
1         0.0093                          0.64
2         0.0055                          0.74

# the windows within 3 Mb of the gene, which is the part we draw
near3 <- scaf[abs(scaf$mid - gene_mid) < 3e6, ]

# the 95th percentile of fdM on the whole scaffold, for each set of samples (the dashed line)
q95 <- aggregate(fdM ~ samples, data = scaf, FUN = quantile, probs = 0.95)

ggplot(near3, aes(mid / 1e6, fdM)) +
  # shaded band: the 400 kb around the gene
  annotate("rect", xmin = (gene_mid - block / 2) / 1e6, xmax = (gene_mid + block / 2) / 1e6,
           ymin = -Inf, ymax = Inf, fill = "#7570b3", alpha = 0.15) +
  geom_hline(yintercept = 0, colour = "grey70") +
  # dashed line: the 95th percentile of windows on the scaffold
  geom_hline(data = q95, aes(yintercept = fdM), linetype = "dashed", colour = "grey40") +
  # purple line: the gene
  geom_vline(xintercept = gene_mid / 1e6, colour = "#7570b3", linewidth = 0.8) +
  geom_line(colour = "grey30") +
  facet_wrap(~ samples) +
  labs(x = paste(locus$scaffold, "position (Mb)"), y = expression(f[dM]~"(car, new, rup)"))

Two line plots of fdM along six megabases of scaffold 1086 around the F3'5'H gene, for all samples and without the Castle Lake plants. Within the shaded 400 kb region around the gene fdM stays well below the dashed 95th-percentile line. A tall spike lies about 0.6 megabases to the right of the gene. — fdM for (*cardwellii*, *newberryi*, *rupicola*) in 500-SNP windows within 3 Mb of *F3′5′H* (purple line; shaded: the 400 kb around it). Dashed: the 95th percentile of windows on the scaffold.

None of the values really stand out and the 400kb around the gene are not differentiated from the background. Probably safe to conclude that there is no support for a shared allele via introgression.

Results rejecting hypothesis two

If the magenta allele moved between the two species, newberryi and rupicola should be close relatives in a tree of F3′5′H itself. We will just look at the gene tree for F3′5′H. Following the paper, we collapse branches with low support, since a single gene of about 1.5 kb has few informative sites.

# read the gene tree and root it on the two outgroups (which are each other's closest relatives)
gt <- read.tree("data/trees/f3p5h_genetree.treefile")
gt <- root(gt, c("lya_44", "mon_61"), resolve.root = TRUE, edgelabel = TRUE)

# the sample IDs of the newberryi plants and of the rupicola plants
new <- names(sp)[sp == "P. newberryi"]
rup <- names(sp)[sp == "P. rupicola"]

# is each species a clade (monophyletic) in this gene tree?
newberryi_monophyletic <- is.monophyletic(gt, new)
rupicola_monophyletic  <- is.monophyletic(gt, rup)

# the smallest clade that holds all eight newberryi and rupicola plants: how many plants is it?
mrca_node <- getMRCA(gt, c(new, rup))
smallest_clade <- extract.clade(gt, mrca_node)

c(newberryi_monophyletic = newberryi_monophyletic,
  rupicola_monophyletic  = rupicola_monophyletic,
  tips_in_smallest_clade_holding_all_new_and_rup = length(smallest_clade$tip.label))

                        newberryi_monophyletic 
                                             1 
                         rupicola_monophyletic 
                                             0 
tips_in_smallest_clade_holding_all_new_and_rup 
                                            16

# collapse branches with no SH-aLRT support, as in the paper.
# Node labels look like "85.2/85" (SH-aLRT/UFBoot): take the number before the "/"
shalrt <- suppressWarnings(as.numeric(sub("/.*", "", gt$node.label)))

# internal nodes are numbered after the tips, so add the number of tips
no_support_nodes <- which(shalrt == 0) + Ntip(gt)

# gt$edge lists the branches; its second column is the node at the end of each branch.
# Set the length of every branch leading to an unsupported node to 0...
gt_collapsed <- gt
is_no_support_branch <- gt$edge[, 2] %in% no_support_nodes
gt_collapsed$edge.length[is_no_support_branch] <- 0

# ...and merge the zero-length branches into multifurcations
gt_collapsed <- di2multi(gt_collapsed, tol = 1e-8)

gt_l <- ladderize(gt_collapsed)

# node labels, leaving the one called "Root" blank
node_labels <- gt_l$node.label
node_labels[node_labels == "Root"] <- ""

tip_colours <- species_cols[sp[gt_l$tip.label]]
plot(gt_l, tip.color = tip_colours, cex = 0.9, underscore = TRUE, no.margin = TRUE)
nodelabels(node_labels, frame = "none", adj = c(1.1, -0.4), cex = 0.7)
species_legend()
add.scale.bar()

A gene tree of 18 plants with most branches collapsed into multifurcations. The four newberryi plants form a group, but the rupicola plants are scattered: three sit with cardwellii and a fruticosus plant, and one sits near the newberryi group. — ML gene tree of *F3′5′H*, rooted on the outgroups, with low-support branches collapsed. Tip colour is the species; node labels are SH-aLRT/UFBoot.

Q9: Do the P. newberryi and P. rupicola plants group together in the gene tree? Why would hypothesis 2 predict this? Is something different about the Castle Lake plant rup_86?

Results supporting hypothesis 1

That will be covered during the discussion!

9. Takeaways

Named species are hypotheses the genomes can test. All five species form clades in the NJ, ML and ASTRAL trees, but the davidsonii–fruticosus relationship is the weakest (70% bootstrap, q1 = 0.39).
Species relationships can have high discordance. A high proportion of gene trees reflected relationships different from the species tree. This was likely due to historical gene flow between present-day species.
Hybrid-zone samples can affect tests of introgression. Two Castle Lake plants turned D = −0.037 (Z = 12.6) into D = 0.003.
Trees and introgression tests can be used to evaluate hypotheses of gene family evolution and speciation.

References

Stone BW, Wessinger CA (2024) Ecological diversification in an adaptive radiation of plants: the role of de novo mutation and introgression. Mol Biol Evol 41:msae007.
Malinsky M, Challis RJ, Tyers AM, et al. (2015) Genomic islands of speciation separate cichlid ecomorphs in an East African crater lake. Science 350:1493–1498.
Martin SH, Davey JW, Jiggins CD (2015) Evaluating the use of ABBA–BABA statistics to locate introgressed loci. Mol Biol Evol 32:244–257.
Malinsky M, Matschiner M, Svardal H (2021) Dsuite — fast D-statistics and related admixture evidence from VCF files. Mol Ecol Resour 21:584–595.
Cruickshank TE, Hahn MW (2014) Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow. Mol Ecol 23:3133–3157.

--- title: "Species delimitation and introgression in *Penstemon*" subtitle: "Speciation Genomics · PB 495/595 · Dasanthera exercise" format: html: toc: true toc-depth: 2 code-fold: false code-tools: true embed-resources: true fig-width: 7 fig-height: 5 --- # 0. Technical details This notebook uses data from whole genomes of 18 plants across seven species of *Penstemon* from [Stone & Wessinger (2024)](https://doi.org/10.1093/molbev/msae007). The system involves a recent radiation of bee-pollinated and hummingbird-pollinated species. - Create a new directory for your work and navigate to it in R with `setwd("your_directory")` - All of the inputs are in `data/`, which you can download from [drive](https://drive.google.com/drive/folders/13_gjS4YGlX1xHpnncTqHwb06wNHo3bqW?usp=sharing) and place in your working directory. - A few R packages will be used (some new some old): `ape`, `ggplot2`, `dplyr`, `tidyr`, `patchwork`, `ggrepel`, and `maps`. - The code is a bit dense in places. Go for it yourself if you want to understand the details, but I will lead the discussion today. - Because some of the analyses are dependant on software outside of R that may not easily work on your computers, those outputs have been provided for us to assess. I will run IQ-TREE and ASTRAL to show interactively to provide some discussion. - Some analyses are done with separate software, but we recreate the D-stats ourselves in R. - There are more questions than usual, but we will just attempt what we can in the time. A little bit of setup to help us along the way. Ensure your packages are installed and loaded. ```{r} #| label: setup #| message: false #| warning: false library(ape) # trees and DNA sequences library(ggplot2) # plots library(dplyr) # filtering and summarising tables library(tidyr) # reshaping tables library(patchwork) # putting several plots side by side library(ggrepel) # labels that don't overlap theme_set(theme_bw(base_size = 11)) # one colour per species, used throughout species_cols <- c(`P. cardwellii` = "#1b9e77", `P. davidsonii` = "#d95f02", `P. fruticosus` = "#7570b3", `P. newberryi` = "#e7298a", `P. rupicola` = "#66a61e", `P. lyallii` = "grey55", `P. montanus` = "grey20") # a colour legend for the tree plots (which are drawn with base R) species_legend <- function(where = "bottomleft") { legend(where, legend = names(species_cols), text.col = species_cols, text.font = 3, bty = "n", cex = 0.8, inset = c(0, 0.08)) } ``` # 1. The plants and the data *Dasanthera* penstemons are woody shrubs of rocky slopes from the Sierra Nevada to the northern Rockies. Most are bee-pollinated with blue-violet flowers. Two, *P. newberryi* and *P. rupicola*, have switched to hummingbird pollination and magenta flowers, and each forms a long-standing **hybrid zone** with the bee-pollinated *P. davidsonii*. Stone & Wessinger sequenced 18 plants from different localities: several each of *P. cardwellii*, *P. davidsonii*, *P. fruticosus*, *P. newberryi* and *P. rupicola*, plus one each of *P. lyallii* and *P. montanus* as outgroups. ```{r} #| label: samples # read the table of plants: one row per plant samples <- read.csv("data/samples.csv") # a lookup from sample ID to species, e.g. sp["car_28"] gives "P. cardwellii" sp <- setNames(samples$species, samples$sample_id) # show the columns we need samples[, c("sample_id", "species", "locality", "hybrid_zone", "role")] ``` The `species` column is the taxonomic assignment from the paper (its Table S1), and it is what we use to colour and group plants from here on. The three-letter prefix on each sample ID is just a collector's code, and it doesn't always match: `dav_118` is assigned to *P. fruticosus*. Two plants, `new_85` and `rup_86`, come from the same place, the **Castle Lake** hybrid zone. Keep an eye on them. ```{r} #| label: sample-map #| eval: false # keep only the 16 ingroup plants; the outgroups are left off the map ingroup_samples <- filter(samples, role == "ingroup") # outlines of the US states, used as a grey background (needs the maps package) us_states <- map_data("state") # nudge points by up to 0.25 degrees so plants from the same site don't hide each other jit <- position_jitter(width = 0.25, height = 0.25, seed = 1) # the map's limits: the range of the plants plus 1.5 degrees of margin lon_limits <- range(ingroup_samples$longitude) + c(-1.5, 1.5) lat_limits <- range(ingroup_samples$latitude) + c(-1.5, 1.5) ggplot() + # layer 1: the state outlines geom_polygon(data = us_states, aes(long, lat, group = group), fill = "grey90", colour = "grey65", linewidth = 0.3) + # layer 2: one point per plant, coloured by species (triangles = Castle Lake) geom_point(data = ingroup_samples, aes(longitude, latitude, colour = species, shape = hybrid_zone), size = 3, position = jit) + # layer 3: the sample ID next to each point geom_text_repel(data = ingroup_samples, aes(longitude, latitude, label = sample_id), size = 3, show.legend = FALSE, seed = 1, position = jit) + scale_colour_manual(values = species_cols) + scale_shape_manual(values = c(no = 16, yes = 17)) + coord_quickmap(xlim = lon_limits, ylim = lat_limits) + labs(x = "longitude", y = "latitude", colour = "species", shape = "Castle Lake hybrid zone") ``` ![Collection sites of the 16 ingroup plants on a grey map of the western US (outgroups not shown; points are jittered slightly so plants from the same site don't hide each other). Colour is the species; triangles are the two plants from the Castle Lake hybrid zone.](figures/sample-map.png){fig-alt="A map of the western United States in grey with state outlines, showing 16 collection sites from northern California to Washington and east to southeast Oregon, coloured by species. Two triangles, one magenta and one green, sit next to each other at Castle Lake in northern California."} `data/genotypes.csv.gz` has one row per plant and one column per SNP, thinned and pruned for linkage. Each value is the number of non-reference allele copies (0, 1 or 2), or `NA` if the genotype wasn't called. ```{r} #| label: genotypes # read the genotype table: first column is the plant, the other columns are SNPs geno <- read.csv("data/genotypes.csv.gz", check.names = FALSE) # turn it into a matrix G with plants as rows and SNPs as columns G <- as.matrix(geno[, -1]) rownames(G) <- geno$sample_id dim(G) # number of plants and number of SNPs # Have a look at the first 5 snps from the 18 individuals G[1:18,1:5] ``` # 2. PCA PCA again! It is an easy and first check to look at genetic similarities among closely related species. Since we are reading in the data as a table rather than from the VCF, we need to handle missing genotypes ourselves. We *impute* the missing genotypes by filling them with the SNP's mean, since PCA can't use `NA`. ```{r} #| label: pca # the ingroup plants only (the outgroups are removed) ingroup <- samples$sample_id[samples$role == "ingroup"] Gi <- G[ingroup, ] # drop SNPs that don't vary among these plants, since they carry no information snp_variance <- apply(Gi, 2, var, na.rm = TRUE) Gi <- Gi[, snp_variance > 0] # fill each missing genotype (NA) with the mean genotype of that SNP for (j in 1:ncol(Gi)) { is_missing <- is.na(Gi[, j]) Gi[is_missing, j] <- mean(Gi[, j], na.rm = TRUE) } # run the PCA: each individual gets a score on each principal component (PC) pca <- prcomp(Gi) # the percent of the variance each PC explains pve <- round(100 * pca$sdev^2 / sum(pca$sdev^2), 1) # put the scores for the first three PCs in a table, next to the sample information pc_scores <- data.frame(sample_id = rownames(pca$x), pca$x[, 1:3]) pcs <- left_join(pc_scores, samples, by = "sample_id") pve[1:3] ``` ```{r} #| label: pca-plot #| eval: false #| fig-width: 10 #| fig-height: 4.5 # a function that draws one PCA panel: PC "x" against PC "y" pc_panel <- function(x, y) { # axis labels that include the percent variance, e.g. "PC1 (15.9%)" x_number <- as.integer(sub("PC", "", x)) # "PC1" -> 1 y_number <- as.integer(sub("PC", "", y)) x_label <- sprintf("%s (%.1f%%)", x, pve[x_number]) y_label <- sprintf("%s (%.1f%%)", y, pve[y_number]) ggplot(pcs, aes(.data[[x]], .data[[y]], colour = species)) + geom_point(aes(shape = hybrid_zone), size = 3) + geom_text_repel(aes(label = sample_id), size = 2.8, show.legend = FALSE, seed = 1) + scale_colour_manual(values = species_cols) + scale_shape_manual(values = c(no = 16, yes = 17)) + labs(x = x_label, y = y_label, colour = "species", shape = "Castle Lake") } # two panels side by side, sharing one legend pc_panel("PC1", "PC2") + pc_panel("PC2", "PC3") + plot_layout(guides = "collect") ``` ![PCA of the 16 ingroup plants. Colour is the species.](figures/pca-plot.png){fig-alt="Two PCA scatter plots of 16 plants coloured by species. Left, PC1 against PC2: cardwellii and newberryi on the left, rupicola top right, davidsonii and fruticosus bottom right with dav_118 among the davidsonii plants. Right, PC2 against PC3: the two fruticosus plants, dav_118 and fru_106, are both low on PC3 while the three davidsonii plants are high."} **Q1: Look at both PCA panels. Are there a couple of species that like to stick close together?** # 3. Trees A tree is consistent with how we think about evolution - individuals share common ancestors. A distance based tree does not utilize an explicit evolutionary model, but it can still get the tree right under some conditions. They are computationally easy and feasible for very large datasets. Here, we calculate pairwise distances (**p-distances**) between all individuals. The distance between two individuals is the average difference over the SNPs both were called at (missing data does not count as a difference), divided by 2. The **Neighbour-joining (NJ)** algorithm can build a tree from these distances, and we root it on *P. montanus*. ```{r} #| label: nj # a function that computes the genotype distance between every pair of plants gdist <- function(X) { n <- nrow(X) # number of plants # start with a matrix of zeros: one row and one column per plant D <- matrix(0, n, n, dimnames = list(rownames(X), rownames(X))) # go through every pair of plants (i, j) for (i in 1:(n - 1)) { for (j in (i + 1):n) { # difference in genotype (0, 1 or 2) at each SNP, ignoring SNPs where either plant is missing differences <- abs(X[i, ] - X[j, ]) # average over SNPs, then divide by 2: roughly the proportion of allele copies that differ d <- mean(differences, na.rm = TRUE) / 2 D[i, j] <- d D[j, i] <- d } } as.dist(D) } # build the neighbour-joining tree, then root it on the outgroup P. montanus distance_matrix <- gdist(G) nj_tree <- nj(distance_matrix) nj_tree <- root(nj_tree, "mon_61", resolve.root = TRUE) ``` ```{r} #| label: nj-plot #| eval: false #| fig-height: 5.5 # ladderize() just tidies the order of the branches for drawing nj_l <- ladderize(nj_tree) # colour each tip label by the species of that plant tip_colours <- species_cols[sp[nj_l$tip.label]] plot(nj_l, tip.color = tip_colours, cex = 0.9, underscore = TRUE, no.margin = TRUE) species_legend() add.scale.bar() ``` ![Neighbour-joining tree from genotype distances, rooted on *P. montanus*. Tip colour is the species.](figures/nj-plot.png){fig-alt="A neighbour-joining tree of 18 plants with tips coloured by species. The tree has two main groups: davidsonii plus fruticosus, and rupicola plus cardwellii plus newberryi. Every species forms a clade, and the two fruticosus plants, dav_118 and fru_106, are sisters."} **Q2: Are there notable clades or relationships in the tree? Does each accepted species form its own clade?** ## A maximum-likelihood tree with IQ-TREE A maximum-likelihood (ML) tree uses a model of sequence evolution. We could do this for NJ too, but it is common to assess uncertainty in an ML framework using **bootstrap replicates** too. `data/snps.phy` is an alignment of 14,488 variable SNPs, with heterozygous sites written as IUPAC ambiguity codes. Because it only has ** variable sites**, the model needs an ascertainment-bias correction (`+ASC`), otherwise branch lengths would be overestimated. We will not worry about the ascertainment-bias correction, but we can estimate the ML tree with a model and look at some of the components we discussed in class. Below is the command I used: ```bash iqtree3 -s data/snps.phy -st DNA -m MFP+ASC -B 1000 -T 4 --prefix penstemon_tree ``` This output already exists for you in the `/data/trees` directory. ```{r} #| label: ml # provide the path to the tree file treefile <- "data/trees/iqtree_reference.treefile" # read the ML tree and root it on P. montanus # (edgelabel = TRUE keeps the bootstrap values on the right branches when re-rooting) ml <- read.tree(treefile) ml <- root(ml, "mon_61", resolve.root = TRUE, edgelabel = TRUE) # We can check the newly estimated tree against the author's tree authors <- read.tree("data/trees/concatML_authors.tre") authors <- root(authors, "mon_61", resolve.root = TRUE) # The Robinson-Foulds distance can be used to assess if trees are different. A 0 means the topology is the same. dist.topo(unroot(ml), unroot(authors)) dist.topo(unroot(ml), unroot(nj_tree)) ``` ```{r} #| label: ml-plot #| eval: false #| fig-height: 5.5 #Let's go ahead and plot the tree, and we can look at the branch length and node support ml_l <- ladderize(ml) # the node labels are bootstrap values stored as text; turn them into numbers # (the root has no value, so it becomes NA and is left blank) support <- suppressWarnings(as.numeric(ml_l$node.label)) support_labels <- ifelse(is.na(support), "", support) tip_colours <- species_cols[sp[ml_l$tip.label]] plot(ml_l, tip.color = tip_colours, cex = 0.9, underscore = TRUE, no.margin = TRUE) nodelabels(support_labels, frame = "none", adj = c(1.2, -0.4), cex = 0.7) species_legend() add.scale.bar() ``` ![Maximum-likelihood tree from the SNP alignment, rooted on *P. montanus*. Tip colour is the species; node labels are ultrafast bootstrap support.](figures/ml-plot.png){fig-alt="A maximum-likelihood tree of 18 plants with bootstrap values at the nodes, identical in topology to the neighbour-joining tree. The two fruticosus plants are sisters with 100 percent support; the node joining them to the davidsonii clade has 70 percent."} Every species clade has 100% bootstrap support. The weakest node is the clade of *P. fruticosus* and *P. davidsonii*. For these types of bootstrap values, we consider anything 95% and above strong support. Thus, there is some uncertainty in the data for this relationship, but the species hypotheses are well-supported. ## Do we need a model? p-distances against ML distances If sequences had diverged a lot more, some sites will have changed more than once and p-distances would underestimate the true amount of change. We call this scenario **saturation**. An ML tree's branch lengths come from a model that corrects for those hidden changes, so the **ML distance** between two individuals should climb above the p-distance as divergence grows. You can plot the p-distance against the ML distance to assess if your data is saturated. ```{r} #| label: ml-vs-p # read the alignment that IQ-TREE used aln <- read.dna("data/snps.phy", format = "sequential") # p-distance: the proportion of sites that differ between each pair of plants # (pairwise.deletion = TRUE ignores sites with missing data in either plant) p_dist <- dist.dna(aln, model = "raw", pairwise.deletion = TRUE) p_dist <- as.matrix(p_dist) # ML distance: add up the branch lengths along the path between two plants in the ML tree ml_dist <- cophenetic(ml) ml_dist <- ml_dist[rownames(p_dist), rownames(p_dist)] # same plant order as p_dist # list every pair of plants once, using the lower triangle of the matrix pair <- which(lower.tri(p_dist), arr.ind = TRUE) # row and column number of each pair dd <- data.frame(a = rownames(p_dist)[pair[, 1]], b = colnames(p_dist)[pair[, 2]], p = p_dist[pair], ml = ml_dist[pair]) # look up the species of each plant in the pair dd$species_a <- unname(sp[dd$a]) dd$species_b <- unname(sp[dd$b]) # label each pair: the later lines overwrite the earlier ones where they apply outgroups <- c("P. lyallii", "P. montanus") dd$type <- "different species" dd$type[dd$species_a == dd$species_b] <- "same species" dd$type[dd$species_a %in% outgroups | dd$species_b %in% outgroups] <- "involves an outgroup" # summarise each type of pair; the ratio ml / p shows how much the model adds dd |> group_by(type) |> summarise(pairs = n(), median_p = round(median(p), 3), median_ml = round(median(ml), 3), median_ratio = round(median(ml / p), 2)) # correlation between the two distances cor(dd$p, dd$ml) ``` ```{r} #| label: ml-vs-p-plot #| eval: false #| fig-width: 6 #| fig-height: 5.5 # each point is one pair of individuals # the one-to-one dashed line shows the expectation where distances are equal ggplot(dd, aes(p, ml, colour = type)) + geom_abline(slope = 1, intercept = 0, linetype = "dashed") + geom_point(alpha = 0.7, size = 2) + coord_equal(xlim = c(0, 0.45), ylim = c(0, 0.45)) + labs(x = "p-distance (proportion of differing sites)", y = "ML distance (patristic, from the IQ-TREE tree)", colour = NULL) ``` ![Pairwise p-distance against ML (patristic) distance for all 153 pairs of plants. The dashed line is the one-to-one line.](figures/ml-vs-p-plot.png){fig-alt="A scatter plot of ML distance against p-distance for all pairs of plants, with a dashed 1:1 line. Points rise steadily with p-distance and mostly sit slightly above the line for pairs from different species and below it for pairs within a species, with no sign of the ML distance shooting upwards at large p-distances."} The relationship is close to linear (r = 0.94). Thus, no saturation. The NJ tree is often fine for a few closely-related species, but estimating trees from multiple approaches is usually necessary for publication. **Q3: When would you expect a model-based distance to matter much more than it does here?** ## The coalescent (ASTRAL) tree The ML tree above is a **concatenation** tree. It treats all 14,488 SNPs as if they shared one history. But, different parts of the genome can have different histories, so the authors also built a **coalescent** species tree. They inferred a tree for each of 28,255 10 kb windows and used ASTRAL to find the species tree that agrees with the most gene-tree **quartets**. This builds off of our discussion about variation in relationships for the three-taxon trees. The author's provided the ASTRAL result as a *NEXUS* file that has some additional annotations. This includes the propotion of gene trees that recover a given node in the species tree. If the genome carried one history, these values (called q1 here) would be near 1. If gene trees carried no signal about the branch, q1 would be 1/3, as the three possible arrangements would be equally common. ```{r} #| label: astral # read the NEXUS file as lines of text and pick out the line that holds the tree nex <- readLines("data/trees/astral_10kb_authors.tre") tree_line <- grep("TREE 1 =", nex, value = TRUE) s <- sub("^\\s*TREE 1 = ", "", tree_line) # keep only the tree itself # step 1: rename the plants, e.g. 'car 28-5 S435' becomes car_28 s <- gsub("'([a-z]{3}) ([0-9]+)-[0-9]+ S[0-9]+'", "\\1_\\2", s) # step 2: make q1 the node label, e.g. "):0.5[&q1=0.6,...]" becomes ")0.6:0.5" s <- gsub("\\)(:[0-9.eE-]+)\\[&q1=([0-9.eE-]+),[^]]*\\]", ")\\2\\1", s) # step 3: delete the other annotations, which are the text inside square brackets s <- gsub("\\[[^]]*\\]", "", s) # read the cleaned text as a tree, and root it on P. montanus astral <- read.tree(text = s) astral <- root(astral, "mon_61", resolve.root = TRUE, edgelabel = TRUE) # the node labels are q1 values stored as text; turn them into numbers q1 <- suppressWarnings(as.numeric(astral$node.label)) # do the trees have the same topology? 0 means yes dist.topo(unroot(astral), unroot(ml)) dist.topo(unroot(astral), unroot(nj_tree)) summary(q1) ``` ```{r} #| label: astral-plot #| eval: false #| fig-height: 5.5 astral_l <- ladderize(astral) # node labels: q1 to two decimals (the root has no value, so it is left blank) q1_l <- suppressWarnings(as.numeric(astral_l$node.label)) q1_labels <- ifelse(is.na(q1_l), "", sprintf("%.2f", q1_l)) tip_colours <- species_cols[sp[astral_l$tip.label]] # use.edge.length = FALSE draws the tree without branch lengths (a cladogram) plot(astral_l, tip.color = tip_colours, cex = 0.9, use.edge.length = TRUE, underscore = TRUE, no.margin = TRUE) nodelabels(q1_labels, frame = "none", adj = c(1.2, -0.4), cex = 0.7) species_legend() ``` ![ASTRAL species tree from 10 kb window trees (authors' result), drawn without branch lengths. Tip colour is the species; node labels are q1, the fraction of gene-tree quartets supporting the branch.](figures/astral-plot.png){fig-alt="A species tree of 18 plants drawn as a cladogram, with q1 quartet support values at the nodes. The topology is the same as the ML tree. Support for most nodes is between about 0.4 and 0.65, with a higher value of 0.84 for the node uniting everything except the outgroups."} **Q5: The ML tree gave 100% bootstrap support to every species clade. All of the species tree estimation approaches match. The q1 values across nodes imply variation across the genome though. How can both be true?** A note on bootstrap support - it asks a very specific question: how much uncertainty is there due to sampling error? When alignments get long, like this one, it becomes less likely to get disagreements among bootstrap replicates. The coalescent tree is based on gene tree variation though. Most nodes have q1 between 0.4 and 0.65, against 1/3 for no signal. The *fruticosus* node (`dav_118`, `fru_106`) has q1 = 0.52, and the node joining *fruticosus* to *davidsonii* is among the weakest at 0.39. # 4. Tests of gene flow:Genome-wide D across the species The **ABBA-BABA test** asks whether two non-sister species share more derived alleles than incomplete lineage sorting alone would give. Take three species related as ((P1, P2), P3) and an outgroup O carrying the ancestral allele A. At an **ABBA** site, P2 and P3 share the derived allele B. At a **BABA** site, P1 and P3 share it. Incomplete lineage sorting produces both equally often, and gene flow between P2 and P3 produces extra ABBA. Significance comes from a **block jackknife** along the genome, reported as a Z-score, with |Z| > 3 the usual threshold. An external program, *Dsuite*, ran this test for every trio of the five species, with *P. montanus* as outgroup. It orders each trio so that D is positive. Groups in the Dsuite files are labelled with short species codes: `car` = *P. cardwellii*, `dav` = *P. davidsonii* (three plants), `fru` = *P. fruticosus* (two plants, including `dav_118`), `new` = *P. newberryi* and `rup` = *P. rupicola*. ```{r} #| label: dsuite-species # Dsuite's results: one row per trio of species dtab <- read.csv("data/dsuite_trios_species.csv", check.names = FALSE) # keep the columns we need, round them, and sort by Z-score (largest first) dtab_short <- data.frame(P1 = dtab$P1, P2 = dtab$P2, P3 = dtab$P3, D = round(dtab$Dstatistic, 3), Z = round(dtab$`Z-score`, 1), f4_ratio = round(dtab$`f4-ratio`, 3)) dtab_short[order(-dtab_short$Z), ] ``` **Q6: How many of the ten trios show significant gene flow? Which pair of species stands out most? Does the result for *P. newberryi* and *P. rupicola* surprise you?** Eight of the ten trios have Z > 3. The strongest is *newberryi* sharing more derived alleles with *davidsonii* than its sister *cardwellii* does (D = 0.088, Z = 18.8), the *davidsonii*–*newberryi* hybrid zone. The last row, car–new–rup, shows nothing (D = 0.003, Z = 0.9), which we take apart next. # 5. ABBA-BABA and fd by hand ## The equations With allele frequencies instead of single genomes, each site contributes a fractional ABBA and BABA. Let $p_1, p_2, p_3$ be the derived-allele frequencies in P1, P2 and P3. When the outgroup is ancestral ($p_O = 0$): $$\text{ABBA} = (1 - p_1)\,p_2\,p_3 \qquad \text{BABA} = p_1\,(1 - p_2)\,p_3$$ $$D = \frac{\sum \text{ABBA} - \sum \text{BABA}}{\sum \text{ABBA} + \sum \text{BABA}}$$ D is approximately 0 without gene flow and positive with P2–P3 gene flow. But in small windows it behaves erratically due to sampling error. **fd** (Martin et al. 2015) keeps the same numerator and divides by the ABBA − BABA you would get if P2 and P3 were randomly mating mixed. At each site, the donor frequency $p_D$ stands in for both P2 and P3 and is whichever has the higher derived-allele frequency: $$f_d = \frac{\sum \left[(1-p_1)\,p_2\,p_3 \;-\; p_1\,(1-p_2)\,p_3\right]} {\sum \left[(1-p_1)\,p_D\,p_D \;-\; p_1\,(1-p_D)\,p_D\right]}, \qquad p_D = \max(p_2,\,p_3)$$ fd only makes sense when D > 0. **fdM** (Malinsky et al. 2015) fixes that by using the same numerator but switching the denominator with the sign of D, so it is symmetric around zero: positive means P2–P3 gene flow, negative means P1–P3. $$f_{dM} = \frac{\sum \left[(1-p_1)\,p_2\,p_3 \;-\; p_1\,(1-p_2)\,p_3\right]}{\text{denominator}}, \qquad \text{denominator} = \begin{cases} \sum \left[(1-p_1)\,p_D\,p_D - p_1\,(1-p_D)\,p_D\right], \; p_D=\max(p_2,p_3) & D \ge 0 \\[6pt] -\sum \left[(1-p_D)\,p_2\,p_D - p_D\,(1-p_2)\,p_D\right], \; p_D=\max(p_1,p_3) & D < 0 \end{cases}$$ Where the outgroup is polymorphic ($p_O > 0$), Dsuite weights each term by $(1-p_O)$ and adds the same site counted with the roles of A and B swapped, weighted by $p_O$. That's what the code below does, and the extra terms vanish when $p_O = 0$. ## D by hand `data/trio_af_scaffold_1086.csv.gz` has one row per SNP on the largest scaffold (2.1 million). Each column is the derived-allele frequency in one group: | column | group | |---|---| | `p_car` | *P. cardwellii* (P1) | | `p_new` | *P. newberryi* (P2) | | `p_rup` | *P. rupicola* (P3) | | `p_new_noCL`, `p_rup_noCL` | the same two species **without the Castle Lake plants** (`NA` where no other plant was called) | | `p_out` | the outgroup, *P. montanus* | First, a worked example for a single SNP. Suppose P1 carries only the ancestral allele (frequency of the derived allele p1 = 0), and P2 and P3 both carry only the derived allele (p2 = 1, p3 = 1). The outgroup is ancestral. This site is a perfect ABBA site: ```{r} #| label: abba-example ex_p1 <- 0 ex_p2 <- 1 ex_p3 <- 1 (1 - ex_p1) * ex_p2 * ex_p3 # ABBA: P1 ancestral, P2 and P3 derived. 1 ex_p1 * (1 - ex_p2) * ex_p3 # BABA: P1 and P3 derived, P2 ancestral. 0 ``` Real SNPs have in-between frequencies (for example p2 = 0.5, because one of two plants carries the derived allele), so each SNP contributes a fraction of an ABBA and a fraction of a BABA. We add these up over all SNPs on the scaffold. The function below does this for every SNP at once. ```{r} #| label: abba-baba # read the table: one row per SNP, with the derived-allele frequency in each group af <- read.csv("data/trio_af_scaffold_1086.csv.gz") nrow(af) # number of SNPs head(af) # the first few rows # a function that gives the ABBA and BABA contribution of every SNP. # p1, p2, p3 and pO are columns of derived-allele frequencies in P1, P2, P3 and the outgroup. abba_baba <- function(p1, p2, p3, pO) { # the usual case, outgroup ancestral: P2 and P3 share the derived allele (ABBA), # or P1 and P3 share it (BABA). (1 - pO) is the chance the outgroup is ancestral. abba_main <- (1 - p1) * p2 * p3 * (1 - pO) baba_main <- p1 * (1 - p2) * p3 * (1 - pO) # if the outgroup is itself polymorphic (pO > 0), Dsuite also counts the site with the # roles of the two alleles swapped, weighted by pO. These terms are 0 when pO = 0. abba_swapped <- p1 * (1 - p2) * (1 - p3) * pO baba_swapped <- (1 - p1) * p2 * (1 - p3) * pO # return a table with one row per SNP data.frame(abba = abba_main + abba_swapped, baba = baba_main + baba_swapped) } # a function that turns the per-SNP table into D: add up ABBA and BABA over all SNPs # (na.rm = TRUE skips SNPs with missing values) calc_D <- function(sites) { abba_total <- sum(sites$abba, na.rm = TRUE) baba_total <- sum(sites$baba, na.rm = TRUE) (abba_total - baba_total) / (abba_total + baba_total) } # trio ((cardwellii, newberryi), rupicola), with all plants... all_samples <- abba_baba(af$p_car, af$p_new, af$p_rup, af$p_out) # ...and again without the two Castle Lake plants no_castle <- abba_baba(af$p_car, af$p_new_noCL, af$p_rup_noCL, af$p_out) c(all_samples = calc_D(all_samples), without_Castle_Lake = calc_D(no_castle)) ``` We can check our by-hand numbers against Dsuite on the same scaffold. Dsuite lists the trio as `new car rup`, because it always puts first whichever sister makes D positive, so its D has the opposite sign to yours. ```{r} #| label: abba-baba-check # Dsuite's own result for the same scaffold dsuite_check <- read.csv("data/dsuite_car_new_rup_scaffold_1086.csv", check.names = FALSE) dsuite_check[, c("samples", "P1", "P2", "P3", "Dstatistic", "Z-score", "ABBA", "BABA")] ``` **Q8: Removing the two hybrid zone plants flipped our D stat to be significant. Any thoughts on why including these could mask a signal of historical gene flow with this approach?** ## Some weaknesses about D in window-based tests To find *where* introgression happened, you might compute D in windows. In a low-diversity window ABBA + BABA is tiny, so D can become very noisy. Here we compute D and fd in 100 kb windows, without the Castle Lake plants. Since fd only works when D > 0, we swap P1 and P2 so the trio becomes ((*newberryi*, *cardwellii*), *rupicola*). We also keep only sites where the outgroup is homozygous. ```{r} #| label: windows-by-hand # step 1: keep the SNPs where the outgroup is ancestral (p_out is 0), and where both # newberryi and rupicola still have a called plant once Castle Lake is removed sites <- filter(af, p_out == 0, !is.na(p_new_noCL), !is.na(p_rup_noCL)) # step 2: name the three groups. P1 and P2 are swapped, so P1 = newberryi, P2 = cardwellii p1 <- sites$p_new_noCL p2 <- sites$p_car p3 <- sites$p_rup_noCL # step 3: the donor frequency pD. At each SNP it is the higher of p2 and p3. # pmax() takes the larger of two values at each position pD <- pmax(p2, p3) # step 4: three numbers for every SNP abba <- (1 - p1) * p2 * p3 # fraction of an ABBA site baba <- p1 * (1 - p2) * p3 # fraction of a BABA site fd_max <- (1 - p1) * pD * pD - p1 * (1 - pD) * pD # ABBA - BABA if P2 and P3 were completely mixed # step 5: a table of these, with the 100 kb window each SNP falls in # (a SNP at position 250,000 is in window 2, and so on) site_table <- data.frame(window = floor(sites$pos / 1e5), abba = abba, baba = baba, fd_max = fd_max) # step 6: add up the SNPs within each window win <- site_table |> group_by(window) |> summarise(n_sites = n(), abba_total = sum(abba), baba_total = sum(baba), fd_max_total = sum(fd_max)) # step 7: D and fd for each window. They share the same numerator (ABBA - BABA); # D divides by ABBA + BABA, fd divides by the maximum possible ABBA - BABA win$numerator <- win$abba_total - win$baba_total win$denom <- win$abba_total + win$baba_total win$D <- win$numerator / win$denom win$fd <- win$numerator / win$fd_max_total summary(win[, c("n_sites", "D", "fd")]) ``` ```{r} #| label: d-vs-fd #| eval: false #| fig-width: 10 #| fig-height: 4 # put D and fd in a single column called "value", with a column "statistic" saying which # is which, so ggplot can draw one panel for each win_long <- pivot_longer(win, cols = c(D, fd), names_to = "statistic", values_to = "value") # each dot is one 100 kb window; the x axis is the amount of ABBA + BABA in the window p <- ggplot(win_long, aes(x=statistic, y=value, fill=statistic)) + geom_violin() + geom_jitter(height = 0, width = 0.2, aes(color = n_sites)) + scale_color_viridis_c(option = "plasma") print(p) ``` ![Distribution of D and fd in 100 kb windows on scaffold 1086 for ((*newberryi*, *cardwellii*), *rupicola*), without the Castle Lake plants.](figures/Distribution_D_fd.png){fig-alt="Violin plots showing the distribution of D and fd in 100 kb windows, with points colored by the number of sites."} The main takeaway is that people are always playing with this type of statistic to make it more interpretable and robust. In this case, fd is about small genomic windows. But if you grasp the concept of D, it should translate to other scenarios. # 8. Flower colour: one origin or two? Both hummingbird-pollinated species, *P. newberryi* and *P. rupicola*, have magenta flowers, while their bee-pollinated relatives are blue-violet. The blue colour comes from delphinidin-type anthocyanin pigments, which need the enzyme **F3′5′H**. Magenta flowers contain only the simpler pink cyanidin pigments, so they are consistent with a loss of F3′5′H function. The paper will show some excepetions. In our trees *P. newberryi* is sister to *P. cardwellii*. *P. rupicola* sits outside them, so the two magenta species are not sisters. The paper considered three hhypotheses for the loss of F3′5′H: 1. **Two independent losses**, one in each species. 2. **One loss in the common ancestor**, followed by reversals in the species that make blue pigments again. A broken gene is hard to repair, so this seems unlikely. 3. **One loss shared between *newberryi* and *rupicola***, by incomplete lineage sorting or by introgression (adaptive introgression of the magenta allele). ## Results rejecting hypothesis three Could magenta flowers have a single origin shared between *P. newberryi* and *P. rupicola* due to gene flow? We can test this by looking at fdM around the **F3′5′H** gene. If the magenta allele was shared by introgression, we would expect a peak at the **F3′5′H** gene. The values are Dsuite fdM for (*cardwellii*, *newberryi*, *rupicola*) in 500-SNP windows sliding by 250 SNPs, for all samples and without the two Castle Lake plants. The dashed line is the 95th percentile of windows on the scaffold, which we call "unusually high". Following the paper, we look at the 400 kb around the gene, and also rank that 400 kb against every other 400 kb block on the scaffold. ```{r} #| label: f3p5h # where the gene is: its scaffold, start and end, and its midpoint locus <- read.csv("data/f3p5h_locus.csv") gene_mid <- (locus$start + locus$end) / 2 # fdM in windows for the trio (car, new, rup), for both sets of samples nr <- read.csv("data/fdm_car_new_rup_500SNP.csv.gz") # keep the windows on the scaffold that holds the gene, and add each window's midpoint scaf <- filter(nr, scaffold == locus$scaffold) scaf$mid <- (scaf$start + scaf$end) / 2 # the paper looked at 400 kb around the gene block <- 4e5 # the start of every 400 kb block along the scaffold, sliding along 50 kb at a time block_starts <- seq(min(scaf$mid), max(scaf$mid) - block, by = 5e4) # the two sets of samples: "all" and "no Castle Lake" sample_sets <- unique(scaf$samples) # an empty table that we fill in, one row per sample set results <- data.frame() for (set in sample_sets) { # the windows for this sample set d <- scaf[scaf$samples == set, ] # fdM in the windows within 200 kb either side of the gene near_gene <- d$fdM[abs(d$mid - gene_mid) <= block / 2] # the mean fdM of every 400 kb block on the scaffold block_means <- numeric(length(block_starts)) for (b in seq_along(block_starts)) { in_block <- d$mid >= block_starts[b] & d$mid < block_starts[b] + block block_means[b] <- mean(d$fdM[in_block]) } # how does the gene's neighbourhood compare with the rest of the scaffold? one_row <- data.frame( samples = set, highest_window_near_gene = max(near_gene), scaffold_95th_percentile_window = unname(quantile(d$fdM, 0.95)), mean_near_gene = mean(near_gene), percentile_among_400kb_blocks = mean(block_means <= mean(near_gene))) results <- rbind(results, one_row) } # round the numbers (everything except the first column) for display results[, -1] <- signif(results[, -1], 2) results ``` ```{r} #| label: f3p5h-plot #| eval: false #| fig-width: 9 #| fig-height: 4 # the windows within 3 Mb of the gene, which is the part we draw near3 <- scaf[abs(scaf$mid - gene_mid) < 3e6, ] # the 95th percentile of fdM on the whole scaffold, for each set of samples (the dashed line) q95 <- aggregate(fdM ~ samples, data = scaf, FUN = quantile, probs = 0.95) ggplot(near3, aes(mid / 1e6, fdM)) + # shaded band: the 400 kb around the gene annotate("rect", xmin = (gene_mid - block / 2) / 1e6, xmax = (gene_mid + block / 2) / 1e6, ymin = -Inf, ymax = Inf, fill = "#7570b3", alpha = 0.15) + geom_hline(yintercept = 0, colour = "grey70") + # dashed line: the 95th percentile of windows on the scaffold geom_hline(data = q95, aes(yintercept = fdM), linetype = "dashed", colour = "grey40") + # purple line: the gene geom_vline(xintercept = gene_mid / 1e6, colour = "#7570b3", linewidth = 0.8) + geom_line(colour = "grey30") + facet_wrap(~ samples) + labs(x = paste(locus$scaffold, "position (Mb)"), y = expression(f[dM]~"(car, new, rup)")) ``` ![fdM for (*cardwellii*, *newberryi*, *rupicola*) in 500-SNP windows within 3 Mb of *F3′5′H* (purple line; shaded: the 400 kb around it). Dashed: the 95th percentile of windows on the scaffold.](figures/f3p5h-plot.png){fig-alt="Two line plots of fdM along six megabases of scaffold 1086 around the F3'5'H gene, for all samples and without the Castle Lake plants. Within the shaded 400 kb region around the gene fdM stays well below the dashed 95th-percentile line. A tall spike lies about 0.6 megabases to the right of the gene."} None of the values really stand out and the 400kb around the gene are not differentiated from the background. Probably safe to conclude that there is no support for a shared allele via introgression. ## Results rejecting hypothesis two If the magenta allele moved between the two species, *newberryi* and *rupicola* should be close relatives in a tree of *F3′5′H* itself. We will just look at the gene tree for *F3′5′H*. Following the paper, we collapse branches with low support, since a single gene of about 1.5 kb has few informative sites. ```{r} #| label: f3p5h-tree # read the gene tree and root it on the two outgroups (which are each other's closest relatives) gt <- read.tree("data/trees/f3p5h_genetree.treefile") gt <- root(gt, c("lya_44", "mon_61"), resolve.root = TRUE, edgelabel = TRUE) # the sample IDs of the newberryi plants and of the rupicola plants new <- names(sp)[sp == "P. newberryi"] rup <- names(sp)[sp == "P. rupicola"] # is each species a clade (monophyletic) in this gene tree? newberryi_monophyletic <- is.monophyletic(gt, new) rupicola_monophyletic <- is.monophyletic(gt, rup) # the smallest clade that holds all eight newberryi and rupicola plants: how many plants is it? mrca_node <- getMRCA(gt, c(new, rup)) smallest_clade <- extract.clade(gt, mrca_node) c(newberryi_monophyletic = newberryi_monophyletic, rupicola_monophyletic = rupicola_monophyletic, tips_in_smallest_clade_holding_all_new_and_rup = length(smallest_clade$tip.label)) # collapse branches with no SH-aLRT support, as in the paper. # Node labels look like "85.2/85" (SH-aLRT/UFBoot): take the number before the "/" shalrt <- suppressWarnings(as.numeric(sub("/.*", "", gt$node.label))) # internal nodes are numbered after the tips, so add the number of tips no_support_nodes <- which(shalrt == 0) + Ntip(gt) # gt$edge lists the branches; its second column is the node at the end of each branch. # Set the length of every branch leading to an unsupported node to 0... gt_collapsed <- gt is_no_support_branch <- gt$edge[, 2] %in% no_support_nodes gt_collapsed$edge.length[is_no_support_branch] <- 0 # ...and merge the zero-length branches into multifurcations gt_collapsed <- di2multi(gt_collapsed, tol = 1e-8) ``` ```{r} #| label: f3p5h-tree-plot #| eval: false #| fig-height: 5 gt_l <- ladderize(gt_collapsed) # node labels, leaving the one called "Root" blank node_labels <- gt_l$node.label node_labels[node_labels == "Root"] <- "" tip_colours <- species_cols[sp[gt_l$tip.label]] plot(gt_l, tip.color = tip_colours, cex = 0.9, underscore = TRUE, no.margin = TRUE) nodelabels(node_labels, frame = "none", adj = c(1.1, -0.4), cex = 0.7) species_legend() add.scale.bar() ``` ![ML gene tree of *F3′5′H*, rooted on the outgroups, with low-support branches collapsed. Tip colour is the species; node labels are SH-aLRT/UFBoot.](figures/f3p5h-tree-plot.png){fig-alt="A gene tree of 18 plants with most branches collapsed into multifurcations. The four newberryi plants form a group, but the rupicola plants are scattered: three sit with cardwellii and a fruticosus plant, and one sits near the newberryi group."} **Q9: Do the *P. newberryi* and *P. rupicola* plants group together in the gene tree? Why would hypothesis 2 predict this? Is something different about the Castle Lake plant `rup_86`?** ## Results supporting hypothesis 1 That will be covered during the discussion! # 9. Takeaways 1. **Named species are hypotheses the genomes can test.** All five species form clades in the NJ, ML and ASTRAL trees, but the *davidsonii*–*fruticosus* relationship is the weakest (70% bootstrap, q1 = 0.39). 2. **Species relationships can have high discordance.** A high proportion of gene trees reflected relationships different from the species tree. This was likely due to historical gene flow between present-day species. 3. **Hybrid-zone samples can affect tests of introgression.** Two Castle Lake plants turned D = −0.037 (Z = 12.6) into D = 0.003. 4. **Trees and introgression tests can be used to evaluate hypotheses of gene family evolution and speciation.** ## References - Stone BW, Wessinger CA (2024) Ecological diversification in an adaptive radiation of plants: the role of de novo mutation and introgression. *Mol Biol Evol* 41:msae007. - Malinsky M, Challis RJ, Tyers AM, et al. (2015) Genomic islands of speciation separate cichlid ecomorphs in an East African crater lake. *Science* 350:1493–1498. - Martin SH, Davey JW, Jiggins CD (2015) Evaluating the use of ABBA–BABA statistics to locate introgressed loci. *Mol Biol Evol* 32:244–257. - Malinsky M, Matschiner M, Svardal H (2021) Dsuite — fast D-statistics and related admixture evidence from VCF files. *Mol Ecol Resour* 21:584–595. - Cruickshank TE, Hahn MW (2014) Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow. *Mol Ecol* 23:3133–3157.