Phylogenetic community assembly with the Fish Tree of Life

Jonathan Chang

Here’s a quick example to show how we could use the fishtree package to conduct some phylogenetic community analyses. First, we load fishtree and ensure that the other packages that we need are installed.

library(ape)
library(fishtree)
requireNamespace("rfishbase")
#> Loading required namespace: rfishbase
requireNamespace("picante")
#> Loading required namespace: picante
requireNamespace("geiger")
#> Loading required namespace: geiger
#> Registered S3 method overwritten by 'geiger':
#>   method            from
#>   unique.multiPhylo ape

Next we’ll start downloading some data from rfishbase. We’ll be seeing if reef-associated ray-finned fish species are clustered or overdispersed in the Atlantic, Pacific, and Indian Oceans.

# Get reef-associated species from the `species` table
reef_species <- rfishbase::species(fields = c("Species", "DemersPelag"))
reef_species <- reef_species[reef_species$DemersPelag == "reef-associated", ]

# Get native and endemic species from the Atlantic, Pacific, and Indian Oceans
eco <- rfishbase::ecosystem(species_list = reef_species$Species)
valid_idx <- eco$Status %in% c("native", "endemic") & eco$EcosystemName %in% c("Atlantic Ocean", "Pacific Ocean", "Indian Ocean") 
eco <- eco[valid_idx, c("Species", "EcosystemName")]

# Retrieve the phylogeny of only native reef species across all three oceans.
phy <- fishtree_phylogeny(species = eco$Species)
#> Warning: Requested 5788 but only found 3997 species. Missing names:
#> * Conger esculentus
#> * Lutjanus guttatus
#> * Lutjanus inermis
#> * Lutjanus jordani
#> * Lutjanus lemniscatus
#> * ...and 1786 others

We’ll have to clean up the data in a few ways before sending it to picante for analysis. First, we’ll need to convert our species-by-site data frame into a presence-absence matrix. We’ll use base::table for this, and use unclass to convert the table into a standard matrix object.

sample_matrix <- unclass(table(eco))
dimnames(sample_matrix)$Species <- gsub(" ", "_", dimnames(sample_matrix)$Species, fixed = TRUE)

Next, we’ll use geiger::name.check to ensure the tip labels of the phylogeny and the rows of the data matrix match each other.

nc <- geiger::name.check(phy, sample_matrix)
sample_matrix <- sample_matrix[!rownames(sample_matrix) %in% nc$data_not_tree, ]

Finally, we’ll generate the cophenetic matrix based on the phylogeny, and transpose the presence-absence matrix since picante likes its columns to be species and its rows to be sites.

cophen <- cophenetic(phy)
sample_matrix <- t(sample_matrix)

We’ll run ses.mpd and ses.mntd with only 100 iterations, to speed up the analysis. For a real analysis you would likely increase this to 1000, and possibly test other null models if your datasets have e.g., abundance information.

picante::ses.mpd(sample_matrix, cophen, null.model = "taxa.labels", runs = 99)
#>                ntaxa  mpd.obs mpd.rand.mean mpd.rand.sd mpd.obs.rank
#> Atlantic Ocean   637 238.3167      232.2084   1.8066282          100
#> Indian Ocean    1222 233.5305      231.8798   1.2929589           86
#> Pacific Ocean   1511 232.1380      231.8682   0.7284014           65
#>                mpd.obs.z mpd.obs.p runs
#> Atlantic Ocean 3.3810650      1.00   99
#> Indian Ocean   1.2766983      0.86   99
#> Pacific Ocean  0.3702948      0.65   99
picante::ses.mntd(sample_matrix, cophen, null.model = "taxa.labels", runs = 99)
#>                ntaxa mntd.obs mntd.rand.mean mntd.rand.sd mntd.obs.rank
#> Atlantic Ocean   637 41.62559       47.83711    1.4396703             1
#> Indian Ocean    1222 34.98367       37.89510    0.6658603             1
#> Pacific Ocean   1511 34.00154       34.90489    0.4914791             4
#>                mntd.obs.z mntd.obs.p runs
#> Atlantic Ocean  -4.314548       0.01   99
#> Indian Ocean    -4.372430       0.01   99
#> Pacific Ocean   -1.838023       0.04   99

The Atlantic and Indian Oceans are overdispersed using the MPD metric, and all three oceans are clustered under the MNTD metric. MNTD is thought to be more sensitive to patterns closer to the root of the tree. We can confirm these patterns visually:

plot(phy, show.tip.label = FALSE, no.margin = TRUE)
obj <- get("last_plot.phylo", .PlotPhyloEnv)

matr <- t(sample_matrix)[phy$tip.label, ]
xx <- obj$xx[1:obj$Ntip]
yy <- obj$yy[1:obj$Ntip]
cols <- c("#1b9e77", "#d95f02", "#7570b3")
for (ii in 1:ncol(matr)) {
  present_idx <- matr[, ii] == 1
  points(xx[present_idx] + ii, yy[present_idx], col = cols[ii], cex = 0.1)
}