Uyeda et al. 2017

Uyeda et al. did something similar in their 2017 paper The evolution of energetic scaling across the vertebrate tree of life. They stitched together fish, amphibian, squamate, bird, and mammal phylogenies together to visualize metabolic rate across all vertebrates:

At first, I was very hopeful that I would be able to download this supertree and prune it to the taxa in my analysis since the authors were using the same clade specific phylogenies that I was. The phylogeny is available on Data Dryad (https://datadryad.org/resource/doi:10.5061/dryad.3c6d2). Unfortunately, after downloading that phylogeny and pruning it to include species I used in my analysis, I ended up with approximately 15% of the species I analyzed in the resulting tree.

The original phylogenies

I was working with four separate phylogenies: amphibians, squamates, birds, and mammals. For amphibians, I used a congruified time-tree from the PhyloOrchard package (O’Meara et al. 2013) that was constructed using the Alfaro et al. timetree of gnathostomes (Alfaro et al. 2009) as the reference and the Pyron and Wiens amphibian phylogeny as the target (Pyron and Wiens 2011).

length(amphibiantree$tip.label)
## [1] 2871
plot(amphibiantree, type = "fan", show.tip.label = FALSE)

For squamates, I used the Zheng and Wiens time-calibrated phylogeny (Zheng and Wiens 2016).

length(squamatetree$tip.label)
## [1] 378
plot(squamatetree, type = "fan", show.tip.label = FALSE)

For birds, I used the Jetz phylogeny (Jetz et al. 2012).

length(birdtree$tip.label)
## [1] 9993
plot(birdtree, type = "fan", show.tip.label = FALSE)

For mammals, I used the supertree from Fritz et al. 2009.

length(mammaltree$tip.label)
## [1] 5020
plot(mammaltree, type = "fan", show.tip.label = FALSE)

Reading in trees

The first step of the process is reading in the individual phylogenies you want to stitch together. This step is straightforward, with one exception: you cannot have species that are present in multiple of the individual trees. For example, my squamate phylogeny included Gallus gallus (red junglefowl) and Dromaius novaehollandia (emu). Since these species were also present in my bird phylogeny, I got the following error: Found matching tips in 'subtree' and 'phy'. To solve this problem, I just removed these tips from the squamate tree:

squamatetree <- drop.tip(phy = squamatetree, tip = c("Gallus_gallus", "Dromaius_novaehollandiae"))

I ended up with the following list of trees and corresponding tip labels:

tree_list <- list(amphib=amphibiantree, birds=birdtree, squam=squamatetree, mamm=mammaltree)
class(tree_list) <- "multiPhylo"

Make a tree with orders

In Uyeda et al. (2017), the authors were creating a phylogeny for all vertebrates, but for my analyses I was only examining tetrapods, so I didn’t have a fish phylogeny to include. The original code from Uyeda et al. to create a tree with the 5 vertebrate orders looks like this:

tip.labels <- c("fish", "amphib", "squam", "birds", "mamm")

## Make a tree with just orders:
edge <- matrix(c(9, 4,
  9, 3,
  8, 5,
  8, 9,
  7, 8,
  7, 2,
  6, 7,
  6, 1), byrow=TRUE, ncol=2)
## Dates from Timetree of life (timetree.org)
edge.length <- c(274.9, 274.9, 324.5, 324.5-274.9, 382.9-324.5, 382.9, 454.6-382.9 , 454.6)
Nnode <- 4
ordertree <- list(edge=edge, Nnode=Nnode, tip.label=tip.labels, edge.length=edge.length)
class(ordertree) <- 'phylo'
plot(ordertree)

To visualize the results, we can add tip labels, node labels, and edge labels to the tree with the branch lengths:

plot(ordertree)
tiplabels()
nodelabels()
edgelabels(ordertree$edge.length, bg="black", col="white", font=2)

# remove "fish" from tip.labels:
tip.labels <- c("amphib", "squam", "birds", "mamm")

edge <- matrix(c(7, 3,
  7, 2,
  6, 4,
  6, 7,
  5, 6,
  5, 1), byrow=TRUE, ncol=2)

edge.length <- c(274.9, 274.9, 324.5, 324.5-274.9, 382.9-324.5, 382.9)

Nnode <- 3

ordertree <- list(edge=edge, Nnode=Nnode, tip.label=tip.labels, edge.length=edge.length)
class(ordertree) <- 'phylo'
plot(ordertree)
edgelabels(ordertree$edge.length, bg="black", col="white", font=2)

Node dates

…Not so fast. I ended up with an additional problem I needed to solve before grafting the trees together. I ran into the error 'split_age' is inconsistent with edge lengths in 'phy', which means that the earliest node in one of my individual phylogenies was older than the node age I gave in edge.length. By using debug, I was able to tell that the error occurred when I added the squamate tree. The oldest node in my squamate tree was 277.8 million years ago, but I had set the divergence time between birds and squamates at 274.9 mya, so R was having problems. The species causing the problem was the tuatara, which is the only surviving member of its order.

I had a couple of choices: either delete the tuatara from the squamate phylogeny or increase the age of the last common ancestor of birds and squamates when I created the vector edge.length. I chose to do the latter because why get rid of such a cool animal!

I went to http://timetree.org/ to see if I could find a reasonable range of estimates for this node. According to the website, which allows you to search for the divergence time between any two taxa, the estimated divergence of birds and squamates occurred 280 mya.

So I ended up with the following code and order tree:

tip.labels <- c("amphib", "squam", "birds", "mamm")
edge <- matrix(c(7, 3,
  7, 2,
  6, 4,
  6, 7,
  5, 6,
  5, 1), byrow=TRUE, ncol=2)
edge.length <- c(280, 280, 324.5, 324.5-274.9, 382.9-324.5, 382.9)
Nnode <- 3
ordertree <- list(edge=edge, Nnode=Nnode, tip.label=tip.labels, edge.length=edge.length)
class(ordertree) <- 'phylo'
plot(ordertree)
edgelabels(ordertree$edge.length, bg="black", col="white", font=2)

tree_list <- list(amphib=amphibiantree, birds=birdtree, squam=squamatetree, mamm=mammaltree)
class(tree_list) <- "multiPhylo"

Grafting the trees

The final step is grafting the individual trees onto the order tree in the proper place.

#Add taxonomic information to tree
otax <- data.frame("Class"= ordertree$tip.label, "Superclass"=c(rep("Tetrapoda",2)))
rownames(otax) <- ordertree$tip.label
classtree <- nodelabel.phylo(ordertree, otax, ncores=1)

res <- glomogram.phylo(classtree, tree_list)
plot(res, type = "fan", show.tip.label = FALSE)

Voila - a tree with 18262 species of tetrapods!