In order to begin exploring evolutionary history in the context of biodiversity analyses, we will first explore taxonomic and then phylogenetic diversity analyses.
Introduction to Taxonomic Diversity Analyses
In the first section of this lab, we will examine the patterns of average body size across 198 families of animals (kingdom Animalia) from a the 2012 article “Can clade age alone explain the relationship between body size and diversity?.”(Etienne et al. 2012)
Load in “Etienne_etal_metfam_data.txt” as df animalia, noting that this is a different filetype than we have used before. Explore the file outside of Rstudio. Now explore the read.table() help page to try to figure out which parameter values to specify to load the file in properly:
animalia <-read.table("Lab_9_data/Etienne_etal_metfam_data.txt", header =TRUE, sep ="\t", as.is =TRUE)head(animalia)
Take a moment to get familiar with the data. Describe the dataset:
dim(animalia)
[1] 198 11
The dataframe consists of 198 family rows, each with various summary measures for a different animal family, including the average body size of species within that family ($bodysize), the crown age of the family ($age), and its approximate number of species ($richness).
Crown Age:
Age of the most recent common ancestor of all species in a clade
Stem Age:
Time since an ancestral lineage diverged from its sister clade
Plot family mean body size against species richness:
ggplot(animalia, aes(x = bodysize, y = richness)) +geom_point() +labs(title ="Animal family richness vs average body size",x ="Average Body Size",y ="Species Richness")
We will now plot the relationship between log-transformed body size and log-transformed richness:
ggplot(animalia, aes(x =log10(bodysize), y =log10(richness))) +geom_point() +labs(title ="Animal family log richness vs log average body size",x ="Log10(Average Body Size)",y ="Log10(Species Richness)")
Which plot is more informative (un-transformed or log-transformed) and why? Think beyond visual analysis but also biological significance.
Log-transformed
Un-transformed is squished to one side on the x-axis so cannot distinguish differences among really small-sized families. Log transformation spreads the values more evenly along the x-axis.
Log transformation also makes more biological sense: a 100g difference in body size has a profound ecological significance if we are comparing two animals that weigh less than 200g, but is trivial if we are comparing two animals that weigh more than 100kg.
Another way to think about this is that differences in log-space are proportional: so a 10% change to a 10kg organism and a 10% change to a 100kg organism is recorded as the same amount of change in log space
On the basis of this plot of average body size against species richness across all animal families, what conclusion would you draw about the association between these two variables?
There is a negative relationship between the average body size of a family and its richness, though with large variation in species richness among families.
Taxonomic bias is a common issue in biodiversity studies. Taxonomic bias occurs when the number of species in a group is not proportional to the number of species in the real world. For example, there are more species of arthropods than chordates, but the number of species in a dataset may not reflect this.
Create a stacked geom_bar() plot to check for taxonomic bias in our dataset to the class level by plotting by phylum and filling by class:
With two intuitive lines of code, we are able to generate a useful plot that helps us to see the bias in our dataset.
What bias do you observe?
Over-representation of chordates and specifically mammals (~5.4K species) versus arthropods (>1M species), mollusks (~100K species), echinodermata (10k), bryozoa (~5k)
Using the dplyr package to create a new column in our animalia df called is_chordate. Use the ifelse() function to label if from phylum (label as chordata) or not (label as nonchordata)
animalia <- animalia %>%mutate(is_chordate =ifelse(phylum =="Chordata", "chordata", "nonchordata"))head(animalia)
Now repeat your log-transformed plot of body size and richness, but color by is_chordate
ggplot(animalia, aes(x =log10(bodysize), y =log10(richness), color = is_chordate)) +geom_point() +labs(title ="Animal family log richness vs log average body size",x ="Log10(Average Body Size)",y ="Log10(Species Richness)",color ="Phylum")
Does the plot with chordates and non-chordates distinguished change your interpretation of the association between body size and species richness? If so, how?
The negative relationship is now less clear; the chordates and non-chordates seem to cluster into two separate groups on the plot. The non-chordates all tend to have small family-level average body sizes and high family richness, and vice vera for chordates.
Also, the negative relationship between mean body size and species richness is much less clear (and almost disappears), especially in the Chordata (red points).
Subset for just mammals and plot log-transformed bodysize by log transformed richness using ggplot2.
Tip: you can combine a dplyr and ggplot() pipe together to avoid creating another new dataframe.
animalia %>%filter(class =='Mammalia') %>%ggplot(aes(x =log10(bodysize), y =log10(richness))) +geom_point() +labs(title ="Log richness vs. bodysize of mammals")
Are there any noteworthy or consistent patterns of variation within this group?
The pattern is weaker: there is low richness and large average body sizes overall. The negative relationship between body size and species richness is less clear in Mammalia than in the earlier plot of all Animalia families.
To dive deeper into this question of the relationship between richness and bodysize, let’s look at individual species instead of family averages.
We will work with the panTHERIA dataset(Jones et al. 2009), a comprehensive data set of all known extant and recently extinct mammals, compiled from various literature sources.
Read in panTHERIA dataset provided to you "PanTHERIA_1-0_WR05_Aug2008.txt", assign to object pantheria, use the function dplyr::rename() to change the body mass column name to bodysize, and explore the dataset, especially regarding bodysize.
Some options to explore could include checking the range() of bodysize values, check for NA values with is.na(), check for non-numerical element is is.nan(), etc.
There are negative numbers in our dataset which does not make sense. Warning that these negative numbers disappeared with our log-transformed plot, so it is good that we did not rely on this plot alone.
It looks like the number “-999” was used a placeholder for entries that do not have a bodysize value instead of NA or a non-numeric value.
Let’s clean up our dataset by filtering out these values
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Does the body-size histogram for mammals offer evidence for a positive association between body size and species richness within this taxon? What features of the shape of the histogram might help us understand the body size/species richness association?
The histogram shows that there are many small mammal species. This is consistent with an association between high speciation rate and small body size.
The histogram is heavily skewed to the right (longer tail to the right), with a long upper tail (most species are small; there are a few really big species).
Howevever, interestingly, the very smallest mammals are not the most diverse. The first 2 size classes have lower species diversity than the next two size classes.
Does this trend hold for mammalian orders?
Plot a histogram of bodysize faceted by every order in our pantheria dataset. Set scales = "free" noting what happens to the y- and x-axes between plots when we change that parameter.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
What do you observe?
Many orders are right-skewed in log space, with some exceptions like Artiodactyla, Diprotodontia, and Primates.
Side question: what is dangerous about setting our parameter to scales = "free"?
It can be deceiving because it does not allow us to readily compare values across plots
One hypothesis for the high amount of species with small body size is that small body size promotes more rapid diversification. If this is true, then we might expect to see a negative association if we plot mean body size of animal families against their rate of diversification.
A simple way to summarize the average diversification rate over the lifetime of a clade is to divide the log of extant species richness (N) by its estimated age (t).
For the animal dataset, we can plot \(\frac{log(N)}{t}\) against mean body size.
Calculate divrate For each family within the animalia df using the mutate() function.
animalia <- animalia %>%mutate(divrate =log(richness) / age)head(animalia)
Let’s generate log-transformed plots body size vs. species richness and plot body size vs. diversification rate, coloring by Chordata vs. nonchordata. Assign these plots to p1 and p2.
(Note that we have already taken the log of richness to calculate divrate so we do not take the log of divrate)
p1 <-ggplot(animalia, aes(x =log10(bodysize), y =log10(richness), color = is_chordate)) +geom_point() +labs(title ="Animal family log richness vs log average body size",x ="Log10(Average Body Size)",y ="Log10(Species Richness)",color ="Phylum")p2 <-ggplot(animalia, aes(x =log10(bodysize), y = divrate, color = is_chordate)) +geom_point() +labs(title ="Animal diversification vs log average body size",x ="Log10(Average Body Size)",y ="Diversification Rate",color ="Phylum")
Patchwork
The package patchwork(Pedersen 2022) allows you to combine multiple plots into a single plot using intuitive syntax:
Using patchwork display our plots p1 and p2 side by side for comparison
p1 / p2
On the basis of these plots, can you conclude that families of smaller species have diversified more rapidly? If not, what other explanations could there be for the difference in mean species richness between chordate and non-chordate families?
No, we cannot conclude that families of smaller species have diversified more rapidly. There is no strong relationship in the 2nd plot between body size and diversification rate.
Other explanations?
Small-bodied families are more diverse, but they also tend to be older so they have had more time to accumulate more species.
Perhaps Chordate families need more “room” because they are bigger, but also need a greater variety of habitat types than smaller animals. Smaller animals specialize on finer subdivisions of habitat, so a given area can support more niches for smaller species. This concept is referred to as the “fractal nature of habitat surfaces.” The speciation (or diversification) rates could be the same, but the standing diversity is just higher for the non-chordates.
Now that we have visualized a pattern relevant to the body size-diversification rate hypothesis, we will apply some simple statistical tests. There are various tests we could apply. Based on the plots we might conclude, for example, that we are justified in grouping families into chordates and non-chordates, and testing for a difference between these two groups in mean diversification rate of families.
First, create a df of just chordates and another of just non-chordates
chordates <- animalia %>%filter(is_chordate =="chordata")nonchordates <- animalia %>%filter(is_chordate =="nonchordata")
Let’s plot mass vs. diversification rate and add a simple regression line using a linear model with geom_smooth().
ggplot(data = chordates, aes(x =log10(bodysize), y = divrate)) +geom_point() +geom_smooth(method ="lm", se =FALSE) +labs(title ="Chordates")
`geom_smooth()` using formula = 'y ~ x'
ggplot(data = nonchordates, aes(x =log10(bodysize), y = divrate)) +geom_point() +geom_smooth(method ="lm", se =FALSE) +labs(title ="Not chordates")
`geom_smooth()` using formula = 'y ~ x'
After visually investigating this trend in the plot, let’s check the output of the linear model (note that the order of independent/dependent variables is different in this function than the one above)
Call:
lm(formula = nonchordates$divrate ~ log10(nonchordates$bodysize))
Residuals:
Min 1Q Median 3Q Max
-0.051727 -0.018285 -0.004638 0.009594 0.136293
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.040810 0.010817 3.773 0.000338 ***
log10(nonchordates$bodysize) -0.005188 0.002500 -2.075 0.041699 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03491 on 69 degrees of freedom
Multiple R-squared: 0.05875, Adjusted R-squared: 0.0451
F-statistic: 4.306 on 1 and 69 DF, p-value: 0.0417
Residuals: The difference (distance) between the actual values and the regression line.
As a simple step, we can use these values to determine if there are outliers that we should consider eliminating if min or max values are much larger than the median.
We can also perform separate statistics on the residuals themselves as a best practice.
A lower “residual standard error (RSE)” overall indicates a better fit
Coefficients:
Intercept: The predicted value of the dependent variable when the independent variable is equal to zero
Slope: The estimated effect of the independent variable on the dependent variable
The standard error can be used to indicate the precision of these estimates
The “Pr(>|t|)” column gives us the p-value for the intercept and slope
R-squared:
A “multiple R-squared” value closer to 1 indicates a better fit.
The “adjusted R-squared” value can be negative because it is adjusted for the number of independent variables in the model.
With this information, does this test support the hypothesis that families of small mean body size have diversified faster than those of larger mean body size? What about for chordates vs. nonchordates?
Visually there is a small negative trend for body size vs. diversification rate for nonchordates.
The p-values and even spread of our outliers indicate that our results likely did not arise by chance.
However, the slope estimate is actually quite low (-0.005188) and the R2 is low, so there is a poor fit, so there is not a strong trend
However, the intercept appears to be a good predictor of the value of the dependent variable as the independent variable approaches zero.
Consider this: could it be that the relationship we see between richness vs. body size is because older families differ in average body size and species richness compared to younger families?
Explore this by plotting this relationship and generating a linear model summary.
ggplot(data = animalia, aes(x =log10(age), y =log10(richness))) +geom_point() +geom_smooth(method ="lm", se =FALSE)
Call:
lm(formula = log10(animalia$age) ~ log10(animalia$bodysize))
Residuals:
Min 1Q Median 3Q Max
-0.90436 -0.15704 -0.00856 0.15726 0.73932
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.689145 0.019547 86.41 <2e-16 ***
log10(animalia$bodysize) -0.075195 0.006706 -11.21 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2462 on 196 degrees of freedom
Multiple R-squared: 0.3908, Adjusted R-squared: 0.3877
F-statistic: 125.8 on 1 and 196 DF, p-value: < 2.2e-16
What do you discover? Do the patterns make biological sense given what you know about how diversification rate is measured (number of species / clade age) and about the evolution of life on earth? Why or why not?
Based on the plots, it does indeed look like older families have higher species richness than younger families
Again, this could simply be because there was more time to accumulate species richness.
Older families also tend to have smaller body sizes than younger families
These patterns make sense in the context of evolutionary history. Life started out microscopically and body size generally increased over evolutionary time alongside complexity. Therefore the trends we observed in chordates vs. nonchordates may just be related to the age of these groups.
Similarly, we can test hypotheses for specific groups. Following the Social/Machiavellian Hypothesis, increased brain size and capacity observed in mammal species like Homo sapiens was selected for to optimize performance in larger social groups. This is one of many competing hypotheses. Weaning age is related to brain size as altricial species with extended ontogeny focus early life resources on brain development. Create a plot and generate a linear model summary to test this hypothesis.
Call:
lm(formula = log10(primates$group_size) ~ log10(primates$weaning_age))
Residuals:
Min 1Q Median 3Q Max
-1.5870 -0.2215 0.1008 0.2775 0.8745
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.8864 0.4675 -1.896 0.064244 .
log10(primates$weaning_age) 0.8144 0.1915 4.253 0.000102 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4716 on 46 degrees of freedom
Multiple R-squared: 0.2822, Adjusted R-squared: 0.2666
F-statistic: 18.09 on 1 and 46 DF, p-value: 0.0001024
Do you find support for this hypothesis?
We find some support for the Social/Machiavellian Hypothesis with a statistically significant, positive relationship between group size and weaning age.
The large range of residuals may imply outliers should be investigated to see if there is support for their removal. This may improve the fit following the R-squared value.
Introduction to Phylogenetic Diversity Analyses
In this section, based on the online workshop by Kembel(Kembel 2024), we will explore phylogeny in the context of grassland plant biodiversity in Alberta, Canada. We will examine their functional traits and phylogenetic relationships. This analysis is based on the data described in Kembel and Cahill (2011), focusing on the independent evolution of leaf and root traits among temperate grassland plant communities.
Load the grassland.community.csv dataset as the community matrix comm with header = TRUE and row.names = 1. Do the same for the metadata, loading it into the metadata df. Explore the dataset and its metadata.
comm <-read.csv("Lab_9_data/grassland.community.csv", header =TRUE, row.names =1)head(comm)
Every cell indicates the percentage coverage of a particular species within a sample.
It’s important to note that several multivariate analytical methods, like vegan, are influenced by the overall abundance present in a sample. Therefore, it may be advisable to transform these absolute abundance measurements into relative abundance values. This adjustment can be efficiently carried out using a specific function available in the vegan package:
comm <-decostand(comm, method ="total")apply(comm, 1, sum)
Trait data include information about the leaf and root traits of each species, which are crucial for understanding functional diversity within the community. We’ll load a phylogeny for these species as well, which helps in understanding the evolutionary relationships and phylogenetic diversity within the community. Load in trait data as the df traits from the species.traits.csv file and the phylogenetic data as the df phy from the Newick file grassland.phylogeny.newick
Understanding Phylogenetic Dataset Formats
Phylogenetic datasets come in various formats, each with its specific structure and use case. Familiarizing yourself with these formats is useful for efficiently managing and analyzing phylogenetic data. Here are some of the most common types:
Newick: A compact way to represent phylogenetic trees using parentheses and commas. Despite its simplicity, Newick format captures the tree structure, including branch lengths. It doesn’t support metadata like taxa descriptions or geographical information.
NEXUS: A more versatile format than Newick, NEXUS files can contain multiple blocks of data including trees, sequences, sets of characters, and assumptions about the data. It’s widely used in software packages for evolutionary analysis because of its flexibility and the ability to include rich metadata.
– .tre: Often just a file extension for trees saved in either Newick or NEXUS formats. By itself, “.tre” doesn’t specify a format but is commonly used to denote phylogenetic trees.
– The commonly used Newick or Nexus format files can be imported into R with the read.tree or read.nexus functions
Phylip: Originally designed for the PHYLIP software package, this format is primarily used for sequence data but can also represent phylogenetic trees. It has strict formatting rules for sequence names and alignment, which can be a limitation.
FASTA: While primarily known for representing sequences, FASTA files are sometimes used in phylogenetic analysis to supply the sequence data from which trees are inferred. FASTA files are simple, with each sequence prefaced by a ‘>’ symbol followed by an identifier line.
Each format has its advantages and is supported by different software tools for phylogenetic analysis.
Our phylogenetic data structure is encapsulated in a phylo object type. Comprehensive documentation on the phylo format is accessible on the ape package’s official site (ape homepage). Intrinsically, a phylo object acts as a specialized list within R, incorporating various attributes such as labels for tips and the lengths of branches. This structured format, as delineated by the ape package, enables R to adeptly perform summarization and generate visual representations of phylogenetic trees.
Among other packages, the picante package is commonly used for phylogenetic biodiversity analyses. The function match.phylo.comm processes both a community matrix and a phylogenetic tree (phylo object), identifying species discrepancies between the two datasets. It then realigns both datasets to ensure species order consistency and mutual presence, yielding modified versions of each dataset that share identical species compositions.
It’s also essential to verify that the arrangement of our community dataset aligns with that of our metadata.
all.equal(rownames(comm), rownames(metadata))
[1] TRUE
If they didn’t match, we could have used our metadata to match the orders:
metadata <- metadata[rownames(comm), ]
Now visualize species richness (alpha diversity) between metadata habitats as a boxplot. The habitat names are in the metadata df while the counts are in our comm community data matrix. Hint: Do you remember the species richness function from vegan?
ggplot(metadata, aes(x = habitat, y =specnumber(comm))) +geom_boxplot() +labs(x ="Habitat", y ="Species Richness", title ="Species richness by habitat")
What observations do you make?
There is higher alpha diversity in fescue habitats
Phylogenetic Diversity
One of the first measures of phylogenetic relatedness within ecological communities was the phylogenetic diversity (PD) index. Faith’s PD(Faith 1992) is defined as the total branch length of a tree including all local community species, optionally also including the root node. The pd() function returns Faith’s PD and species richness (SR).
ggplot(metadata, aes(x = habitat, y = comm.pd$PD)) +geom_boxplot() +labs(x ="Habitat", y ="Faith's Phylogenetic Diversity (PD)")
Hierarchical Clustering
Like we did in our biodiversity metrics lab, let’s create a distance matrix comm.bc.dist using the vegdist() function. Since the Bray-Curtis distance incorporates species abundances, we should use abundance information when calculating our distances from our comm community data matrix.
comm.bc.dist <-vegdist(comm, method ="bray")
Similar to a phylogeny tree structure, we can plot these results as a dendrogram. We will then cluster together communities that are similar using an agglomerative hierarchical clustering algorithm:
comm.bc.clust <-hclust(comm.bc.dist, method ="average")plot(comm.bc.clust, main ="Hierarchical Clustering of Communities")
What do you observe?
It looks like mixedgrass and fescue habitats contain different plant communities
The two main clusters separate fescue samples from the other samples
We will use these distance matrices we produced later
Phylogenetic Beta-Diversity
To assess community dissimilarity in a phylogenetic context, we employ metrics such as UniFrac and PhyloSor for phylogenetic beta-diversity, alongside comdist and comdistnt for quantifying mean pairwise distance (MPD) and mean nearest taxon distance (MNTD) between communities. These metrics illuminate the evolutionary and functional relationships driving community assembly across landscapes.
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = comm.bc.dist, ydis = comm.mntd.traits.dist)
Mantel statistic r: 0.9499
Significance: 0.001
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.082 0.112 0.152 0.199
Permutation: free
Number of permutations: 999
What do you observe?
The phylogenetic MNTD and community composition show a strong positive correlation (Mantel statistic r = 0.8597), indicating that communities with similar species compositions tend to be more closely related phylogenetically. This strong association suggests that phylogenetic relatedness is a key factor in community assembly.
The significant correlations (p < 0.001) in both tests underline the importance of evolutionary histories and functional adaptations in the assembly of ecological communities. These patterns reflect the influence of both phylogenetic lineage and functional traits on the ecological strategies and environmental adaptations of species within communities.
The correlation between functional trait MNTD and community composition is even stronger (Mantel statistic r = 0.9499), emphasizing the crucial role of functional traits in structuring communities. This implies that species within similar communities not only share a closer phylogenetic relationship but also possess similar functional traits.
These findings highlight the complex interplay between phylogenetic relatedness and functional trait similarities in influencing community structure and suggest a convergence in ecological strategies among closely related and functionally similar species. Big picture, traits in this case are a good proxy for phylogeny and vice versa.
Phylogeny/Trait-Based Ordinations
Remember that ordination methods, such as Non-metric Multidimensional Scaling (NMDS) that we used in the biodiversity metrics lab, provide a powerful tool for visualizing complex patterns of similarity and dissimilarity among communities. These methods can be applied to phylogenetic metrics as well. By ordinating communities, we can discern underlying gradients of environmental variation, biotic interactions, or historical contingency shaping community structure.
Like in the biodiversity metrics lab, we will use the metaMDS() function from the vegan package
comm.mntd.mds <-metaMDS(comm.mntd.dist)
Run 0 stress 0.116569
Run 1 stress 0.116569
... Procrustes: rmse 5.247208e-06 max resid 1.302183e-05
... Similar to previous best
Run 2 stress 0.116569
... Procrustes: rmse 5.51446e-06 max resid 1.326615e-05
... Similar to previous best
Run 3 stress 0.133858
Run 4 stress 0.1355939
Run 5 stress 0.1166597
... Procrustes: rmse 0.0532475 max resid 0.2560234
Run 6 stress 0.1171039
Run 7 stress 0.116569
... Procrustes: rmse 4.690297e-06 max resid 1.212965e-05
... Similar to previous best
Run 8 stress 0.1189597
Run 9 stress 0.1166597
... Procrustes: rmse 0.05324795 max resid 0.2560256
Run 10 stress 0.1431935
Run 11 stress 0.1290987
Run 12 stress 0.1278436
Run 13 stress 0.116569
... Procrustes: rmse 4.525123e-06 max resid 1.578564e-05
... Similar to previous best
Run 14 stress 0.1171553
Run 15 stress 0.1269754
Run 16 stress 0.116569
... New best solution
... Procrustes: rmse 2.182259e-06 max resid 4.392028e-06
... Similar to previous best
Run 17 stress 0.1168237
... Procrustes: rmse 0.05400601 max resid 0.2564189
Run 18 stress 0.1278301
Run 19 stress 0.1278301
Run 20 stress 0.1269811
*** Best solution repeated 1 times
ggplot(as.data.frame(scores(comm.mntd.mds)), aes(NMDS1, NMDS2, color = metadata$habitat)) +geom_point() +labs(title ="NMDS Ordination based on Phylogenetic Distances")
What do you observe?
Visualizing trait evolution
Because we have phylogenetic data, we can also visualize evolution of traits from our dataset. Let’s visualize the evolution of leaf area:
The sections presented above offer a very introductory exploration into the integration of phylogenetic perspectives within biodiversity studies. By incorporating phylogenetic diversity, functional trait analyses, and beta-diversity metrics, we begin to appreciate the complexity and depth of ecological relationships and evolutionary histories that underpin community assembly and species distribution. These integrations not only enrich our understanding of biodiversity patterns but also enhance the precision of conservation strategies, highlighting areas of evolutionary significance and functional diversity critical for ecosystem resilience and stability.
References
Etienne, Rampal S, Sara N de Visser, Thijs Janzen, Jeanine L Olsen, Han Olff, and James Rosindell. 2012. “Can Clade Age Alone Explain the Relationship Between Body Size and Diversity?”Interface Focus 2 (2): 170–79.
Faith, Daniel P. 1992. “Conservation Evaluation and Phylogenetic Diversity.”Biological Conservation 61 (1): 1–10.
Jones, Kate E, Jon Bielby, Marcel Cardillo, Susanne A Fritz, Justin O’Dell, C David L Orme, Kamran Safi, et al. 2009. “PanTHERIA: A Species-Level Database of Life History, Ecology, and Geography of Extant and Recently Extinct Mammals: Ecological Archives E090-184.”Ecology 90 (9): 2648–48.
Kembel, Steven. 2024. “R for Biodiversity Analysis Workshop.” Université du Québec à Montréal (UQÀM); https://kembellab.ca/r-workshop/.
Mittelbach, Gary G, and Brian J McGill. 2019. Community Ecology. Oxford University Press.
