Is selfing a dead end oin Solanaceae?

0. Technical details

Similar to last time, create a directory for your work and navigate to it in R with setwd("path/to/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: ape, diversitree, ggplot2 and patchwork. diversitree is the new one, so please install it before class with install.packages("diversitree").
We will fit every BiSSE model and run a short MCMC ourselves. Some of today’s analyses take too long for class (long MCMC chains and the hidden-state models), so those outputs have been provided.
This is a computer activity and a paper discussion in one. The main text of Goldberg et al. (2010) is very short, so we will just have quick look at it while some analyses run and discuss.

1. The dataset and the question

Most flowering plants are hermaphrodites, yet many can’t fertilise themselves. In Solanaceae the style makes an S-RNase that destroys the RNA of growing pollen tubes, and the pollen carries proteins that disable every S-RNase except the one matching its own S allele. So a pollen grain is rejected by any style carrying its allele, including its own parent’s.

Three things about this system:

SI is the ancestral state (Ioerger et al. 1990).
SI is easy to lose through loss-of-function mutations that kill S-RNase activity in the pistil (Broz et al. 2021).
Reversals from SC to SI are unlikely. One the function is lost, there is no selective pressure against more mutations (Igić et al. 2006).

About 41% of Solanaceae retain their ancestral SI state. Goldberg et al. (2010) hypothesized that while being self-compatible could be advantageous in the right context, SI lineages diversify fast enough to replace the lineages that lose it. Stebbins (1957) suggested that selfing was an evolutionary dead end.

Here, we have the dated tree and the SI or SC state of each species. We will reproduce their result with BiSSE, and explore some of the limitations.

Two caveats about this data:

The tree has 356 of about 2,300 accepted species. We correct for that with sampling fractions.
BiSSE needs every species in one of two states. We consider any species with some SI as SI (some species have both), and everything else counts as SC. The SC group includes 25 dioecious species, which have no SI but cannot self either. So, consider “SC” here as “not SI”.

2. Reading the data into R

library(ape)           # reading and handling trees
library(diversitree)   # BiSSE
library(ggplot2)       # plots
library(patchwork)     # putting several plots side by side
theme_set(theme_bw(base_size = 11))

# one colour per breeding system, used in every figure
sys_cols <- c("SI" = "#0072B2", "SI+SC" = "#56B4E9", "SC" = "#D55E00")
two_cols <- sys_cols[c("SI", "SC")]    # the two states BiSSE uses

# the biggest genera, which we name in the figures
big_genera <- c("Solanum", "Nicotiana", "Lycium", "Capsicum", "Physalis", "Nolana")

The tree and the traits are two plain files. The tree is the dated tree from the paper, already ultrametric (every tip is the same distance from the root), which diversitree insists on. The trait file has one row per species, with 0 for SI and 1 for SC.

tree <- read.tree("data/Solanaceae_timetree.tre")
traits <- read.table("data/trait_data.txt", header = TRUE)
head(traits)

tree                                    # 356 tips
is.ultrametric(tree)                    # TRUE, so diversitree will accept it

# BiSSE wants a 0 or a 1 for every species, labelled with its name and in the same
# order as the tips of the tree
states <- traits$state
names(states) <- traits$species
states <- states[tree$tip.label]
table(states)

You should have 135 SI and 221 SC species. We also have a table with the three-way summary: 0 is SI, 1 is SI+SC (the species is polymorphic) and 2 is SC. The two-state state01 column is the same trait as above.

three <- read.table("data/trait_data_3state.txt", header = TRUE)

# the label for each code, in order from 0 to 2. R counts from 1, so we add 1 to each
# code: code 0 picks the first name, code 1 the second, and so on.
system_names <- c("SI", "SI+SC", "SC")
three$system <- system_names[three$state012 + 1]
table(three$system)

# the two-state trait is the same thing with SI+SC folded into SI
all(three$state01 == as.integer(three$state012 == 2))
all(three$state01 == traits$state)

# a lookup for later: give it species names and it returns their breeding systems
system_of <- three$system
names(system_of) <- three$species

Now let’s see how the two states are spread around the tree.

# ladderize() sorts clades by size, which only changes how the tree is drawn
tree_to_draw <- ladderize(tree)
tips <- tree_to_draw$tip.label

# trait.plot() wants the states in a data frame, with species as the row names
ring <- data.frame(system = states[tips], row.names = tips)

# name only the big genera; NA means no label
genus <- sub("_.*", "", tips)           # "Solanum_tuberosum" becomes "Solanum"
genus_label <- ifelse(genus %in% big_genera, genus, NA)

# type = "f" makes it a fan, and w sets the width of the ring
trait.plot(tree_to_draw, ring, cols = list(system = unname(two_cols)),
           class = genus_label, type = "f", w = 1/30, margin = 0.3,
           lab = "", legend = FALSE, cex.lab = 0.9)
legend("topleft", legend = names(two_cols), fill = two_cols, border = NA, bty = "n",
       cex = 0.9)

A circular tree of 356 species with a coloured ring around the tips. Blue self-incompatible and vermillion self-compatible tips are interspersed in most of the tree. Solanum has some long vermillion blocks, Nicotiana is nearly all vermillion, and Nolana is mostly blue. — The 356 species, with each species’ state in the ring around the tips. The root is 36 million years old.

SI and SC are interspersed throughout, which is expected under a multiple loss scenario.

Q1: Does it look like some clades have more SI than others?

3. Fitting BiSSE

BiSSE (Maddison et al. 2007) gives each state its own speciation rate ($\lambda$) and extinction rate ($\mu$), plus rates of change between the states. With SI as state 0 and SC as state 1:

`diversitree`	paper	meaning
`lambda0`, `lambda1`	$\lambda_I$, $\lambda_C$	speciation in SI and SC lineages
`mu0`, `mu1`	$\mu_I$, $\mu_C$	extinction in SI and SC lineages
`q01`	$q_{IC}$	loss of SI
`q10`	$q_{CI}$	regain of SI

make.bisse() returns a likelihood function, and constrain() makes a simpler model from it. The sampling fractions come from the paper’s supplement: 16.2% of SI species and 15.0% of SC species are in the tree.

f <- c(0.162, 0.150)                    # SI first, then SC

# lik is a function: give it the six rates and it returns the likelihood of our data
lik <- make.bisse(tree, states, sampling.f = f)
argnames(lik)                           # the six rates it asks for

# rough values for the search to start from
start <- starting.point.bisse(tree)

Is one set of rates enough for both states? If not, is it speciation that differs, extinction, or both? And should SI be allowed to come back? We can write each idea as a model and fit them one at a time, starting with the simplest. find.mle() climbs the likelihood from our starting values until it can’t improve. First, a small helper to print what we want to see from each fit.

# print the fitted rates, the log likelihood and the AIC (lower AIC is better)
show_fit <- function(fit) {
  print(round(fit$par, 3))
  c(lnL = round(fit$lnLik, 2), AIC = round(AIC(fit), 1))
}

Model 1: same rates. SI and SC speciate and go extinct at the same rates, and SI is lost but never regained. constrain() ties one rate to another: lambda1 ~ lambda0 sets SC’s speciation rate equal to SI’s, so there is one rate less to estimate, and q10 ~ 0 stops SI from coming back.

same_model <- constrain(lik, lambda1 ~ lambda0, mu1 ~ mu0, q10 ~ 0)
argnames(same_model)                    # only three rates are left to estimate

# start[argnames(same_model)] picks out the starting values for just those rates
same_rates <- find.mle(same_model, start[argnames(same_model)])
show_fit(same_rates)

A single speciation rate and a single extinction rate, both near 4/Ma and almost equal. Relaxing the contraints might improve the AIC.

Model 2: speciation differs. SC gets its own speciation rate, but extinction is shared.

speciation_model <- constrain(lik, mu1 ~ mu0, q10 ~ 0)
speciation_differs <- find.mle(speciation_model, start[argnames(speciation_model)])
show_fit(speciation_differs)

Model 3: extinction differs. The reverse: SC gets its own extinction rate and speciation is shared.

extinction_model <- constrain(lik, lambda1 ~ lambda0, q10 ~ 0)
extinction_differs <- find.mle(extinction_model, start[argnames(extinction_model)])
show_fit(extinction_differs)

Model 4: both differ. This is the paper’s model: both rates can differ, and SI never returns. Compare its rates with the first three.

both_model <- constrain(lik, q10 ~ 0)
both_differ <- find.mle(both_model, start[argnames(both_model)])
show_fit(both_differ)

Model 5: SI can come back. The full model with nothing constrained. We do not think this is plausible, but will include it to hopefully make a point about hypothesis testing with models.

si_can_return <- find.mle(lik, start[argnames(lik)])
show_fit(si_can_return)

Now let’s make our own table to compare all five. The delta_AIC is the different between each model AIC and the best AIC. Can you fill out the following table from your results?

Model	Parameters	lnL	AIC	delta AIC
Same rates
Speciation differs
Extinction differs
Both differ
SI can return

Letting rates differ between SI and SC is the best model, but we do not consider this to be plausible. We will ignore that issue for now and focus on the rates estimated by the “both_differ” model.

You can access the parameter estimates for this model by typing

both_differ$par

Can you use these estimates to calcualte the net diversification rates for SI ($r_0$) and SC lineages ($r_1$)? They should be pretty close to the paper’s reported values.

Q2: Which diversification rate is higher, SI (state 0) or SC (state 1)? Did you get a negative number for one of them; what does that mean?

4. Bayesian analyses and testing the dead-end SC hypothesis.

The paper and most studies use a Bayesian analysis, sampling many parameter combinations to assess uncertainty in the estiamtes. Remember, we suggested there were some non-identifiability issues between $\lambda$ and $\mu$. An MCMC wanders through the rates, visiting each combination in proportion to how well it fits the data and how plausible the prior says it is. Bayesian models need prior distribtions on all parameters. The prior here is the paper’s: an exponential with a mean of 1/0.3 on every rate. We will try to talk about priors while your analysis runs Each iteraction of the MCMC calculates the likelihood again, so the computation can take a while. We will run a short one here and compare it with four long chains that I ran ahead of time.

Run this chunk and let it go. It takes about 12 minutes on a laptop, and prints a line every 500 steps, so carry on with the reading below while you wait.

# the paper's model again: SI is lost but never regained
one_way <- constrain(lik, q10 ~ 0)
prior <- make.prior.exponential(0.3)

# Start the chain at the maximum likelihood rates. w is how wide a step the sampler
# tries first, and set.seed() makes your chain repeatable.
set.seed(1)
short <- mcmc(one_way, both_differ$par, nsteps = 5000, w = 1, prior = prior,
              print.every = 500)

write.csv(short, "short_chain.csv", row.names = FALSE)    # save it in case R crashes

Each row is one step: the step number i, the five rates and p, the log of likelihood times prior. Now the trace, which is each rate plotted against step. A chain that has found the right region looks like a fuzzy band that stays level.

rate_names <- c("lambda0", "lambda1", "mu0", "mu1", "q01")

par(mfrow = c(2, 3))                    # two rows and three columns of plots
for (rate in rate_names) {
  plot(short$i, short[[rate]], type = "l", xlab = "step", ylab = rate)
}
plot(short$i, short$p, type = "l", xlab = "step", ylab = "log posterior")
par(mfrow = c(1, 1))                    # back to one plot at a time

The log posterior levels off almost at once, but the rates wander slowly up and down, and λ and μ move together: a pair of larger values fits about as well as a pair of smaller ones, which we will meet again in section 6. The chain starts at the maximum likelihood estimate, which is not a typical draw from the posterior, so we drop the first 500 steps as burn-in. We then work out the same quantities as the table in section 4 and compare with the long chains.

Also we can calculate the proportion of SI lineages expected by the model ($p_{SI}$) from the rates using the formula $p_{SI} = 1 - q_{IC} / (r_I - r_C)$, where $r_I = \lambda_I - \mu_I$ and $r_C = \lambda_C - \mu_C$.

# The chains hold only the five rates in the model, so add the quantities built from
# them. Each line does the sum for every row at once.
add_derived <- function(x) {
  x$r_I <- x$lambda_I - x$mu_I
  x$r_C <- x$lambda_C - x$mu_C
  x$p_SI <- 1 - x$q_IC / (x$r_I - x$r_C)
  x
}

# our short chain, with the burn-in dropped and the columns named as in the paper
mine <- short[short$i > 500, ]
names(mine)[names(mine) %in% rate_names] <- c("lambda_I", "lambda_C", "mu_I", "mu_C", "q_IC")
mine <- add_derived(mine)

# the long chains: 4 chains, with 500 samples kept from each
post <- read.csv("data/bisse_mcmc.csv")      # one row per sample from the chains
post <- add_derived(post)
table(post$chain)

# the same quantities as the table above, in the same order
my_rates <- mine[names(paper)]
post_rates <- post[names(paper)]

post_table <- data.frame(short_mean = sapply(my_rates, mean),
                         long_mean = sapply(post_rates, mean),
                         Table_S1 = paper,
                         short_sd = sapply(my_rates, sd),
                         long_sd = sapply(post_rates, sd))
round(post_table, 3)

Now let’s put the two side by side. The grey density is the longer chains, and the colour is your short one.

# one density for the long chains and one for our chain, for a single quantity
compare_chains <- function(column, x_label) {
  both <- data.frame(value = c(post[[column]], mine[[column]]),
                     chain = rep(c("long chains", "our short chain"),
                                 times = c(nrow(post), nrow(mine))))
  ggplot(both, aes(value, fill = chain)) +
    geom_density(alpha = 0.5, colour = NA) +
    scale_fill_manual(values = c("long chains" = "grey50",
                                 "our short chain" = "#D55E00")) +
    labs(x = x_label, y = "posterior density", fill = NULL)
}

compare_chains("lambda_I", "SI speciation") +
  compare_chains("mu_C", "SC extinction") +
  compare_chains("r_C", "SC net diversification") +
  compare_chains("q_IC", "loss of SI") +
  plot_layout(nrow = 1, guides = "collect")

Why do the speciation and extinction traces wander so slowly? The sampler moves one rate at a time, and some rates are tied tightly together. Let’s plot pairs of parameters against each other, with the long chains in grey and our short chain in colour.

# one point per sample for a pair of parameters
scatter_pair <- function(x_column, y_column, x_label, y_label) {
  ggplot() +
    geom_point(data = post, aes(.data[[x_column]], .data[[y_column]]),
               size = 0.4, alpha = 0.3, colour = "grey60") +
    geom_point(data = mine, aes(.data[[x_column]], .data[[y_column]]),
               size = 0.4, alpha = 0.4, colour = "#D55E00") +
    labs(x = x_label, y = y_label)
}

scatter_pair("lambda_I", "mu_I", "SI speciation", "SI extinction") +
  scatter_pair("lambda_C", "mu_C", "SC speciation", "SC extinction") +
  scatter_pair("r_I", "q_IC", "SI net diversification", "loss of SI") +
  scatter_pair("r_I", "r_C", "SI net diversification", "SC net diversification") +
  plot_layout(nrow = 1)

With the seed above, our 5,000 steps land close to the long chains: the means for the net diversification rates and $p_{SI}$ agree to within about 0.02, and the speciation and extinction means are off by up to 0.15. The scatterplots show why those two are less reliable. Speciation and extinction are highly correlated, and since the chain changes one rate at a time, that parameter space is explored slowly. Note the last panel is a little different, more of a cloud. This is positive as it signals there is some independence between the net diversifiction rates.

If your chain crashed or you ran out of time, short <- read.csv("data/short_chain_example.csv") reads in the one used for this text, and you can carry on from the trace plot.

Means and standard deviations from the long chains match Table S1 closely, so we have reproduced the paper. From here on we use the long chains, in post. To see how the next chunks do with your own chain, replace post by mine instead.

Now, we can get to the core hypothesis test. Do SI lineages diversify faster than SC lineages?

# Draw the posterior of one rate for SI and SC on the same axes. ggplot wants all of
# the values in one column, with a second column saying which state each one is.
two_densities <- function(si_column, sc_column, x_label) {
  long <- data.frame(rate = c(post[[si_column]], post[[sc_column]]),
                     state = rep(c("SI", "SC"), each = nrow(post)))
  long$state <- factor(long$state, levels = c("SI", "SC"))    # SI first in the legend

  ggplot(long, aes(rate, fill = state)) +
    geom_density(alpha = 0.6, colour = NA) +
    scale_fill_manual(values = two_cols) +
    labs(x = x_label, y = "posterior density", fill = NULL)
}

speciation_plot <- two_densities("lambda_I", "lambda_C", "speciation rate (per Myr)")
extinction_plot <- two_densities("mu_I", "mu_C", "extinction rate (per Myr)")
net_div_plot <- two_densities("r_I", "r_C", "net diversification (per Myr)") +
  geom_vline(xintercept = 0, linetype = 2)

# one point per sample, to show how SC's two rates rise and fall together
sc_rates_plot <- ggplot(post, aes(lambda_C, mu_C)) +
  geom_point(size = 0.3, alpha = 0.3, colour = two_cols[["SC"]]) +
  geom_abline(linetype = 2, colour = "grey50") +
  labs(x = "SC speciation", y = "SC extinction")

# patchwork lets us add plots together
speciation_plot + extinction_plot + net_div_plot + sc_rates_plot +
  plot_layout(nrow = 1, guides = "collect")

Four panels. In the first three, blue is SI and vermillion SC: SC speciation is higher than SI, SC extinction is higher than SI, and SI net diversification is positive while SC is negative, with no overlap. The fourth panel plots SC extinction against SC speciation as a tight diagonal band that sits just above the one-to-one line. — Posterior distributions, as in the paper’s Figure 1. Speciation and extinction are uncertain, and net diversification is much tighter. The right panel shows why: λ_C and μ_C rise and fall together. Points above the dashed line have λ < μ.

So, SI species grow at $r_{I}$ and evolve into SC at $q_{IC}$. So SI can only hold on if it out-diversifies SC by more than the rate it is lost. When that is true, the propotion of SI species can then be expressed as:

\[p_{SI} = 1 - \frac{q_{IC}}{r_I - r_C}\]

The posterior mean for $p_{SI}$ from the post object is 35% of species. The observed proportion is 41% across Solanaceae. Pretty close.

Q3: Do these model estimates explain the persistence of SI lineages despite its susceptibility to loss-of-function?

Stebbins suggested there are many SC species today only because SI lineages keep producing new ones.

5. A trait that does nothing

Five years after this paper, the first author published a study suggesting caution in interpreting BiSSE results (Rabosky & Goldberg 2015). They simulated traits with no effect on diversification on real trees, and BiSSE often reported an effect anyway. Real trees are uneven for many reasons, so any trait that happens to sit in a fast clade can get the credit. Could that be happening here?

We can run the same test on this tree, using the model where SI can only be lost. I simulated 300 traits that can only be lost, at a rate that leaves about as many tips in the starting state as there are SI species today. These traits have no connection to diversification. For each one I fit BiSSE as in section 3 and we can retain the difference in net diversification between the two states ($r_{I} - r_{C}$) as our summary statistic. We can check how different our observed trait is from the simulated distribution to find out if the trait-dependent rates are surprising.

neutral <- read.csv("data/neutral_oneway_fits.csv")           # the 300 simulated traits
real <- read.csv("data/neutral_oneway_observed.csv")          # the real trait

neutral$r_gap <- neutral$r_I - neutral$r_C
real$r_gap <- real$r_I - real$r_C

ggplot(neutral, aes(r_gap)) +
  geom_histogram(binwidth = 0.05, boundary = 0, fill = "grey60") +
  geom_vline(xintercept = real$r_gap, colour = two_cols[["SC"]], linewidth = 1) +
  labs(x = "net diversification of SI minus SC (per Ma)", y = "simulated traits")

A histogram of the gap in net diversification for 300 simulated traits, all between about −0.1 and 0.35 and centred near 0.1. The vermillion line for the real trait, at 0.74, lies far to the right of all of them. — The gap in net diversification between the two states for 300 simulated traits with no effect on diversification, which can only lose the first state. The vermillion line is the real trait.

The real trait associated with diversification is very different from a trait that is not.

Q4: There is a real difference in diversification between states. Does that show that SI itself changes diversification? Can you think of other factors that would lead to transitions from SI to SC AND differences in diversification rates?

6. Hidden states

BiSSE imposes a limited hypothesis test. - $H_O$: rates do not vary over the tree - $H_A$: rates do vary over the tree, and that depends on the trait

Beaulieu & O’Meara (2016) aimed to correct this limitation by letting rates vary for another reason, some unknown trait. In the CID models rates depend only on a hidden trait with 2 or 4 classes, and in HiSSE they depend on the observed (e.g. SI vs SC) and a hidden trait.

The bigger models take about half an hour each to fit, so I ran them ahead of time with the paper’s assumptions. One fit looks like this, shown for reading and not for running.

library(hisse)
dat <- data.frame(species = tree$tip.label, state = states)

# the allowed changes: between SI and SC, and between two hidden classes
transitions <- TransMatMakerHiSSE(hidden.traits = 1, make.null = TRUE)
loss_only <- ParDrop(transitions, 1)         # rate 1 is SC -> SI, which we remove

# Rates are given in the order 0A, 1A, 0B, 1B, and equal numbers share a rate. So this
# is CID-2: rates differ between hidden classes A and B, but not between SI and SC.
cid2_rates <- c(1, 1, 2, 2)

hisse(tree, dat, f = c(0.162, 0.150), hidden.states = TRUE, trans.rate = loss_only,
      turnover = cid2_rates, eps = cid2_rates)

hisse_models <- read.csv("data/hisse_models.csv")

# leave out the versions that share one extinction fraction across states
keep <- !grepl("_eps", hisse_models$model)
hisse_models <- hisse_models[keep, ]

# k is the number of parameters, and dAICc is how far each model is behind the best
hisse_table <- hisse_models[, c("model", "k", "lnL", "dAICc")]
hisse_table$lnL <- round(hisse_table$lnL, 1)
hisse_table$dAICc <- round(hisse_table$dAICc, 1)
print(hisse_table, row.names = FALSE)

BiSSE is 246 AICc units behind the best model, and it is even worse than CID-2, where we have two unknown traits (but transition between both states are possible unlike our constrained BiSSE model). The fair comparison is the top two rows. HiSSE and CID-4 have the same number of parameters, and HiSSE wins by 12 AICc units.

Here are the rates in the best model.

hisse_rates <- read.csv("data/hisse_rates.csv")

# keep just the rates from the HiSSE model
best <- hisse_rates[hisse_rates$model == "HiSSE", ]

rate_columns <- c("lambda", "mu", "r")
best[rate_columns] <- round(best[rate_columns], 2)
print(best, row.names = FALSE)

Compare SI with SC within each hidden class. SI has the higher net diversification in both, 0.73 against −0.14 in class A and 0.44 against −0.01 in class B. The core biological result is that SI lineages have higher net diversification compared to SC survives the hidden-class correction.

Can we find the lineages that belong to the strange fast turnover class?

tip_rates <- read.csv("data/hisse_tip_rates.csv")      # one row per species

# state 0 is SI and state 1 is SC, as before
tip_rates$system <- ifelse(tip_rates$state == 0, "SI", "SC")
tip_rates$system <- factor(tip_rates$system, levels = c("SI", "SC"))

# Lump the small genera together. rev() puts the first genus at the top of the plot.
tip_rates$group <- ifelse(tip_rates$genus %in% big_genera, tip_rates$genus, "others")
tip_rates$group <- factor(tip_rates$group, levels = rev(c(big_genera, "others")))

# one row of dots per genus, with the x axis added separately for each panel below
dot_plot <- ggplot(tip_rates, aes(y = group, colour = system)) +
  geom_jitter(height = 0.2, width = 0, size = 1, alpha = 0.6) +
  scale_colour_manual(values = two_cols) +
  labs(y = NULL, colour = NULL)

net_div_panel <- dot_plot + aes(x = net.div) +
  geom_vline(xintercept = 0, colour = "grey60") +
  labs(x = "net diversification (per Myr)")

net_div_panel

The two hidden classes roughly mean “Solanum” and “everything else”: median turnover is 36 per Myr in Solanum and about 2.7 in most other groups. SI and SC separate within every genus, but part of that is by construction, since each tip’s rate depends on its state.

Q5: The fast hidden class is mostly Solanum. What do you think the hidden trait could be? What could you do to test that selection acts on SI itself and not an associated trait?

7. Takeaways

SC is a capacity, not a mating system. The tree records the loss of SI, and selfing is a separate component.
Trees can be used to estimate diversification rates. There was high uncertainty in the speciation and extinction rates, but net diversification was quite reliable, even with a short MCMC chain.
A trait with no effect can’t produce the observed rates under BiSSE. When SI can only be lost, none of 300 simulated random traits were similar to our observed rates.
The effect survives hidden states, but there is more to learn. SI out-diversifies SC in both hidden classes, but something else we do not know is important to Solanum.

--- title: "Is selfing a dead end oin Solanaceae?" format: html: toc: true toc-depth: 2 code-fold: false code-tools: true embed-resources: true fig-width: 7 fig-height: 5 execute: eval: false --- # 0. Technical details - Similar to last time, create a directory for your work and navigate to it in R with `setwd("path/to/your/directory")` - All of the inputs are in `data/`, which you can download from [drive](https://drive.google.com/drive/folders/1o7PvjYZ2LaFoLRj7wdFsCUYw3sxcmabP?usp=sharing) and place in your working directory. - A few R packages will be used: `ape`, `diversitree`, `ggplot2` and `patchwork`. `diversitree` is the new one, so please install it before class with `install.packages("diversitree")`. - We will fit every BiSSE model and run a short MCMC ourselves. Some of today's analyses take too long for class (long MCMC chains and the hidden-state models), so those outputs have been provided. - This is a computer activity and a paper discussion in one. The main text of [Goldberg et al. (2010)](https://doi.org/10.1126/science.1194513) is very short, so we will just have quick look at it while some analyses run and discuss. # 1. The dataset and the question Most flowering plants are hermaphrodites, yet many can't fertilise themselves. In Solanaceae the style makes an S-RNase that destroys the RNA of growing pollen tubes, and the pollen carries proteins that disable every S-RNase except the one matching its own S allele. So a pollen grain is rejected by any style carrying its allele, including its own parent's. Three things about this system: 1. SI is the ancestral state ([Ioerger et al. 1990](https://doi.org/10.1073/pnas.87.24.9732)). 2. SI is easy to lose through loss-of-function mutations that kill S-RNase activity in the pistil ([Broz et al. 2021](https://doi.org/10.3389/fgene.2021.780793)). 3. Reversals from SC to SI are unlikely. One the function is lost, there is no selective pressure against more mutations ([Igić et al. 2006](https://doi.org/10.1073/pnas.0506283103)). About 41% of Solanaceae retain their ancestral SI state. [Goldberg et al. (2010)](https://doi.org/10.1126/science.1194513) hypothesized that while being self-compatible could be advantageous in the right context, SI lineages diversify fast enough to replace the lineages that lose it. [Stebbins (1957)](https://doi.org/10.1086/281999) suggested that selfing was an evolutionary dead end. Here, we have the dated tree and the SI or SC state of each species. We will reproduce their result with BiSSE, and explore some of the limitations. Two caveats about this data: - The tree has 356 of about 2,300 accepted species. We correct for that with sampling fractions. - BiSSE needs every species in one of two states. We consider any species with some SI as SI (some species have both), and everything else counts as SC. The SC group includes 25 dioecious species, which have no SI but cannot self either. So, consider "SC" here as "not SI". # 2. Reading the data into R ```{r} #| label: setup #| eval: false library(ape) # reading and handling trees library(diversitree) # BiSSE library(ggplot2) # plots library(patchwork) # putting several plots side by side theme_set(theme_bw(base_size = 11)) # one colour per breeding system, used in every figure sys_cols <- c("SI" = "#0072B2", "SI+SC" = "#56B4E9", "SC" = "#D55E00") two_cols <- sys_cols[c("SI", "SC")] # the two states BiSSE uses # the biggest genera, which we name in the figures big_genera <- c("Solanum", "Nicotiana", "Lycium", "Capsicum", "Physalis", "Nolana") ``` The tree and the traits are two plain files. The tree is the dated tree from the paper, already ultrametric (every tip is the same distance from the root), which `diversitree` insists on. The trait file has one row per species, with 0 for SI and 1 for SC. ```{r} #| label: read-data #| eval: false tree <- read.tree("data/Solanaceae_timetree.tre") traits <- read.table("data/trait_data.txt", header = TRUE) head(traits) tree # 356 tips is.ultrametric(tree) # TRUE, so diversitree will accept it # BiSSE wants a 0 or a 1 for every species, labelled with its name and in the same # order as the tips of the tree states <- traits$state names(states) <- traits$species states <- states[tree$tip.label] table(states) ``` You should have 135 SI and 221 SC species. We also have a table with the three-way summary: 0 is SI, 1 is SI+SC (the species is polymorphic) and 2 is SC. The two-state `state01` column is the same trait as above. ```{r} #| label: read-three-state #| eval: false three <- read.table("data/trait_data_3state.txt", header = TRUE) # the label for each code, in order from 0 to 2. R counts from 1, so we add 1 to each # code: code 0 picks the first name, code 1 the second, and so on. system_names <- c("SI", "SI+SC", "SC") three$system <- system_names[three$state012 + 1] table(three$system) # the two-state trait is the same thing with SI+SC folded into SI all(three$state01 == as.integer(three$state012 == 2)) all(three$state01 == traits$state) # a lookup for later: give it species names and it returns their breeding systems system_of <- three$system names(system_of) <- three$species ``` Now let's see how the two states are spread around the tree. ```{r} #| label: fan-tree #| eval: false #| fig-width: 8 #| fig-height: 8 # ladderize() sorts clades by size, which only changes how the tree is drawn tree_to_draw <- ladderize(tree) tips <- tree_to_draw$tip.label # trait.plot() wants the states in a data frame, with species as the row names ring <- data.frame(system = states[tips], row.names = tips) # name only the big genera; NA means no label genus <- sub("_.*", "", tips) # "Solanum_tuberosum" becomes "Solanum" genus_label <- ifelse(genus %in% big_genera, genus, NA) # type = "f" makes it a fan, and w sets the width of the ring trait.plot(tree_to_draw, ring, cols = list(system = unname(two_cols)), class = genus_label, type = "f", w = 1/30, margin = 0.3, lab = "", legend = FALSE, cex.lab = 0.9) legend("topleft", legend = names(two_cols), fill = two_cols, border = NA, bty = "n", cex = 0.9) ``` ![The 356 species, with each species' state in the ring around the tips. The root is 36 million years old.](figures/fan-tree.png){fig-alt="A circular tree of 356 species with a coloured ring around the tips. Blue self-incompatible and vermillion self-compatible tips are interspersed in most of the tree. Solanum has some long vermillion blocks, Nicotiana is nearly all vermillion, and Nolana is mostly blue."} SI and SC are interspersed throughout, which is expected under a multiple loss scenario. **Q1: Does it look like some clades have more SI than others?** # 3. Fitting BiSSE BiSSE ([Maddison et al. 2007](https://doi.org/10.1080/10635150701607033)) gives each state its own speciation rate ($\lambda$) and extinction rate ($\mu$), plus rates of change between the states. With SI as state 0 and SC as state 1: | `diversitree` | paper | meaning | |:---|:---|:------| | `lambda0`, `lambda1` | $\lambda_I$, $\lambda_C$ | speciation in SI and SC lineages | | `mu0`, `mu1` | $\mu_I$, $\mu_C$ | extinction in SI and SC lineages | | `q01` | $q_{IC}$ | loss of SI | | `q10` | $q_{CI}$ | regain of SI | `make.bisse()` returns a likelihood function, and `constrain()` makes a simpler model from it. The sampling fractions come from the paper's supplement: 16.2% of SI species and 15.0% of SC species are in the tree. ```{r} #| label: bisse-setup #| eval: false f <- c(0.162, 0.150) # SI first, then SC # lik is a function: give it the six rates and it returns the likelihood of our data lik <- make.bisse(tree, states, sampling.f = f) argnames(lik) # the six rates it asks for # rough values for the search to start from start <- starting.point.bisse(tree) ``` Is one set of rates enough for both states? If not, is it speciation that differs, extinction, or both? And should SI be allowed to come back? We can write each idea as a model and fit them one at a time, starting with the simplest. `find.mle()` climbs the likelihood from our starting values until it can't improve. First, a small helper to print what we want to see from each fit. ```{r} #| label: show-fit #| eval: false # print the fitted rates, the log likelihood and the AIC (lower AIC is better) show_fit <- function(fit) { print(round(fit$par, 3)) c(lnL = round(fit$lnLik, 2), AIC = round(AIC(fit), 1)) } ``` **Model 1: same rates.** SI and SC speciate and go extinct at the same rates, and SI is lost but never regained. `constrain()` ties one rate to another: `lambda1 ~ lambda0` sets SC's speciation rate equal to SI's, so there is one rate less to estimate, and `q10 ~ 0` stops SI from coming back. ```{r} #| label: fit-same #| eval: false same_model <- constrain(lik, lambda1 ~ lambda0, mu1 ~ mu0, q10 ~ 0) argnames(same_model) # only three rates are left to estimate # start[argnames(same_model)] picks out the starting values for just those rates same_rates <- find.mle(same_model, start[argnames(same_model)]) show_fit(same_rates) ``` A single speciation rate and a single extinction rate, both near 4/Ma and almost equal. Relaxing the contraints might improve the AIC. **Model 2: speciation differs.** SC gets its own speciation rate, but extinction is shared. ```{r} #| label: fit-speciation #| eval: false speciation_model <- constrain(lik, mu1 ~ mu0, q10 ~ 0) speciation_differs <- find.mle(speciation_model, start[argnames(speciation_model)]) show_fit(speciation_differs) ``` **Model 3: extinction differs.** The reverse: SC gets its own extinction rate and speciation is shared. ```{r} #| label: fit-extinction #| eval: false extinction_model <- constrain(lik, lambda1 ~ lambda0, q10 ~ 0) extinction_differs <- find.mle(extinction_model, start[argnames(extinction_model)]) show_fit(extinction_differs) ``` **Model 4: both differ.** This is the paper's model: both rates can differ, and SI never returns. Compare its rates with the first three. ```{r} #| label: fit-both #| eval: false both_model <- constrain(lik, q10 ~ 0) both_differ <- find.mle(both_model, start[argnames(both_model)]) show_fit(both_differ) ``` **Model 5: SI can come back.** The full model with nothing constrained. We do not think this is plausible, but will include it to hopefully make a point about hypothesis testing with models. ```{r} #| label: fit-return #| eval: false si_can_return <- find.mle(lik, start[argnames(lik)]) show_fit(si_can_return) ``` Now let's make our own table to compare all five. The delta_AIC is the different between each model AIC and the best AIC. Can you fill out the following table from your results? | Model | Parameters | lnL | AIC | delta AIC | |-------|------------|-----|-----|------------| | Same rates | | | | | | Speciation differs | | | | | | Extinction differs | | | | | | Both differ | | | | | | SI can return | | | | | Letting rates differ between SI and SC is the best model, but we do not consider this to be plausible. We will ignore that issue for now and focus on the rates estimated by the "both_differ" model. You can access the parameter estimates for this model by typing ```{r} both_differ$par ``` Can you use these estimates to calcualte the net diversification rates for SI ($r_0$) and SC lineages ($r_1$)? They should be pretty close to the paper's reported values. **Q2: Which diversification rate is higher, SI (state 0) or SC (state 1)? Did you get a negative number for one of them; what does that mean?** # 4. Bayesian analyses and testing the dead-end SC hypothesis. The paper and most studies use a Bayesian analysis, sampling many parameter combinations to assess uncertainty in the estiamtes. Remember, we suggested there were some non-identifiability issues between $\lambda$ and $\mu$. An MCMC wanders through the rates, visiting each combination in proportion to how well it fits the data and how plausible the prior says it is. Bayesian models need prior distribtions on all parameters. The prior here is the paper's: an exponential with a mean of 1/0.3 on every rate. We will try to talk about priors while your analysis runs Each iteraction of the MCMC calculates the likelihood again, so the computation can take a while. We will run a short one here and compare it with four long chains that I ran ahead of time. Run this chunk and let it go. It takes about 12 minutes on a laptop, and prints a line every 500 steps, so carry on with the reading below while you wait. ```{r} #| label: mcmc-short #| eval: false # the paper's model again: SI is lost but never regained one_way <- constrain(lik, q10 ~ 0) prior <- make.prior.exponential(0.3) # Start the chain at the maximum likelihood rates. w is how wide a step the sampler # tries first, and set.seed() makes your chain repeatable. set.seed(1) short <- mcmc(one_way, both_differ$par, nsteps = 5000, w = 1, prior = prior, print.every = 500) write.csv(short, "short_chain.csv", row.names = FALSE) # save it in case R crashes ``` Each row is one step: the step number `i`, the five rates and `p`, the log of likelihood times prior. Now the **trace**, which is each rate plotted against step. A chain that has found the right region looks like a fuzzy band that stays level. ```{r} #| label: mcmc-trace #| eval: false #| fig-width: 10 #| fig-height: 5 rate_names <- c("lambda0", "lambda1", "mu0", "mu1", "q01") par(mfrow = c(2, 3)) # two rows and three columns of plots for (rate in rate_names) { plot(short$i, short[[rate]], type = "l", xlab = "step", ylab = rate) } plot(short$i, short$p, type = "l", xlab = "step", ylab = "log posterior") par(mfrow = c(1, 1)) # back to one plot at a time ``` The log posterior levels off almost at once, but the rates wander slowly up and down, and λ and μ move together: a pair of larger values fits about as well as a pair of smaller ones, which we will meet again in section 6. The chain starts at the maximum likelihood estimate, which is not a typical draw from the posterior, so we drop the first 500 steps as **burn-in**. We then work out the same quantities as the table in section 4 and compare with the long chains. Also we can calculate the proportion of SI lineages expected by the model ($p_{SI}$) from the rates using the formula $p_{SI} = 1 - q_{IC} / (r_I - r_C)$, where $r_I = \lambda_I - \mu_I$ and $r_C = \lambda_C - \mu_C$. ```{r} #| label: posterior #| eval: false # The chains hold only the five rates in the model, so add the quantities built from # them. Each line does the sum for every row at once. add_derived <- function(x) { x$r_I <- x$lambda_I - x$mu_I x$r_C <- x$lambda_C - x$mu_C x$p_SI <- 1 - x$q_IC / (x$r_I - x$r_C) x } # our short chain, with the burn-in dropped and the columns named as in the paper mine <- short[short$i > 500, ] names(mine)[names(mine) %in% rate_names] <- c("lambda_I", "lambda_C", "mu_I", "mu_C", "q_IC") mine <- add_derived(mine) # the long chains: 4 chains, with 500 samples kept from each post <- read.csv("data/bisse_mcmc.csv") # one row per sample from the chains post <- add_derived(post) table(post$chain) # the same quantities as the table above, in the same order my_rates <- mine[names(paper)] post_rates <- post[names(paper)] post_table <- data.frame(short_mean = sapply(my_rates, mean), long_mean = sapply(post_rates, mean), Table_S1 = paper, short_sd = sapply(my_rates, sd), long_sd = sapply(post_rates, sd)) round(post_table, 3) ``` Now let's put the two side by side. The grey density is the longer chains, and the colour is your short one. ```{r} #| label: mcmc-compare #| eval: false #| fig-width: 10 #| fig-height: 3.5 # one density for the long chains and one for our chain, for a single quantity compare_chains <- function(column, x_label) { both <- data.frame(value = c(post[[column]], mine[[column]]), chain = rep(c("long chains", "our short chain"), times = c(nrow(post), nrow(mine)))) ggplot(both, aes(value, fill = chain)) + geom_density(alpha = 0.5, colour = NA) + scale_fill_manual(values = c("long chains" = "grey50", "our short chain" = "#D55E00")) + labs(x = x_label, y = "posterior density", fill = NULL) } compare_chains("lambda_I", "SI speciation") + compare_chains("mu_C", "SC extinction") + compare_chains("r_C", "SC net diversification") + compare_chains("q_IC", "loss of SI") + plot_layout(nrow = 1, guides = "collect") ``` Why do the speciation and extinction traces wander so slowly? The sampler moves one rate at a time, and some rates are tied tightly together. Let's plot pairs of parameters against each other, with the long chains in grey and our short chain in colour. ```{r} #| label: mcmc-scatter #| eval: false #| fig-width: 10 #| fig-height: 2.8 # one point per sample for a pair of parameters scatter_pair <- function(x_column, y_column, x_label, y_label) { ggplot() + geom_point(data = post, aes(.data[[x_column]], .data[[y_column]]), size = 0.4, alpha = 0.3, colour = "grey60") + geom_point(data = mine, aes(.data[[x_column]], .data[[y_column]]), size = 0.4, alpha = 0.4, colour = "#D55E00") + labs(x = x_label, y = y_label) } scatter_pair("lambda_I", "mu_I", "SI speciation", "SI extinction") + scatter_pair("lambda_C", "mu_C", "SC speciation", "SC extinction") + scatter_pair("r_I", "q_IC", "SI net diversification", "loss of SI") + scatter_pair("r_I", "r_C", "SI net diversification", "SC net diversification") + plot_layout(nrow = 1) ``` With the seed above, our 5,000 steps land close to the long chains: the means for the net diversification rates and $p_{SI}$ agree to within about 0.02, and the speciation and extinction means are off by up to 0.15. The scatterplots show why those two are less reliable. Speciation and extinction are highly correlated, and since the chain changes one rate at a time, that parameter space is explored slowly. Note the last panel is a little different, more of a cloud. This is positive as it signals there is some independence between the net diversifiction rates. If your chain crashed or you ran out of time, `short <- read.csv("data/short_chain_example.csv")` reads in the one used for this text, and you can carry on from the trace plot. Means and standard deviations from the long chains match Table S1 closely, so **we have reproduced the paper**. From here on we use the long chains, in `post`. To see how the next chunks do with your own chain, replace `post` by `mine` instead. Now, we can get to the core hypothesis test. Do SI lineages diversify faster than SC lineages? ```{r} #| label: posterior-plot #| eval: false #| fig-width: 10 #| fig-height: 3.8 # Draw the posterior of one rate for SI and SC on the same axes. ggplot wants all of # the values in one column, with a second column saying which state each one is. two_densities <- function(si_column, sc_column, x_label) { long <- data.frame(rate = c(post[[si_column]], post[[sc_column]]), state = rep(c("SI", "SC"), each = nrow(post))) long$state <- factor(long$state, levels = c("SI", "SC")) # SI first in the legend ggplot(long, aes(rate, fill = state)) + geom_density(alpha = 0.6, colour = NA) + scale_fill_manual(values = two_cols) + labs(x = x_label, y = "posterior density", fill = NULL) } speciation_plot <- two_densities("lambda_I", "lambda_C", "speciation rate (per Myr)") extinction_plot <- two_densities("mu_I", "mu_C", "extinction rate (per Myr)") net_div_plot <- two_densities("r_I", "r_C", "net diversification (per Myr)") + geom_vline(xintercept = 0, linetype = 2) # one point per sample, to show how SC's two rates rise and fall together sc_rates_plot <- ggplot(post, aes(lambda_C, mu_C)) + geom_point(size = 0.3, alpha = 0.3, colour = two_cols[["SC"]]) + geom_abline(linetype = 2, colour = "grey50") + labs(x = "SC speciation", y = "SC extinction") # patchwork lets us add plots together speciation_plot + extinction_plot + net_div_plot + sc_rates_plot + plot_layout(nrow = 1, guides = "collect") ``` ![Posterior distributions, as in the paper's Figure 1. Speciation and extinction are uncertain, and net diversification is much tighter. The right panel shows why: λ~C~ and μ~C~ rise and fall together. Points above the dashed line have λ < μ.](figures/posterior-plot.png){fig-alt="Four panels. In the first three, blue is SI and vermillion SC: SC speciation is higher than SI, SC extinction is higher than SI, and SI net diversification is positive while SC is negative, with no overlap. The fourth panel plots SC extinction against SC speciation as a tight diagonal band that sits just above the one-to-one line."} So, SI species grow at $r_{I}$ and evolve into SC at $q_{IC}$. So SI can only hold on if it out-diversifies SC by more than the rate it is lost. When that is true, the propotion of SI species can then be expressed as: $$p_{SI} = 1 - \frac{q_{IC}}{r_I - r_C}$$ The posterior mean for $p_{SI}$ from the `post` object is 35% of species. The observed proportion is 41% across Solanaceae. Pretty close. **Q3: Do these model estimates explain the persistence of SI lineages despite its susceptibility to loss-of-function?** Stebbins suggested there are many SC species today only because SI lineages keep producing new ones. # 5. A trait that does nothing Five years after this paper, the first author published a study suggesting caution in interpreting BiSSE results ([Rabosky & Goldberg 2015](https://doi.org/10.1093/sysbio/syu131)). They simulated traits with **no effect** on diversification on real trees, and BiSSE often reported an effect anyway. Real trees are uneven for many reasons, so any trait that happens to sit in a fast clade can get the credit. Could that be happening here? We can run the same test on this tree, using the model where SI can only be lost. I simulated 300 traits that can only be lost, at a rate that leaves about as many tips in the starting state as there are SI species today. These traits have no connection to diversification. For each one I fit BiSSE as in section 3 and we can retain the difference in net diversification between the two states ($r_{I} - r_{C}$) as our summary statistic. We can check how different our observed trait is from the simulated distribution to find out if the trait-dependent rates are surprising. ```{r} #| label: neutral-plot #| eval: false #| fig-width: 8 #| fig-height: 3.8 neutral <- read.csv("data/neutral_oneway_fits.csv") # the 300 simulated traits real <- read.csv("data/neutral_oneway_observed.csv") # the real trait neutral$r_gap <- neutral$r_I - neutral$r_C real$r_gap <- real$r_I - real$r_C ggplot(neutral, aes(r_gap)) + geom_histogram(binwidth = 0.05, boundary = 0, fill = "grey60") + geom_vline(xintercept = real$r_gap, colour = two_cols[["SC"]], linewidth = 1) + labs(x = "net diversification of SI minus SC (per Ma)", y = "simulated traits") ``` ![The gap in net diversification between the two states for 300 simulated traits with no effect on diversification, which can only lose the first state. The vermillion line is the real trait.](figures/neutral-plot.png){fig-alt="A histogram of the gap in net diversification for 300 simulated traits, all between about −0.1 and 0.35 and centred near 0.1. The vermillion line for the real trait, at 0.74, lies far to the right of all of them."} The real trait associated with diversification is very different from a trait that is not. **Q4: There is a real difference in diversification between states. Does that show that SI itself changes diversification? Can you think of other factors that would lead to transitions from SI to SC AND differences in diversification rates?** # 6. Hidden states BiSSE imposes a limited hypothesis test. - $H_O$: rates do not vary over the tree - $H_A$: rates do vary over the tree, and that depends on the trait [Beaulieu & O'Meara (2016)](https://doi.org/10.1093/sysbio/syw022) aimed to correct this limitation by letting rates vary for another reason, some unknown trait. In the CID models rates depend only on a hidden trait with 2 or 4 classes, and in HiSSE they depend on the observed (e.g. SI vs SC) *and* a hidden trait. The bigger models take about half an hour each to fit, so I ran them ahead of time with the paper's assumptions. One fit looks like this, shown for reading and not for running. ```r library(hisse) dat <- data.frame(species = tree$tip.label, state = states) # the allowed changes: between SI and SC, and between two hidden classes transitions <- TransMatMakerHiSSE(hidden.traits = 1, make.null = TRUE) loss_only <- ParDrop(transitions, 1) # rate 1 is SC -> SI, which we remove # Rates are given in the order 0A, 1A, 0B, 1B, and equal numbers share a rate. So this # is CID-2: rates differ between hidden classes A and B, but not between SI and SC. cid2_rates <- c(1, 1, 2, 2) hisse(tree, dat, f = c(0.162, 0.150), hidden.states = TRUE, trans.rate = loss_only, turnover = cid2_rates, eps = cid2_rates) ``` ```{r} #| label: hisse-models #| eval: false hisse_models <- read.csv("data/hisse_models.csv") # leave out the versions that share one extinction fraction across states keep <- !grepl("_eps", hisse_models$model) hisse_models <- hisse_models[keep, ] # k is the number of parameters, and dAICc is how far each model is behind the best hisse_table <- hisse_models[, c("model", "k", "lnL", "dAICc")] hisse_table$lnL <- round(hisse_table$lnL, 1) hisse_table$dAICc <- round(hisse_table$dAICc, 1) print(hisse_table, row.names = FALSE) ``` BiSSE is 246 AICc units behind the best model, and it is even worse than CID-2, where we have two unknown traits (but transition between both states are possible unlike our constrained BiSSE model). The fair comparison is the top two rows. HiSSE and CID-4 have the same number of parameters, and HiSSE wins by 12 AICc units. Here are the rates in the best model. ```{r} #| label: hisse-rates #| eval: false hisse_rates <- read.csv("data/hisse_rates.csv") # keep just the rates from the HiSSE model best <- hisse_rates[hisse_rates$model == "HiSSE", ] rate_columns <- c("lambda", "mu", "r") best[rate_columns] <- round(best[rate_columns], 2) print(best, row.names = FALSE) ``` Compare SI with SC **within each hidden class**. SI has the higher net diversification in both, 0.73 against −0.14 in class A and 0.44 against −0.01 in class B. The core biological result is that SI lineages have higher net diversification compared to SC survives the hidden-class correction. Can we find the lineages that belong to the strange fast turnover class? ```{r} #| label: hisse-tips-plot #| eval: false #| fig-width: 10 #| fig-height: 4 tip_rates <- read.csv("data/hisse_tip_rates.csv") # one row per species # state 0 is SI and state 1 is SC, as before tip_rates$system <- ifelse(tip_rates$state == 0, "SI", "SC") tip_rates$system <- factor(tip_rates$system, levels = c("SI", "SC")) # Lump the small genera together. rev() puts the first genus at the top of the plot. tip_rates$group <- ifelse(tip_rates$genus %in% big_genera, tip_rates$genus, "others") tip_rates$group <- factor(tip_rates$group, levels = rev(c(big_genera, "others"))) # one row of dots per genus, with the x axis added separately for each panel below dot_plot <- ggplot(tip_rates, aes(y = group, colour = system)) + geom_jitter(height = 0.2, width = 0, size = 1, alpha = 0.6) + scale_colour_manual(values = two_cols) + labs(y = NULL, colour = NULL) net_div_panel <- dot_plot + aes(x = net.div) + geom_vline(xintercept = 0, colour = "grey60") + labs(x = "net diversification (per Myr)") net_div_panel ```  The two hidden classes roughly mean "*Solanum*" and "everything else": median turnover is 36 per Myr in *Solanum* and about 2.7 in most other groups. SI and SC separate within every genus, but part of that is by construction, since each tip's rate depends on its state. **Q5: The fast hidden class is mostly *Solanum*. What do you think the hidden trait could be? What could you do to test that selection acts on SI itself and not an associated trait?** # 7. Takeaways 1. **SC is a capacity, not a mating system.** The tree records the loss of SI, and selfing is a separate component. 2. **Trees can be used to estimate diversification rates.** There was high uncertainty in the speciation and extinction rates, but net diversification was quite reliable, even with a short MCMC chain. 3. **A trait with no effect can't produce the observed rates under BiSSE.** When SI can only be lost, none of 300 simulated random traits were similar to our observed rates. 4. **The effect survives hidden states, but there is more to learn.** SI out-diversifies SC in both hidden classes, but something else we do not know is important to *Solanum*.

`diversitree`	paper	meaning
`lambda0`, `lambda1`	\(\lambda_I\), \(\lambda_C\)	speciation in SI and SC lineages
`mu0`, `mu1`	\(\mu_I\), \(\mu_C\)	extinction in SI and SC lineages
`q01`	\(q_{IC}\)	loss of SI
`q10`	\(q_{CI}\)	regain of SI