Tutorial for analyses done in ‘Bridging performance and adaptive landscapes to understand long-term functional evolution’
Monique N. Simon and Daniel S. Moen
04 April 2023
- R version and loading packages
- Part I: Performance Landscape Approach (PLA)
- Part II: Adaptive Landscape Approach (ALA)
- Part III: Merged framework
- Estimating adaptive peaks and 95% confidence intervals using bootstrapping
- Comparing performance and adaptive peaks
- Comparing species phenotypes with adaptive peaks
- Testing for the effect of historical versus species-specific factors on phenotypic variation
- Testing the effect of time spent in current environment on species’ deviations from the adaptive peak
- Conclusion
- References
In this tutorial, you can reproduce the analyses performed in the paper ‘Bridging performance and adaptive landscapes to understand long-term functional evolution’. Our main aim was to address whether aquatic and terrestrial lifestyles have influenced the location of performance and adaptive peaks for turtle shells.
We re-analyzed turtle shell data collected by Stayton et al. (2018) for 237 species, the same data also used in Stayton (2019). We performed first the performance landscape approach (PLA, Part I), then the adaptive landscape approach (ALA, Part II), and finally the merged framework: comparing performance and adaptive peaks, then comparing species’ phenotypes with inferred primary adaptive optima (Part III).
R version and loading packages
R.Version()$version.string## [1] "R version 4.2.0 (2022-04-22)"R.Version()$platform## [1] "x86_64-apple-darwin17.0"library(metR)
packageVersion("metR")## [1] '0.14.0'library(geomorph)
packageVersion("geomorph")## [1] '4.0.5'library(phytools)
packageVersion("phytools")## [1] '1.2.0'library(ape)
packageVersion("ape")## [1] '5.6.2'library(tidyverse)
packageVersion("tidyverse")## [1] '1.3.2'library(ouch)
packageVersion("ouch")## [1] '2.18'library(ggplot2)
packageVersion("ggplot2")## [1] '3.4.0'library(viridis)
packageVersion("viridis")## [1] '0.6.2'library(phylolm)
packageVersion("phylolm")## [1] '2.6.2'Part I: Performance Landscape Approach (PLA)
Estimating the performance landscapes for stress, drag, and self-righting ability
We will first read data provided by Stayton (2019) to construct a theoretical morphospace using the principal component (PC) axes that explain the most variance in the shape of turtle shells. Then, we will plot the performance surface for each function known to be relevant in turtle shells: breaking stress, drag, and self-righting ability. Stayton (2019) determined the association between PC scores and function using modeling and experiments.
stayton_data <- read.csv("stayton_data.csv", header = TRUE)
colnames(stayton_data) <- c("PC1", "PC2", "stress", "drag", "self.righting")
head(stayton_data)## PC1 PC2 stress drag self.righting
## 1 -0.15 -0.150 276.10 0.19815 3.6765
## 2 -0.15 -0.125 255.17 0.18320 3.6275
## 3 -0.15 -0.100 239.30 0.17985 3.5202
## 4 -0.15 -0.075 222.68 0.18298 3.4146
## 5 -0.15 -0.050 209.36 0.19181 3.2992
## 6 -0.15 -0.025 193.85 0.20483 3.1750plot(stayton_data$PC1, stayton_data$PC2, xlab = "PC1", ylab = "PC2", main = "Theoretical morphospace")
To plot the performance surfaces, we will use the package metR that can create contour plots using the ggplot framework.
stress_contour <- ggplot(stayton_data, aes(PC1, PC2, z = stress)) + geom_contour_fill(breaks = MakeBreaks(binwidth = 15,
exclude = 0)) + scale_fill_divergent(breaks = MakeBreaks(binwidth = 15,
exclude = 0)) + coord_equal() + theme_bw() + scale_fill_viridis() + guides(fill = guide_legend(title = "Stress")) +
theme(text = element_text(size = 13, color = "black"), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.text = element_text(size = 12,
color = "black"), axis.title = element_text(size = 16, color = "black"),
legend.text = element_text(size = 10), legend.title = element_text(size = 12))
stress_contour
drag_contour <- ggplot(stayton_data, aes(PC1, PC2, z = drag)) + geom_contour_fill() +
coord_equal() + theme_bw() + scale_fill_viridis() + guides(fill = guide_legend(title = "Drag")) +
theme(text = element_text(size = 13, color = "black"), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.text = element_text(size = 12,
color = "black"), axis.title = element_text(size = 16, color = "black"),
legend.text = element_text(size = 10), legend.title = element_text(size = 12))
drag_contour
righting_contour <- ggplot(stayton_data, aes(PC1, PC2, z = self.righting)) +
geom_contour_fill() + coord_equal() + theme_bw() + scale_fill_viridis() +
guides(fill = guide_legend(title = "Self righting")) + theme(text = element_text(size = 13,
color = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text = element_text(size = 12, color = "black"), axis.title = element_text(size = 16,
color = "black"), legend.text = element_text(size = 10), legend.title = element_text(size = 12))
righting_contour
With these performance surfaces, we can see the mapping of the shell phenotypes (PC scores) to performance on each of the three functions. Note that hotter colors actually indicate lower performance because they reflect higher stress, higher drag, and longer self-righting time. We can also confirm that these surfaces are the same as in Stayton (2019).
Estimating the combined performance landscapes for aquatic and terrestrial species
We will now to estimate the combined performance landscape for terrestrial and aquatic species. To do this, we will first re-scale performance to range from 0 (minimum performance) to 1 (maximum performance), as described in Polly et al. (2016). We will also subtract from 1, such that greater values indicate better performance (hotter colors) in the combined landscapes (Stayton 2019). Then, we will use the average weights for stress, drag, and self-righting ability calculated for aquatic and terrestrial species, as reported in Stayton (2019).
str(stayton_data)## 'data.frame': 255 obs. of 5 variables:
## $ PC1 : num -0.15 -0.15 -0.15 -0.15 -0.15 -0.15 -0.15 -0.15 -0.15 -0.15 ...
## $ PC2 : num -0.15 -0.125 -0.1 -0.075 -0.05 -0.025 0 0.025 0.05 0.075 ...
## $ stress : num 276 255 239 223 209 ...
## $ drag : num 0.198 0.183 0.18 0.183 0.192 ...
## $ self.righting: num 3.68 3.63 3.52 3.41 3.3 ...perf_rescaled <- as.data.frame(apply(stayton_data[, 3:5], 2, function(x) 1 -
((x - min(x))/(max(x) - min(x)))))
head(perf_rescaled)## stress drag self.righting
## 1 0.5580185 0.9822614 0.00000000
## 2 0.6252570 0.9967528 0.01971672
## 3 0.6762400 1.0000000 0.06289232
## 4 0.7296325 0.9969660 0.10538387
## 5 0.7724235 0.9884069 0.15181877
## 6 0.8222501 0.9757864 0.20179462str(perf_rescaled)## 'data.frame': 255 obs. of 3 variables:
## $ stress : num 0.558 0.625 0.676 0.73 0.772 ...
## $ drag : num 0.982 0.997 1 0.997 0.988 ...
## $ self.righting: num 0 0.0197 0.0629 0.1054 0.1518 ...perf_aquatic <- (0.47 * perf_rescaled$stress) + (0.38 * perf_rescaled$drag) +
(0.15 * perf_rescaled$self.righting)
combined_aquatic <- cbind(stayton_data[, 1:2], perf_aquatic)
perf_terrestrial <- (0.67 * perf_rescaled$stress) + (0.18 * perf_rescaled$drag) +
(0.16 * perf_rescaled$self.righting)
# Combine datasets
all_perf <- cbind(stayton_data[, 1:2], perf_aquatic, perf_terrestrial) %>%
as_tibble() %>%
pivot_longer(cols = c(perf_aquatic, perf_terrestrial), names_to = "type",
values_to = "Performance") %>%
separate(col = "type", into = c(NA, "type"), sep = "_") %>%
mutate(type = factor(type), type = fct_recode(type, `A. Combined aquatic landscape` = "aquatic",
`B. Combined terrestrial landscape` = "terrestrial"))
all_perf## # A tibble: 510 × 4
## PC1 PC2 type Performance
## <dbl> <dbl> <fct> <dbl>
## 1 -0.15 -0.15 A. Combined aquatic landscape 0.636
## 2 -0.15 -0.15 B. Combined terrestrial landscape 0.551
## 3 -0.15 -0.125 A. Combined aquatic landscape 0.676
## 4 -0.15 -0.125 B. Combined terrestrial landscape 0.601
## 5 -0.15 -0.1 A. Combined aquatic landscape 0.707
## 6 -0.15 -0.1 B. Combined terrestrial landscape 0.643
## 7 -0.15 -0.075 A. Combined aquatic landscape 0.738
## 8 -0.15 -0.075 B. Combined terrestrial landscape 0.685
## 9 -0.15 -0.05 A. Combined aquatic landscape 0.761
## 10 -0.15 -0.05 B. Combined terrestrial landscape 0.720
## # … with 500 more rows# Color breaks
min_perf <- round(min(all_perf$Performance), 2)
max_perf <- 1
mybreaks <- seq.int(min_perf, max_perf, length.out = 12)
# Plot together with facetting to put on same color scale
all_perf %>%
ggplot(aes(PC1, PC2, z = Performance)) + geom_contour_filled(breaks = mybreaks) +
coord_equal() + theme_bw() + guides(fill = guide_legend(title = "Performance")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text = element_text(size = 14, color = "black"), axis.title = element_text(size = 16,
color = "black"), legend.text = element_text(size = 8), legend.title = element_text(size = 10),
strip.text = element_text(size = 12, hjust = 0), strip.background = element_blank()) +
facet_wrap(~type)
Estimating the performance peaks for aquatic and terrestrial species
With the combined landscapes, we can locate the performance peaks, that is, the PC scores that correspond to the highest performance. We will also calculate a region of high performance that is up to 5% lower than peak performance, to later compare with the 95% confidence intervals of inferred primary adaptive peaks using OU modeling (Table A1, available online).
combined_aquatic <- cbind(stayton_data[, 1:2], perf_aquatic)
colnames(combined_aquatic) <- c("PC1", "PC2", "perf_aquatic")
max_comb_aquatic <- combined_aquatic[combined_aquatic$perf_aquatic == max(combined_aquatic$perf_aquatic),
]
max_comb_aquatic #aquatic performance peak## PC1 PC2 perf_aquatic
## 82 -0.025 0 0.8619485max_comb_aquatic$perf_aquatic * 0.95## [1] 0.8188511# region of 5% lower performance
summary(combined_aquatic[which(combined_aquatic$perf_aquatic > 0.818), ])## PC1 PC2 perf_aquatic
## Min. :-0.12500 Min. :-0.02500 Min. :0.8181
## 1st Qu.:-0.07500 1st Qu.: 0.00000 1st Qu.:0.8294
## Median :-0.05000 Median : 0.02500 Median :0.8418
## Mean :-0.05403 Mean : 0.03468 Mean :0.8404
## 3rd Qu.:-0.02500 3rd Qu.: 0.06250 3rd Qu.:0.8513
## Max. : 0.02500 Max. : 0.12500 Max. :0.8619combined_terrestrial <- cbind(stayton_data[, 1:2], perf_terrestrial)
colnames(combined_terrestrial) <- c("PC1", "PC2", "perf_terrestrial")
max_comb_terrestrial <- combined_terrestrial[combined_terrestrial$perf_terrestrial ==
max(combined_terrestrial$perf_terrestrial), ]
max_comb_terrestrial #terrestrial performance peak## PC1 PC2 perf_terrestrial
## 173 0.125 0.025 0.9047935max_comb_terrestrial$perf_terrestrial * 0.95## [1] 0.8595538# region of 5% lower performance
summary(combined_terrestrial[which(combined_terrestrial$perf_terrestrial > 0.859),
])## PC1 PC2 perf_terrestrial
## Min. :-0.05000 Min. :-0.02500 Min. :0.8591
## 1st Qu.: 0.00000 1st Qu.: 0.00000 1st Qu.:0.8669
## Median : 0.07500 Median : 0.02500 Median :0.8739
## Mean : 0.06223 Mean : 0.03723 Mean :0.8754
## 3rd Qu.: 0.12500 3rd Qu.: 0.07500 3rd Qu.:0.8835
## Max. : 0.17500 Max. : 0.10000 Max. :0.9048We can see that the regions of 5% lower performance from the peak are quite large. They range from -0.125 to 0.025 for PC1 and -0.025 to 0.125 for PC2 in the aquatic combined landscape, and from -0.05 to 0.175 for PC1 and -0.025 to 0.100 for PC2 in the terrestrial combined landscape.
Part II: Adaptive Landscape Approach (ALA)
Projecting species phenotypes in shape space
We will begin with the 53 three-dimensional landmarks placed on the shells of 274 species of turtles, provided by Stayton (2019). We will then perform a General Procrustes Analysis to project each species onto the PCs that reflect variation in shell shape. These PC axes are the same as the ones used to estimate performance landscapes. Note that we downloaded the original data from Dryad but corrected some misspellings in the names of some species.
turtle_data <- read.csv("turtle_data.csv", row.names = 1)
length(levels(as.factor(turtle_data$Name)))## [1] 274# Changing the data frame to an array
turtle_array <- arrayspecs(turtle_data[, 8:166], 53, 3, sep = "")
length(turtle_data$Name)## [1] 2722dimnames(turtle_array)[[3]] <- turtle_data$Museum.Number
# had to put individuals names to get specimens' scores
head(turtle_array[, , 1])## X Y Z
## 1 -72.66066 111.0733 -48.95750
## 2 -75.59137 133.1744 -53.16478
## 3 -86.27445 151.4542 -54.64143
## 4 -98.03491 167.3384 -61.81512
## 5 -109.82220 177.7144 -69.63888
## 6 -129.79511 185.4073 -72.50737# 53 landmarks (x, y, z) for the first specimen
# Align all specimens with GPA
proc_turtle <- gpagen(turtle_array)##
## Performing GPA
##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%
##
## Making projections... Finished!# shape variables
# Perform a PCA analysis and get PCs scores of all specimens
dimnames(proc_turtle$coords)[[3]] <- turtle_data$Museum.Number
proc_turtle_pca <- gm.prcomp(proc_turtle$coords)
plot(proc_turtle_pca, main = "PC scores for all specimens of turtles")
From the plot above, we can confirm that the PCs we obtained are the same as the ones reported in Stayton (2019), because they explain the same amount of variance in shape.
We will now calculate the average scores for each species, then use the consensus phylogeny from McLaughlin and Stayton (2016) to keep only the 237 species that we analyzed because they are present in this phylogeny. To do this, we will have to match the shape data with the tip labels in the phylogeny.
# Getting average scores per species
proc_turtle_pca$x <- as.data.frame(proc_turtle_pca$x)
specimens <- turtle_data$Name
proc_turtle_pca$x$specimens <- specimens
length(proc_turtle_pca$x$specimens)## [1] 2722turtle_mean_pc <- aggregate(proc_turtle_pca$x[, 1:2], by = list(proc_turtle_pca$x$specimens),
mean)
head(turtle_mean_pc)## Group.1 Comp1 Comp2
## 1 Acanthochelys_macrocephala 0.02881497 0.010182069
## 2 Acanthochelys_pallidipectoris 0.06379932 -0.030418939
## 3 Acanthochelys_radiata 0.03855873 0.001652305
## 4 Acanthochelys_spixii 0.03111754 -0.006863926
## 5 Aldabrachelys_gigantea -0.08764407 0.041687006
## 6 Astrochelys_radiata -0.14901220 0.059577794turtle_tree <- read.tree("turtle_tree.tre")
turtle_tree #using consensus tree already pruned to have the 237 species##
## Phylogenetic tree with 237 tips and 236 internal nodes.
##
## Tip labels:
## Caretta_caretta, Lepidochelys_kempii, Lepidochelys_olivacea, Eretmochelys_imbricata, Chelonia_mydas, Natator_depressus, ...
##
## Rooted; includes branch lengths.## Need to match species data with species phylogeny
which(is.na(match(levels(as.factor(turtle_data$Name)), turtle_tree$tip.label)))## [1] 15 17 18 19 20 21 23 24 25 26 27 28 29 36 37 65 66 67 86
## [20] 95 97 112 148 164 165 166 167 168 169 174 197 203 204 227 260 271 273# species of turtle data that are not in the phylogeny (n=37)
# Now exclude the extra species to have n=237
species <- turtle_tree$tip.label
# Adding species names again now that some names in the data have been
# changed
specimens <- turtle_data$Name
proc_turtle_pca$x$specimens <- specimens
length(proc_turtle_pca$x$specimens)## [1] 2722turtle_mean_pc <- aggregate(proc_turtle_pca$x[, 1:2], by = list(proc_turtle_pca$x$specimens),
mean)
str(turtle_mean_pc) #274 species## 'data.frame': 274 obs. of 3 variables:
## $ Group.1: chr "Acanthochelys_macrocephala" "Acanthochelys_pallidipectoris" "Acanthochelys_radiata" "Acanthochelys_spixii" ...
## $ Comp1 : num 0.0288 0.0638 0.0386 0.0311 -0.0876 ...
## $ Comp2 : num 0.01018 -0.03042 0.00165 -0.00686 0.04169 ...# Removing extra specimens
proc_turtle_pca$x <- proc_turtle_pca$x[-c(which(is.na(match(proc_turtle_pca$x$specimens,
species)))), ]
# calculating means again using only specimens that correspond to the 237
# species we will analyze:
turtle_mean_pc.1 <- aggregate(proc_turtle_pca$x[, 1:2], by = list(proc_turtle_pca$x$specimens),
mean)
length(levels(as.factor(proc_turtle_pca$x$specimens))) #237## [1] 237## Multiplying all pc1 scores by -1.Signs flipped when comparing with
## Stayton (2019)
turtle_mean_pc.1$Comp1 = turtle_mean_pc.1$Comp1 * -1
# Include habitat info
turtle_data_clean <- turtle_data[-c(which(is.na(match(turtle_data$Name, species)))),
]
habitat <- (aggregate(turtle_data_clean$Habitat, by = list(turtle_data_clean$Name,
turtle_data_clean$Habitat), function(x) x))[, 1:2]
# it is alphabetic order per habitat: change to alphabetic order to match
# species
habitat.1 <- as.data.frame(habitat[order(habitat$Group.1), ])
# substitute 'A' and 'T' by colors blue and brown
habitat.1[c(which(habitat.1$Group.2 == "A")), 2] = "blue"
habitat.1[c(which(habitat.1$Group.2 == "T")), 2] = "brown"
plot(turtle_mean_pc.1$Comp1, turtle_mean_pc.1$Comp2, pch = 21, bg = habitat.1$Group.2,
xlab = "PC1 (21.9%)", ylab = "PC2 (10.6%)")
legend("topright", pch = 21, pt.bg = unique(habitat.1$Group.2), legend = c("aquatic",
"terrestrial"), box.lty = 0)
The plot above shows the average scores for the 237 species we will analyze, separated into aquatic and terrestrial habitats.
Estimating ancestral states
Now we will create a data frame with all the information we need to reconstruct ancestral states of habitat, to then use them in the OU modeling using the R package ouch. From now on, we will run analyses using several functions created in previous work by one of us. These functions are provided in the ‘R_functions_PBZ.R’ script as supplementary material. The phylogeny with the reconstructed ancestral states is presented in Fig. A1 of the online appendix.
# Joining habitat info with PC scores
habitat.1 <- as.data.frame(habitat[order(habitat$Group.1), ])
turtle_dat <- cbind(habitat.1, turtle_mean_pc.1[, 2:3])
# putting in the same order as the phylogeny
turtle_dat <- turtle_dat[match(turtle_tree$tip.label, turtle_dat$Group.1), ]
colnames(turtle_dat) <- c("Species", "Habitat", "PC1", "PC2")
# Reconstructing ancestral states with ace
head(turtle_dat)## Species Habitat PC1 PC2
## 11 Caretta_caretta A -0.06998729 -0.0006996088
## 78 Lepidochelys_kempii A -0.07129763 -0.0191977281
## 79 Lepidochelys_olivacea A -0.03459191 0.0095317479
## 41 Eretmochelys_imbricata A -0.07389981 0.0102543461
## 13 Chelonia_mydas A -0.05316871 -0.0234931753
## 100 Natator_depressus A -0.07356504 -0.0633769270states_ou <- setNames(turtle_dat$Habitat, turtle_dat$Species)
# Put the data in the same order as the phylogeny
states_ou <- states_ou[turtle_tree$tip.label]
aic_ER <- AIC(ace(states_ou, turtle_tree, type = "discrete", model = "ER"))
aic_ER## [1] 120.9997aic_ARD <- AIC(ace(states_ou, turtle_tree, type = "discrete", model = "ARD"))
aic_ARD## [1] 122.9325# Calculating AIC weights
source("R_functions_PBZ.R")
aic_list <- list(ER = aic_ER, ARD = aic_ARD)
aicc_wgt_calc(unlist(aic_list))## ER ARD
## 0.7244039 0.2755961# function to calculate AICc weights
anc_states <- ace(states_ou, turtle_tree, type = "discrete", model = "ER")
anc_states$rates## [1] 0.005006582head(anc_states$lik.anc)## A T
## 238 0.7737137 0.2262863052
## 239 0.8900989 0.1099010564
## 240 0.8902528 0.1097471985
## 241 0.9942198 0.0057802000
## 242 0.9941277 0.0058723000
## 243 0.9997557 0.0002442976# use this as internal states
# vector of ancestral states (get states with higher support for each
# node)
turtle_tree$Nnode## [1] 236x = anc_states$lik.anc[, 1] > 0.51 #higher support for A (TRUE)
length(anc_states$lik.anc[, 1][x]) #178 nodes aquatic## [1] 178y = anc_states$lik.anc[, 2] > 0.51 #higher support for T (TRUE)
length(anc_states$lik.anc[, 2][y]) #58 nodes terrestrial## [1] 58head(anc_states$lik.anc[, 1][x])## 238 239 240 241 242 243
## 0.7737137 0.8900989 0.8902528 0.9942198 0.9941277 0.9997557# names are nodes and numbers are support for A replace support for letter
# A
str(anc_states$lik.anc[, 1][x])## Named num [1:178] 0.774 0.89 0.89 0.994 0.994 ...
## - attr(*, "names")= chr [1:178] "238" "239" "240" "241" ...as.numeric(names(anc_states$lik.anc[, 1][x]))## [1] 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
## [19] 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
## [37] 274 275 276 277 278 279 280 285 286 287 288 289 290 291 292 293 294 295
## [55] 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313
## [73] 314 315 316 318 319 320 321 322 323 324 325 326 327 328 329 335 336 337
## [91] 381 382 383 384 385 386 388 393 394 395 396 397 398 399 400 401 402 403
## [109] 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421
## [127] 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
## [145] 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457
## [163] 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473nodes.1 = c(rep("A", 178))
nodes.2 = c(rep("T", 58))
nodes_states <- c(nodes.1, nodes.2)
nodes_names <- c(as.numeric(names(anc_states$lik.anc[, 1][x])), as.numeric(names(anc_states$lik.anc[,
2][y])))
node_states <- nodes_states[order(nodes_names)]
## Internal states in the phylogeny
turtle_tree$node.label = node_states
turtle_tree##
## Phylogenetic tree with 237 tips and 236 internal nodes.
##
## Tip labels:
## Caretta_caretta, Lepidochelys_kempii, Lepidochelys_olivacea, Eretmochelys_imbricata, Chelonia_mydas, Natator_depressus, ...
## Node labels:
## A, A, A, A, A, A, ...
##
## Rooted; includes branch lengths.We can now see that the object turtle_tree has node
labels (A for aquatic, and T for terrestrial).
Comparing models of phenotypic evolution
The next step is to set up the data frame to be used in ouch, so that we can compare different multivariate models: Brownian motion, single optimum Ornstein-Uhlenbeck (OU), and 2-peak OU (aquatic and terrestrial).
## Set up data for ouch
source("R_functions_PBZ.R")
turtle_ouch_OU2 <- data.frame(turtle_dat[, 2:4])
rownames(turtle_ouch_OU2) <- turtle_dat$Species
# Now we add the phylogeny and internal-state estimates
head(rownames(turtle_ouch_OU2) == turtle_tree$tip.label) #true## [1] TRUE TRUE TRUE TRUE TRUE TRUEturtle_ouch_OU2 <- divers_to_ouch(turtle_tree, turtle_tree$node.label, turtle_ouch_OU2) #function to construct the data frame to use in 'ouch' analyses## Loading required package: geigermax(node.depth.edgelength(turtle_tree)) #length of phylogeny = 227.77 mya## [1] 227.7755# Plotting the internal states
ouch_tree <- with(turtle_ouch_OU2, ouchtree(nodes, ancestors, times, labels))
plot(ouch_tree, regimes = turtle_ouch_OU2["model"], palette = c("blue", "brown"),
margin = 0.15, labels = TRUE)
# different colors represent different selective regimes
# Set ouchtree for all analyses
ouch_tree <- with(turtle_ouch_OU2, ouchtree(nodes, ancestors, times, labels))
## OUCH analysis: comparing multivariate models (a) Brownian motion for
## PC1 and PC2 scores at once
brown_ouch_turtle <- brown(turtle_ouch_OU2[c("PC1", "PC2")], ouch_tree)
# (b) Single-optimum OU
set.seed(10) #setting seed so that analysis can be reproduced
turtle_ouch_OU2$global <- as.factor("global") # Set single regime for OU1
OU1_ouch_turtle <- hansen(turtle_ouch_OU2[c("PC1", "PC2")], ouch_tree, regimes = turtle_ouch_OU2["global"],
sqrt.alpha = c(1, 0, 1), sigma = c(1, 0, 1), maxit = 1000)
# (c) Two optima OU2: changed sigma values to closer to 0 to converge
# better
OU2_ouch_turtle <- hansen(turtle_ouch_OU2[c("PC1", "PC2")], ouch_tree, regimes = turtle_ouch_OU2["model"],
sqrt.alpha = c(2, 0, 2), sigma = c(0.001, 0, 0.001), maxit = 1000)
# Put output objects into a list to more easily summarize results
all_out <- list(brown_ouch_turtle, OU1_ouch_turtle, OU2_ouch_turtle)
names(all_out) <- c("BM", "OU1", "OU2")
# Summarize model fits
ouch_sum <- map(all_out, summary)
turtle_fits <- tibble(model = names(all_out), k = map_dbl(ouch_sum, ~.$dof),
loglik = map_dbl(ouch_sum, ~.$loglik), AICc = map_dbl(ouch_sum, ~.$aic.c),
wgt = aicc_wgt_calc(AICc))
turtle_fits## # A tibble: 3 × 5
## model k loglik AICc wgt
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 BM 5 1120. -2229. 3.54e-12
## 2 OU1 8 1139. -2262. 3.97e- 5
## 3 OU2 10 1151. -2282. 1.00e+ 0From the table above, we can conclude that the most-supported model is the 2-peak OU model, describing two divergent optima for aquatic and terrestrial turtles. We can now move forward and apply the merged framework to the turtle shell phenotypes.
Part III: Merged framework
Estimating adaptive peaks and 95% confidence intervals using bootstrapping
We will first calculate 95% confidence intervals for the OU model parameters we just inferred for aquatic and terrestrial species. We can do this by performing a parametric bootstrap analysis using the ouch package.
# adaptive peaks for aquatic and terrestrial species
ouch_sum$OU2$optima## $PC1
## A T
## -0.04126470 0.07602188
##
## $PC2
## A T
## 0.008458624 -0.013919114boot_OU2 <- bootstrap(OU2_ouch_turtle, nboot = 3)
# increase nboot=1000 to get similar results as the paper (Table A1),
# minus slight differences due to the random aspects of simulation
boot_CI <- t(sapply(boot_OU2, quantile, probs = c(0.025, 0.975)))
head(boot_CI)## 2.5% 97.5%
## alpha1 3.656151957 4.915617509
## alpha2 -0.597278119 0.936948869
## alpha3 -0.597278119 0.936948869
## alpha4 4.573345791 5.618456248
## sigma.squared1 0.013013418 0.016545763
## sigma.squared2 0.006863161 0.007460083Comparing performance and adaptive peaks
Now we can compare the locations of performance peaks (point estimate and region of 5% lower highest performance) and the 95% CI of adaptive peaks. This comparison is shown in Table A1 when running 1,000 bootstrap iterations and in Figure 2 of our paper.
# 95% CI of aquatic adaptive peak
boot_CI[c(9, 11), ]## 2.5% 97.5%
## optima.PC1.A -0.053574777 -0.050095505
## optima.PC2.A 0.003521865 0.005573743# 95% CI of terrestrial adaptive peak
boot_CI[c(10, 12), ]## 2.5% 97.5%
## optima.PC1.T 0.03276136 0.067258335
## optima.PC2.T -0.06689918 -0.007268892# point estimates of performance peaks
max_comb_aquatic## PC1 PC2 perf_aquatic
## 82 -0.025 0 0.8619485max_comb_terrestrial## PC1 PC2 perf_terrestrial
## 173 0.125 0.025 0.9047935# min and max PC scores of the region with 5% lower performance
summary(combined_aquatic[which(combined_aquatic$perf_aquatic > 0.818), ])## PC1 PC2 perf_aquatic
## Min. :-0.12500 Min. :-0.02500 Min. :0.8181
## 1st Qu.:-0.07500 1st Qu.: 0.00000 1st Qu.:0.8294
## Median :-0.05000 Median : 0.02500 Median :0.8418
## Mean :-0.05403 Mean : 0.03468 Mean :0.8404
## 3rd Qu.:-0.02500 3rd Qu.: 0.06250 3rd Qu.:0.8513
## Max. : 0.02500 Max. : 0.12500 Max. :0.8619summary(combined_terrestrial[which(combined_terrestrial$perf_terrestrial > 0.859),
])## PC1 PC2 perf_terrestrial
## Min. :-0.05000 Min. :-0.02500 Min. :0.8591
## 1st Qu.: 0.00000 1st Qu.: 0.00000 1st Qu.:0.8669
## Median : 0.07500 Median : 0.02500 Median :0.8739
## Mean : 0.06223 Mean : 0.03723 Mean :0.8754
## 3rd Qu.: 0.12500 3rd Qu.: 0.07500 3rd Qu.:0.8835
## Max. : 0.17500 Max. : 0.10000 Max. :0.9048With the 1,000 iterations, we would see that the point estimate for the aquatic performance peak matches the 95% CI of the adaptive peak, but this is not the case for the terrestrial performance peak. However, when considering the region of 5% lower performance from the peak, the performance and adaptive peaks do overlap for terrestrial species.
Comparing species phenotypes with adaptive peaks
The next step in the merged framework is to calculate the distances of species’ phenotypes to the inferred adaptive peaks. We will first reproduce Figure 3A here to plot all species’ phenotypes and the aquatic and terrestrial adaptive peaks.
Figure 3A: Species’ phenotypes and adaptive peaks
Fig.3A <- ggplot(turtle_dat, aes(x = PC1, y = PC2, fill = Habitat, shape = Habitat)) +
geom_point(size = 2, color = "black") + scale_shape_manual(values = c(21,
24)) + scale_fill_manual(values = c("gray70", "black")) + scale_x_continuous(name = "PC1 (21.9%)",
breaks = seq(-0.1, 0.15, 0.05)) + scale_y_continuous(name = "PC2 (10.6%)",
breaks = seq(-0.1, 0.15, 0.05)) + coord_fixed() + theme_bw() + theme(text = element_text(size = 13,
color = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text = element_text(size = 12, color = "black"), axis.title = element_text(size = 16,
color = "black"), legend.position = c(0.9, 0.85), legend.text = element_text(size = 10),
legend.title = element_text(size = 12), legend.background = element_rect(fill = "white",
color = "black"))
# Adding the adaptive peaks
PC1 <- c(-0.041, 0.076)
PC2 <- c(0.008, -0.0139)
Habitat <- c("A", "T", "A", "T")
opt_df <- data.frame(PC1, PC2, Habitat)
Fig.3A_final <- Fig.3A + geom_point(data = opt_df, aes(x = PC1, y = PC2, fill = Habitat),
shape = 21, size = 6)
Fig.3A_final
From the plot above, we can see that several species’ phenotypes are far from the adaptive peak of their current environment. We can also see that some terrestrial species are closer to the aquatic adaptive peak than to their own peak. This suggests that some terrestrial lineages have only recently changed from an aquatic ancestral environment to the terrestrial environment (as can be seen in Fig. A2)
Testing for the effect of historical versus species-specific factors on phenotypic variation
We can now test for the effects of historical factors versus species-specific factors on phenotypic variation. We expect a high proportion of variation explained by history for lineages that changed environments more recently (i.e., terrestrial lineages).
We will perform this test separately for aquatic and terrestrial species, first including the family Testudinidae and then excluding this family. Results of these tests are shown in Table 2 of our paper.
source("R_functions_PBZ.R")
turtle_optima <- t(as.data.frame(ouch_sum$OU2$optima))
turtle_optima## A T
## PC1 -0.041264703 0.07602188
## PC2 0.008458624 -0.01391911# Separating aquatic and terrestrial lineages aquatic:
aquatic_dat <- turtle_dat[turtle_dat$Habitat == "A", 2:4]
aquatic_opt <- as.data.frame(turtle_optima[, 1])
colnames(aquatic_opt) <- "A"
hist_test_A <- history.test(aquatic_dat, aquatic_opt, c("PC1", "PC2"))
# function to run the test of historical vs species-specific factors
hist_test_A## hist.var specific.var total.var prop.hist prop.specific
## PC1 2.519244e-04 0.0009823335 0.001234258 0.20411003 0.7958900
## PC2 5.840659e-05 0.0019059855 0.001964392 0.02973265 0.9702673
## Total 3.103310e-04 0.0028883191 0.003198650 0.09701937 0.9029806# terrestrial
terrestrial_dat <- turtle_dat[turtle_dat$Habitat == "T", 2:4]
terrestrial_opt <- as.data.frame(turtle_optima[, 2])
colnames(terrestrial_opt) <- "T"
hist_test_T <- history.test(terrestrial_dat, terrestrial_opt, c("PC1", "PC2"))
hist_test_T## hist.var specific.var total.var prop.hist prop.specific
## PC1 -9.253154e-06 0.002798929 0.002789676 -0.003316928 1.0033169
## PC2 1.567038e-04 0.000878833 0.001035537 0.151326170 0.8486738
## Total 1.474507e-04 0.003677762 0.003825213 0.038547048 0.9614530# terrestrial excluding Testudinidae
turtle_dat.1 <- turtle_dat[-c(115, 99, 100, 109, 105, 106, 110, 111, 104, 107,
108, 137:140, 102, 103, 114, 124, 129:131, 116:123, 136, 141, 142, 125:127,
112, 113, 101, 128, 132:135), ]
terrestrial_dat.1 <- turtle_dat.1[turtle_dat.1$Habitat == "T", 2:4] #n=24
hist_test_T.1 <- history.test(terrestrial_dat.1, terrestrial_opt, c("PC1", "PC2"))
hist_test_T.1## hist.var specific.var total.var prop.hist prop.specific
## PC1 3.901390e-03 0.0007217576 0.0046231475 0.84388177 0.1561182
## PC2 -3.696362e-05 0.0009263081 0.0008893445 -0.04156277 1.0415628
## Total 3.864426e-03 0.0016480658 0.0055124920 0.70103072 0.2989693From the results above, we can see that the proportion of phenotypic variation explained by history is only high for terrestrial species when removing Testudinidae. This is because Testudinidae lineages have spent more time in their current environment compared with other terrestrial lineages.
For aquatic lineages, much more phenotypic variation is explained by species-specific factors, such as other (untested) ecological factors and genetic constraints.
Testing the effect of time spent in current environment on species’ deviations from the adaptive peak
We can now specifically test if lineages that have spent less time in their current environment show a higher deviation from their adaptive peak. To do this, we will perform a phylogenetic generalized least-squares (PGLS) regression of distance from optimum on time spent in current environment.
## Make node list (i.e., nodes where changes in microhabitat occurred;
## determined in advance) Note that these are class 'phylo' node orders,
## NOT 'ouch' order Jump intermediary nodes and associate descendants
## arising from the node that changed Node 274 that originates
## Testudnidae: use Aldabrachelys, Astrochelys, Centrochelys, Chelonoidis,
## Chersina, Geochelone, Gopherus, Homopus, Indotestudo, Kinixys,
## Malachocersus, Manouria, Psammobates, Pyxis, Stigmochelys, Testudo
nodes <- c(286, 280, 298, 320, 329, 327, 315, 334, 275, 274, 388, 390, 384,
450) #13 changes of environment
name.list <- list("Cuora_mouhotii", c("Cuora_mccordi", "Cuora_picturata", "Cuora_bourreti",
"Cuora_galbinifrons", "Cuora_flavomarginata"), "Vijayachelys_silvatica",
"Leucocephalon_yuwonoi", "Heosemys_depressa", "Heosemys_spinosa", c("Geoemyda_japonica",
"Geoemyda_spengleri"), c("Rhinoclemmys_funerea", "Rhinoclemmys_melanosterna",
"Rhinoclemmys_punctularia", "Rhinoclemmys_diademata"), c("Rhinoclemmys_areolata",
"Rhinoclemmys_nasuta", "Rhinoclemmys_rubida", "Rhinoclemmys_annulata",
"Rhinoclemmys_pulcherrima"), c("Aldabrachelys", "Astrochelys", "Centrochelys",
"Chelonoidis", "Chersina", "Geochelone", "Gopherus", "Homopus", "Indotestudo",
"Kinixys", "Malacochersus", "Manouria", "Psammobates", "Pyxis", "Stigmochelys",
"Testudo"), c("Terrapene_mexicana", "Terrapene_carolina", "Terrapene_ornata",
"Terrapene_nelsoni"), "Terrapene_coahuila", c("Glyptemys_muhlenbergii",
"Glyptemys_insculpta"), "Platemys_platycephala")
# Get habitat states for all species
habitat <- turtle_dat$Habitat
names(habitat) <- turtle_dat$Species
age.mat <- as.data.frame(matrix(nrow = length(turtle_tree$tip.label), ncol = 2))
rownames(age.mat) <- turtle_dat$Species
age.mat[, 2] <- "A"
# just to start filling in the habitats
for (i in seq(nodes)) {
foo <- age.regime(x = nodes[i], tr = turtle_tree, name = name.list[[i]],
dat = habitat, check.regimes = FALSE)
age.mat[rownames(foo), ] <- foo
}
age.mat[, 2] <- as.factor(age.mat[, 2])
all(age.mat[, 2] == habitat) #TRUE## [1] TRUElevels(as.factor(age.mat$V1))## [1] "6.970332804" "8.028824792" "8.178711937" "8.905407712" "9.3549541895"
## [6] "12.373138447" "14.088181072" "14.11416946" "14.52109394" "18.73936313"
## [11] "19.15355694" "22.164540147" "36.97324032" "59.902911403"# 14 different time lengths for aquatic lineages that switched to
# terrestrial or secondary transitions from terrestrial back to aquatic
## Now calculating the multivariate distance from optima
turtle_optima <- as.data.frame(ouch_sum$OU2$optima)
turtle_optima## PC1 PC2
## A -0.04126470 0.008458624
## T 0.07602188 -0.013919114# use it like this for arg optima
turtle_data <- turtle_dat[, 2:4]
# Use it
distances <- dist.from.opt(dat = turtle_data, model = "Habitat", optima = turtle_optima,
cont.char = paste0("PC", 1:2))
# function to calculate species' phenotypes distances
head(distances)## 11 78 79 41 13 100
## 0.030147312 0.040827083 0.006758528 0.032684473 0.034097257 0.078763303# Append to previous matrix
age.mat$dist <- distances
age.mat <- age.mat[, c(1, 3, 2)]
colnames(age.mat) <- c("time", "dist", "habitat")
# Remove all aquatic taxa that have a history of 100% water(they have 'NA'
# for time values)
age.mat <- age.mat[complete.cases(age.mat), ]
red.tr <- drop.tip(turtle_tree, name.check(turtle_tree, age.mat)$tree_not_data)
age.mat <- age.mat[red.tr$tip.label, ]
# rownames(age.mat) == red.tr$tip.label #TRUE
## Using BM with phylolm package
head(age.mat)## time dist habitat
## Cuora_mccordi 9.354954 0.04186560 T
## Cuora_picturata 9.354954 0.03213078 T
## Cuora_bourreti 9.354954 0.03249231 T
## Cuora_galbinifrons 9.354954 0.02695018 T
## Cuora_flavomarginata 9.354954 0.04786849 T
## Cuora_mouhotii 6.970333 0.04793729 Tlog.dat <- cbind(time = log(age.mat[, 1]), age.mat[, 2:3])
# Need to use absolute distances
head(log.dat$dist)## [1] 0.04186560 0.03213078 0.03249231 0.02695018 0.04786849 0.04793729# data frame must have row names
phy.log.lm <- summary(phylolm(dist ~ time, log.dat, red.tr, model = "BM"))
phy.log.lm##
## Call:
## phylolm(formula = dist ~ time, data = log.dat, phy = red.tr,
## model = "BM")
##
## AIC logLik
## -331.6 168.8
##
## Raw residuals:
## Min 1Q Median 3Q Max
## -0.08071 -0.06128 -0.04230 -0.01770 0.04783
##
## Mean tip height: 227.7755
## Parameter estimate(s) using ML:
## sigma2: 3.462345e-05
##
## Coefficients:
## Estimate StdErr t.value p.value
## (Intercept) 0.0997085 0.0686145 1.4532 0.1506
## time -0.0019733 0.0107082 -0.1843 0.8543
##
## R-squared: 0.0004781 Adjusted R-squared: -0.0136# Same test excluding Testudinidae
# Excluding Testudinidae: only include rows that have time < 59
age.mat.1 <- age.mat[age.mat$time < 59, ] #n=29
log.dat.1 <- cbind(time = log(age.mat.1[, 1]), age.mat.1[, 2:3])
phy.log.lm.1 <- summary(phylolm(dist ~ time, log.dat.1, red.tr, model = "BM"))## Warning in phylolm(dist ~ time, log.dat.1, red.tr, model = "BM"): will drop from
## the tree 44 taxa with missing dataphy.log.lm.1##
## Call:
## phylolm(formula = dist ~ time, data = log.dat.1, phy = red.tr,
## model = "BM")
##
## AIC logLik
## -145.53 75.76
##
## Raw residuals:
## Min 1Q Median 3Q Max
## -0.07433 -0.05341 -0.03218 -0.01612 0.03370
##
## Mean tip height: 227.7755
## Parameter estimate(s) using ML:
## sigma2: 1.766727e-05
##
## Coefficients:
## Estimate StdErr t.value p.value
## (Intercept) 0.0965134 0.0501587 1.9242 0.06493 .
## time 0.0021319 0.0084593 0.2520 0.80293
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-squared: 0.002347 Adjusted R-squared: -0.0346We can see from the results above that variation in distance from the optimum is not explained by time spent in current environment, even when excluding Testudinidae. These results are shown in Table A2 for multivariate distances.
Figure 3B: PGLS of time spent in current environment on distance from optimum
We can plot these results as shown in Figure 3B of the paper.
# plot in log scale but use original time as labels in axis
colnames(age.mat)[3] <- "Habitat"
Fig.3B <- ggplot(data = age.mat, aes(x = time, y = dist, fill = Habitat, shape = Habitat,
size = Habitat)) + geom_point(color = "black") + scale_size_manual(values = c(4,
3)) + scale_shape_manual(values = c(21, 24)) + scale_x_continuous(trans = "log",
breaks = c(7, 14, 28, 56)) + scale_fill_manual(values = c("gray70", "black")) +
labs(x = "Time in current environment (mya)", y = "Multivariate distance to optimum") +
ylim(0, 0.15) + theme_bw() + theme(text = element_text(size = 15), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.text = element_text(size = 12,
color = "black"), axis.title = element_text(size = 16, color = "black"),
legend.position = "none")
Fig.3B
Conclusion
From all the analyses performed, we can conclude that performance and adaptive peaks match for turtle shells, indicating that the three functions measured (stress, drag, and self-righting ability) are relevant for fitness. Moreover, stress has a higher weight for terrestrial species and the same is true for drag for aquatic species, as one might expect. We can also infer that most species’ phenotypes deviate from inferred adaptive peaks, but probably because of constraints that are not generally related to an influence of past adaptation to ancestral environments.
References
McLaughlin C.J. and C.T. Stayton. 2016. Convergent evolution provides evidence of similar radiations in shell shape in the turtle families Emydidae and Geoemydidae. Herpetologica 72:120–129.
Polly P.D., C.T. Stayton, E.R. Dumont, S.E. Pierce, E.J. Rayfield, and K.D. Angielczyk. 2016. Combining geometric morphometrics and finite element analysis with evolutionary modeling: towards a synthesis. J. Vertebr. Paleontol. 36:e1111225.
Stayton C.T. 2019. Performance in three shell functions predicts the phenotypic distribution of hard‐shelled turtles. Evolution 73:720–734.
Stayton C.T., L.F. O’Connor and N.M. Nisivoccia. 2018. The influence of multiple functional demands on morphological diversification: a test on turtle shells. Evolution 72:1933–1949.