Data analysis for:

Why are telomeres the length that they are? Insight from a phylogenetic comparative analysis

Workflow of the analyses produced by Dylan Padilla, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.


Library 

library(AICcmodavg)
library(ape)
library(caper)
library(car)
library(coda)
library(extrafont)
library(geiger)
library(graph)
library(kableExtra)
library(MuMIn)
library(nlme)
library(pbapply)
library(phylopath)
library(phytools)
library(plotrix)
library(rphylopic)
library(scales)
library(sensiPhy)
library(shape)
library(xtable)
R.version
>                _                           
> platform       aarch64-apple-darwin20      
> arch           aarch64                     
> os             darwin20                    
> system         aarch64, darwin20           
> status                                     
> major          4                           
> minor          3.3                         
> year           2024                        
> month          02                          
> day            29                          
> svn rev        86002                       
> language       R                           
> version.string R version 4.3.3 (2024-02-29)
> nickname       Angel Food Cake
sessionInfo()
> R version 4.3.3 (2024-02-29)
> Platform: aarch64-apple-darwin20 (64-bit)
> Running under: macOS 15.1
> 
> Matrix products: default
> BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
> LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
> 
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
> 
> time zone: America/Phoenix
> tzcode source: internal
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base     
> 
> other attached packages:
>  [1] rmarkdown_2.26      xtable_1.8-4        shape_1.4.6.1       sensiPhy_0.8.5      ggplot2_3.5.1      
>  [6] phylolm_2.6.2       scales_1.3.0        rphylopic_1.4.0     plotrix_3.8-4       phylopath_1.2.1    
> [11] pbapply_1.7-2       nlme_3.1-164        MuMIn_1.47.5        kableExtra_1.4.0    graph_1.80.0       
> [16] BiocGenerics_0.48.1 geiger_2.0.11       phytools_2.1-1      maps_3.4.2          extrafont_0.19     
> [21] coda_0.19-4.1       car_3.1-2           carData_3.0-5       caper_1.0.3         mvtnorm_1.2-4      
> [26] MASS_7.3-60.0.1     ape_5.8             AICcmodavg_2.3-3   
> 
> loaded via a namespace (and not attached):
>  [1] mnormt_2.1.1            gridExtra_2.3           phangorn_2.11.1         rlang_1.1.3            
>  [5] magrittr_2.0.3          compiler_4.3.3          png_0.1-8               systemfonts_1.0.6      
>  [9] vctrs_0.6.5             combinat_0.0-8          quadprog_1.5-8          stringr_1.5.1          
> [13] pkgconfig_2.0.3         fastmap_1.1.1           labeling_0.4.3          ggraph_2.2.1           
> [17] utf8_1.2.4              subplex_1.8             deSolve_1.40            purrr_1.0.2            
> [21] xfun_0.43               cachem_1.0.8            clusterGeneration_1.3.8 jsonlite_1.8.8         
> [25] tweenr_2.0.3            jpeg_0.1-10             VGAM_1.1-10             parallel_4.3.3         
> [29] R6_2.5.1                bslib_0.7.0             stringi_1.8.3           parallelly_1.37.1      
> [33] extrafontdb_1.0         jquerylib_0.1.4         numDeriv_2016.8-1.1     Rcpp_1.0.12            
> [37] iterators_1.0.14        knitr_1.45              future.apply_1.11.2     optimParallel_1.0-2    
> [41] ggm_2.5.1               base64enc_0.1-3         Matrix_1.6-5            splines_4.3.3          
> [45] igraph_2.0.3            tidyselect_1.2.1        yaml_2.3.8              viridis_0.6.5          
> [49] rstudioapi_0.16.0       abind_1.4-5             doParallel_1.0.17       codetools_0.2-19       
> [53] curl_5.2.1              listenv_0.9.1           lattice_0.22-5          tibble_3.2.1           
> [57] withr_3.0.0             evaluate_0.23           future_1.33.2           survival_3.5-8         
> [61] polyclip_1.10-6         xml2_1.3.6              BiocManager_1.30.22     pillar_1.9.0           
> [65] foreach_1.5.2           stats4_4.3.3            generics_0.1.3          grImport2_0.3-1        
> [69] munsell_0.5.1           globals_0.16.3          glue_1.7.0              unmarked_1.4.1         
> [73] scatterplot3d_0.3-44    tools_4.3.3             graphlayouts_1.1.1      XML_3.99-0.16.1        
> [77] tidygraph_1.3.1         fastmatch_1.1-4         grid_4.3.3              tidyr_1.3.1            
> [81] Rttf2pt1_1.3.12         colorspace_2.1-0        ggforce_0.4.2           cli_3.6.2              
> [85] rsvg_2.6.0              fansi_1.0.6             expm_0.999-9            viridisLite_0.4.2      
> [89] svglite_2.1.3           dplyr_1.1.4             gtable_0.3.5            sass_0.4.9             
> [93] digest_0.6.35           ggrepel_0.9.5           farver_2.1.1            memoise_2.0.1          
> [97] htmltools_0.5.8.1       lifecycle_1.0.4         httr_1.4.7
set.seed(80)

## Dataset

data <- read.csv('data.csv')
str(data)
> 'data.frame': 122 obs. of  26 variables:
>  $ Common.name               : chr  "Barn swallow " "Northern goshawks " "Orange-winged Amazon parrot " "Red tail hawk  " ...
>  $ Scientific_name           : chr  "Hirundo_rustica" "Accipiter_gentilis" "Amazona_amazonica" "Buteo_jamaicensis" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Aves" "Aves" "Aves" "Aves" ...
>  $ Order                     : chr  "Passeriformes" "Accipitriformes" "Psittaciformes" "Accipitriformes" ...
>  $ Family                    : chr  "Hirundinidae" "Accipitridae" "Psittacidae" "Accipitridae" ...
>  $ Genus                     : chr  "Hirundo" "Accipiter" "Amazona" "Buteo" ...
>  $ Species                   : chr  "rustica" "gentilis" "amazonica" "jamaicensis" ...
>  $ Endo_ectotherm            : chr  "endo" "endo" "endo" "endo" ...
>  $ Adult_mass_grams          : num  18.3 1044 370 1362 980.6 ...
>  $ Lifespan_years            : num  16 22 30 30.7 24.9 13.4 6.4 15 50 17 ...
>  $ Average_Telomere_Length_kb: num  8.7 1.3 11.8 13.3 42.4 3.2 23.4 3.44 9.9 30.2 ...
>  $ Telomerase_activity       : int  NA NA NA NA NA NA NA NA NA NA ...
>  $ Tissue.type.for.TA        : chr  NA NA NA NA ...
>  $ Source.for.TA             : chr  NA NA NA NA ...
>  $ Tissue.type.for.TL        : chr  "Blood" "Blood" "Blood" "Blood" ...
>  $ Tissue_coded              : chr  "A" "A" "A" "A" ...
>  $ Methodology.for.TL        : chr  "TRF (combination)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
>  $ method                    : chr  "A" "B" "B" "B" ...
>  $ Source.for.TL             : chr  "Caprioli 2012; Foote et al. 2013; Tricola et al. 2018; Criscuolo et al. 2021" "Delany et al. 2000; Criscuolo et al. 2021" "Delany et al. 2000; Criscuolo et al. 2021" "Delany et al. 2000" ...
>  $ Source.for.Lifespan       : chr  "AnAge" "AnAge" "AnAge" "AnAge" ...
>  $ Source.for.Mass           : chr  "AnAge" "AnAge" "AnAge" "AnAge" ...
>  $ Captivity                 : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Notes                     : chr  "" "" "" "" ...
head(data)
>                    Common.name     Scientific_name    Domain Kingdom   Phylum Class           Order
> 1                Barn swallow      Hirundo_rustica Eukaryota Metazoa Chordata  Aves   Passeriformes
> 2           Northern goshawks   Accipiter_gentilis Eukaryota Metazoa Chordata  Aves Accipitriformes
> 3 Orange-winged Amazon parrot    Amazona_amazonica Eukaryota Metazoa Chordata  Aves  Psittaciformes
> 4              Red tail hawk     Buteo_jamaicensis Eukaryota Metazoa Chordata  Aves Accipitriformes
> 5             Swainson's hawk      Buteo_swainsoni Eukaryota Metazoa Chordata  Aves Accipitriformes
> 6             Golden pheasant  Chrysolophus_pictus Eukaryota Metazoa Chordata  Aves     Galliformes
>         Family        Genus     Species Endo_ectotherm Adult_mass_grams Lifespan_years
> 1 Hirundinidae      Hirundo     rustica           endo            18.30           16.0
> 2 Accipitridae    Accipiter    gentilis           endo          1044.00           22.0
> 3  Psittacidae      Amazona   amazonica           endo           370.00           30.0
> 4 Accipitridae        Buteo jamaicensis           endo          1362.00           30.7
> 5 Accipitridae        Buteo   swainsoni           endo           980.64           24.9
> 6  Phasianidae Chrysolophus      pictus           endo           625.00           13.4
>   Average_Telomere_Length_kb Telomerase_activity Tissue.type.for.TA Source.for.TA Tissue.type.for.TL
> 1                        8.7                  NA               <NA>          <NA>              Blood
> 2                        1.3                  NA               <NA>          <NA>              Blood
> 3                       11.8                  NA               <NA>          <NA>              Blood
> 4                       13.3                  NA               <NA>          <NA>              Blood
> 5                       42.4                  NA               <NA>          <NA>              Blood
> 6                        3.2                  NA               <NA>          <NA>              Blood
>   Tissue_coded Methodology.for.TL method
> 1            A  TRF (combination)      A
> 2            A    TRF (denatured)      B
> 3            A    TRF (denatured)      B
> 4            A    TRF (denatured)      B
> 5            A    TRF (denatured)      B
> 6            A    TRF (denatured)      B
>                                                                  Source.for.TL Source.for.Lifespan
> 1 Caprioli 2012; Foote et al. 2013; Tricola et al. 2018; Criscuolo et al. 2021               AnAge
> 2                                    Delany et al. 2000; Criscuolo et al. 2021               AnAge
> 3                                    Delany et al. 2000; Criscuolo et al. 2021               AnAge
> 4                                                           Delany et al. 2000               AnAge
> 5                                    Delany et al. 2000; Criscuolo et al. 2021               AnAge
> 6                                   Delaney et al. 2000; Criscuolo et al. 2021               AnAge
>   Source.for.Mass Captivity Notes
> 1           AnAge         0      
> 2           AnAge         0      
> 3           AnAge         0      
> 4           AnAge         0      
> 5           AnAge         0      
> 6           AnAge         0
unique(data$Family)
>  [1] "Hirundinidae"      "Accipitridae"      "Psittacidae"       "Phasianidae"       "Odontophoridae"   
>  [6] "Moronidae"         "Diomedeidae"       "Falconidae"        "Procellaridae"     "Laridae"          
> [11] "Strigopidae"       "Apodidae"          "Carangidae"        "Cyprinidae"        "Lutjanidae"       
> [16] "Mullidae"          "Sparidae"          "Percichthyidae"    "Nothobranchiidae"  "Gadidae"          
> [21] "Platycephalidae"   "Salmonidae"        "Neoceratodontidae" "Urolophidae"       "Callorhinchidae"  
> [26] "Muridae"           "Chinchillidae"     "Cricetidae"        "Myocastoridae"     "Sciuridae"        
> [31] "Otariidae"         "Castoridae"        "Atelidae"          "Hyaenidae"         "Aplodontiidae"    
> [36] "Balaenidae"        "Hominidae"         "Pteropodidae"      "Cebidae"           "Rhinocerotidae"   
> [41] "Giraffidae"        "Caviidae"          "Equidae"           "Myrmecophagidae"   "Tapiridae"        
> [46] "Ursidae"           "Chlamyphoridae"    "Eschrichtiidae"    "Elephantidae"      "Cervidae"         
> [51] "Procaviidae"       "Tupaiidae"         "Heterocephalidae"  "Cercopithecidae"   "Choloepodidae"    
> [56] "Delphinidae"       "Lemuridae"         "Ochotonidae"       "Ailuridae"         "Molossidae"       
> [61] "Vespertilionidae"  "Didelphidae"       "Macroscelididae"   "Erinaceidae"       "Leporidae"        
> [66] "Felidae"           "Tenrecidae"        "Alligatoridae"     "Lacertidae"        "Pythonidae"       
> [71] "Emydidae"          "Estrildidae"       "Corvidae"          "Anatidae"          "Scolopacidae"     
> [76] "Alcidae"           "Fregatidae"        "Haematopotidae"    "Hydrobatidae"      "Paridae"          
> [81] "Emberizidae"       "Spheniscidae"      "Sternidae"         "Pelecani"          "Turdidae"         
> [86] "Vireonidae"        "Poeciliidae"       "Testudinoidea"
dat <- data
names(dat)
>  [1] "Common.name"                "Scientific_name"            "Domain"                    
>  [4] "Kingdom"                    "Phylum"                     "Class"                     
>  [7] "Order"                      "Family"                     "Genus"                     
> [10] "Species"                    "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "Lifespan_years"             "Average_Telomere_Length_kb" "Telomerase_activity"       
> [16] "Tissue.type.for.TA"         "Source.for.TA"              "Tissue.type.for.TL"        
> [19] "Tissue_coded"               "Methodology.for.TL"         "method"                    
> [22] "Source.for.TL"              "Source.for.Lifespan"        "Source.for.Mass"           
> [25] "Captivity"                  "Notes"
unique(is.na(dat$log_mass))
> logical(0)
## Trees

##write.table(data$Scientific_name, 'spp.tree.08-03-24.txt', row.names = FALSE, col.names = FALSE)

full_data_tree <- read.tree("spp.tree.08-03-24.nwk")
is.ultrametric(full_data_tree)
> [1] FALSE
full_data_tree <- force.ultrametric(full_data_tree)
> ***************************************************************
> *                          Note:                              *
> *    force.ultrametric does not include a formal method to    *
> *    ultrametricize a tree & should only be used to coerce    *
> *   a phylogeny that fails is.ultrametric due to rounding --  *
> *    not as a substitute for formal rate-smoothing methods.   *
> ***************************************************************
is.ultrametric(full_data_tree)
> [1] TRUE
full_data_tree
> 
> Phylogenetic tree with 138 tips and 137 internal nodes.
> 
> Tip labels:
>   Urolophus_paucimaculatus, Callorhinchus_milii, Oncorhynchus_mykiss, Lutjanus_argentimaculatus, Acanthopagrus_schlegelii, Dicentrarchus_labrax, ...
> Node labels:
>   , '14', '12', '48', '57', '49', ...
> 
> Rooted; includes branch lengths.
rownames(dat) <- dat$Scientific_name


## Full_tree

check <- name.check(full_data_tree, dat)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
pruned_data_tree <- drop.tip(full_data_tree, rm_phy)
pruned_dat <- subset(dat, subset = dat$Scientific_name %in% pruned_data_tree$tip, select = names(dat))
str(pruned_dat)
> 'data.frame': 122 obs. of  26 variables:
>  $ Common.name               : chr  "Barn swallow " "Northern goshawks " "Orange-winged Amazon parrot " "Red tail hawk  " ...
>  $ Scientific_name           : chr  "Hirundo_rustica" "Accipiter_gentilis" "Amazona_amazonica" "Buteo_jamaicensis" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Aves" "Aves" "Aves" "Aves" ...
>  $ Order                     : chr  "Passeriformes" "Accipitriformes" "Psittaciformes" "Accipitriformes" ...
>  $ Family                    : chr  "Hirundinidae" "Accipitridae" "Psittacidae" "Accipitridae" ...
>  $ Genus                     : chr  "Hirundo" "Accipiter" "Amazona" "Buteo" ...
>  $ Species                   : chr  "rustica" "gentilis" "amazonica" "jamaicensis" ...
>  $ Endo_ectotherm            : chr  "endo" "endo" "endo" "endo" ...
>  $ Adult_mass_grams          : num  18.3 1044 370 1362 980.6 ...
>  $ Lifespan_years            : num  16 22 30 30.7 24.9 13.4 6.4 15 50 17 ...
>  $ Average_Telomere_Length_kb: num  8.7 1.3 11.8 13.3 42.4 3.2 23.4 3.44 9.9 30.2 ...
>  $ Telomerase_activity       : int  NA NA NA NA NA NA NA NA NA NA ...
>  $ Tissue.type.for.TA        : chr  NA NA NA NA ...
>  $ Source.for.TA             : chr  NA NA NA NA ...
>  $ Tissue.type.for.TL        : chr  "Blood" "Blood" "Blood" "Blood" ...
>  $ Tissue_coded              : chr  "A" "A" "A" "A" ...
>  $ Methodology.for.TL        : chr  "TRF (combination)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
>  $ method                    : chr  "A" "B" "B" "B" ...
>  $ Source.for.TL             : chr  "Caprioli 2012; Foote et al. 2013; Tricola et al. 2018; Criscuolo et al. 2021" "Delany et al. 2000; Criscuolo et al. 2021" "Delany et al. 2000; Criscuolo et al. 2021" "Delany et al. 2000" ...
>  $ Source.for.Lifespan       : chr  "AnAge" "AnAge" "AnAge" "AnAge" ...
>  $ Source.for.Mass           : chr  "AnAge" "AnAge" "AnAge" "AnAge" ...
>  $ Captivity                 : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Notes                     : chr  "" "" "" "" ...
head(pruned_dat)
>                                      Common.name     Scientific_name    Domain Kingdom   Phylum Class
> Hirundo_rustica                    Barn swallow      Hirundo_rustica Eukaryota Metazoa Chordata  Aves
> Accipiter_gentilis            Northern goshawks   Accipiter_gentilis Eukaryota Metazoa Chordata  Aves
> Amazona_amazonica   Orange-winged Amazon parrot    Amazona_amazonica Eukaryota Metazoa Chordata  Aves
> Buteo_jamaicensis                Red tail hawk     Buteo_jamaicensis Eukaryota Metazoa Chordata  Aves
> Buteo_swainsoni                 Swainson's hawk      Buteo_swainsoni Eukaryota Metazoa Chordata  Aves
> Chrysolophus_pictus             Golden pheasant  Chrysolophus_pictus Eukaryota Metazoa Chordata  Aves
>                               Order       Family        Genus     Species Endo_ectotherm
> Hirundo_rustica       Passeriformes Hirundinidae      Hirundo     rustica           endo
> Accipiter_gentilis  Accipitriformes Accipitridae    Accipiter    gentilis           endo
> Amazona_amazonica    Psittaciformes  Psittacidae      Amazona   amazonica           endo
> Buteo_jamaicensis   Accipitriformes Accipitridae        Buteo jamaicensis           endo
> Buteo_swainsoni     Accipitriformes Accipitridae        Buteo   swainsoni           endo
> Chrysolophus_pictus     Galliformes  Phasianidae Chrysolophus      pictus           endo
>                     Adult_mass_grams Lifespan_years Average_Telomere_Length_kb Telomerase_activity
> Hirundo_rustica                18.30           16.0                        8.7                  NA
> Accipiter_gentilis           1044.00           22.0                        1.3                  NA
> Amazona_amazonica             370.00           30.0                       11.8                  NA
> Buteo_jamaicensis            1362.00           30.7                       13.3                  NA
> Buteo_swainsoni               980.64           24.9                       42.4                  NA
> Chrysolophus_pictus           625.00           13.4                        3.2                  NA
>                     Tissue.type.for.TA Source.for.TA Tissue.type.for.TL Tissue_coded Methodology.for.TL
> Hirundo_rustica                   <NA>          <NA>              Blood            A  TRF (combination)
> Accipiter_gentilis                <NA>          <NA>              Blood            A    TRF (denatured)
> Amazona_amazonica                 <NA>          <NA>              Blood            A    TRF (denatured)
> Buteo_jamaicensis                 <NA>          <NA>              Blood            A    TRF (denatured)
> Buteo_swainsoni                   <NA>          <NA>              Blood            A    TRF (denatured)
> Chrysolophus_pictus               <NA>          <NA>              Blood            A    TRF (denatured)
>                     method                                                                Source.for.TL
> Hirundo_rustica          A Caprioli 2012; Foote et al. 2013; Tricola et al. 2018; Criscuolo et al. 2021
> Accipiter_gentilis       B                                    Delany et al. 2000; Criscuolo et al. 2021
> Amazona_amazonica        B                                    Delany et al. 2000; Criscuolo et al. 2021
> Buteo_jamaicensis        B                                                           Delany et al. 2000
> Buteo_swainsoni          B                                    Delany et al. 2000; Criscuolo et al. 2021
> Chrysolophus_pictus      B                                   Delaney et al. 2000; Criscuolo et al. 2021
>                     Source.for.Lifespan Source.for.Mass Captivity Notes
> Hirundo_rustica                   AnAge           AnAge         0      
> Accipiter_gentilis                AnAge           AnAge         0      
> Amazona_amazonica                 AnAge           AnAge         0      
> Buteo_jamaicensis                 AnAge           AnAge         0      
> Buteo_swainsoni                   AnAge           AnAge         0      
> Chrysolophus_pictus               AnAge           AnAge         0
pruned_data_tree
> 
> Phylogenetic tree with 122 tips and 121 internal nodes.
> 
> Tip labels:
>   Urolophus_paucimaculatus, Callorhinchus_milii, Oncorhynchus_mykiss, Lutjanus_argentimaculatus, Acanthopagrus_schlegelii, Dicentrarchus_labrax, ...
> Node labels:
>   , '14', '12', '48', '57', '49', ...
> 
> Rooted; includes branch lengths.
name.check(pruned_data_tree, pruned_dat)
> [1] "OK"
## Checking the distribution of lifespan

plot(NA, xlim = c(0, 200), ylim = c(0, 60), type = 'n', las = 1, xlab = '', ylab = '', axes = FALSE)
grid()

par(new = TRUE)
hist(pruned_dat$Lifespan_years, main = "Raw variable", las = 1, xlab = 'Lifespan (years)')
box()

plot(NA, xlim = c(0, 200), ylim = c(0, 60), type = 'n', las = 1, xlab = '', ylab = '', axes = FALSE)
grid()

par(new = TRUE)
hist(log1p(pruned_dat$Lifespan_years), main = "Log-transformed", xlab = 'Log lifespan (years)', las = 1)
box()

## Log-transforming for a better fit

pruned_dat$log.lifespan <- log1p(pruned_dat$Lifespan_years)
pruned_dat$log.mass <- log1p(pruned_dat$Adult_mass_grams)
## Defining models


str(pruned_dat)
> 'data.frame': 122 obs. of  28 variables:
>  $ Common.name               : chr  "Barn swallow " "Northern goshawks " "Orange-winged Amazon parrot " "Red tail hawk  " ...
>  $ Scientific_name           : chr  "Hirundo_rustica" "Accipiter_gentilis" "Amazona_amazonica" "Buteo_jamaicensis" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Aves" "Aves" "Aves" "Aves" ...
>  $ Order                     : chr  "Passeriformes" "Accipitriformes" "Psittaciformes" "Accipitriformes" ...
>  $ Family                    : chr  "Hirundinidae" "Accipitridae" "Psittacidae" "Accipitridae" ...
>  $ Genus                     : chr  "Hirundo" "Accipiter" "Amazona" "Buteo" ...
>  $ Species                   : chr  "rustica" "gentilis" "amazonica" "jamaicensis" ...
>  $ Endo_ectotherm            : chr  "endo" "endo" "endo" "endo" ...
>  $ Adult_mass_grams          : num  18.3 1044 370 1362 980.6 ...
>  $ Lifespan_years            : num  16 22 30 30.7 24.9 13.4 6.4 15 50 17 ...
>  $ Average_Telomere_Length_kb: num  8.7 1.3 11.8 13.3 42.4 3.2 23.4 3.44 9.9 30.2 ...
>  $ Telomerase_activity       : int  NA NA NA NA NA NA NA NA NA NA ...
>  $ Tissue.type.for.TA        : chr  NA NA NA NA ...
>  $ Source.for.TA             : chr  NA NA NA NA ...
>  $ Tissue.type.for.TL        : chr  "Blood" "Blood" "Blood" "Blood" ...
>  $ Tissue_coded              : chr  "A" "A" "A" "A" ...
>  $ Methodology.for.TL        : chr  "TRF (combination)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
>  $ method                    : chr  "A" "B" "B" "B" ...
>  $ Source.for.TL             : chr  "Caprioli 2012; Foote et al. 2013; Tricola et al. 2018; Criscuolo et al. 2021" "Delany et al. 2000; Criscuolo et al. 2021" "Delany et al. 2000; Criscuolo et al. 2021" "Delany et al. 2000" ...
>  $ Source.for.Lifespan       : chr  "AnAge" "AnAge" "AnAge" "AnAge" ...
>  $ Source.for.Mass           : chr  "AnAge" "AnAge" "AnAge" "AnAge" ...
>  $ Captivity                 : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Notes                     : chr  "" "" "" "" ...
>  $ log.lifespan              : num  2.83 3.14 3.43 3.46 3.25 ...
>  $ log.mass                  : num  2.96 6.95 5.92 7.22 6.89 ...
path.mod <- pruned_dat
names(path.mod)
>  [1] "Common.name"                "Scientific_name"            "Domain"                    
>  [4] "Kingdom"                    "Phylum"                     "Class"                     
>  [7] "Order"                      "Family"                     "Genus"                     
> [10] "Species"                    "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "Lifespan_years"             "Average_Telomere_Length_kb" "Telomerase_activity"       
> [16] "Tissue.type.for.TA"         "Source.for.TA"              "Tissue.type.for.TL"        
> [19] "Tissue_coded"               "Methodology.for.TL"         "method"                    
> [22] "Source.for.TL"              "Source.for.Lifespan"        "Source.for.Mass"           
> [25] "Captivity"                  "Notes"                      "log.lifespan"              
> [28] "log.mass"
path.mod$Average_Telomere_Length_kb <- log(path.mod$Average_Telomere_Length_kb)

names(path.mod)[c(14, 11, 28, 27)] <- c('TL', 'M', 'BM', 'LS')

models <- define_model_set(
    "null" = c(TL ~ TL, M ~ LS, M ~ BM, LS ~ BM),
    "(a)" = c(TL ~ LS, M ~ LS, M ~ BM, BM ~ LS),
    "(b)" = c(TL ~ BM, M ~ LS, M ~ BM, BM ~ LS),
    "(c)" = c(TL ~ M, M ~ LS, M ~ BM, BM ~ LS),
    "(d)" = c(TL ~ BM + LS + M),
    "(e)" = c(TL ~ LS),
    "(f)" = c(TL ~ BM),
    "(g)" = c(TL ~ M))


#png("figureS1-revised.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS1-revised.pdf")

plot_model_set(models, nrow = 4) + ggplot2::scale_y_continuous(expand = c(0.1, 0)) + theme(plot.margin=unit(c(0.5, 3, 0.5, 3), "cm"))

#dev.off()

mo <- phylo_path(models, path.mod, pruned_data_tree)
s <- summary(mo)
s
>      model k q     C        p CICc delta_CICc        l        w
> (b)    (b) 2 8  4.90 2.98e-01 22.2       0.00 1.00e+00 6.42e-01
> (a)    (a) 2 8  6.14 1.89e-01 23.4       1.25 5.37e-01 3.44e-01
> null  null 3 7 16.13 1.31e-02 31.1       8.94 1.14e-02 7.35e-03
> (c)    (c) 2 8 14.08 7.06e-03 31.3       9.18 1.02e-02 6.53e-03
> (d)    (d) 3 7 48.20 1.08e-08 63.2      41.01 1.25e-09 7.99e-10
> (f)    (f) 5 5 52.90 7.77e-08 63.4      41.25 1.11e-09 7.09e-10
> (e)    (e) 5 5 54.48 3.95e-08 65.0      42.82 5.02e-10 3.22e-10
> (g)    (g) 5 5 70.06 4.31e-11 80.6      58.41 2.08e-13 1.33e-13
theme_set(theme_bw())


## Conditional independence tested

mo$d_sep$'(b)'
> # A tibble: 2 × 4
>   d_sep                p phylo_par model    
>   <chr>            <dbl>     <dbl> <list>   
> 1 TL ~ BM + LS     0.302     0.538 <phylolm>
> 2 TL ~ LS + BM + M 0.286     0.529 <phylolm>
bar.pl <- tapply(s$w, list(s$model), print)
> [1] 0.344336
> [1] 0.6417914
> [1] 0.006525674
> [1] 7.99095e-10
> [1] 3.224665e-10
> [1] 7.094764e-10
> [1] 1.332557e-13
> [1] 0.007346913
p.val <- s[names(bar.pl), ]
knitr::kable(s, digits = 3, caption = 'Competing models of the evolution of telomere length across vertebrates', align = 'l', row.names = FALSE, table.attr = "style='width:30%;'")
Competing models of the evolution of telomere length across vertebrates
model k q C p CICc delta_CICc l w
(b) 2 8 4.898 0.298 22.172 0.000 1.000 0.642
(a) 2 8 6.143 0.189 23.418 1.245 0.537 0.344
null 3 7 16.130 0.013 31.112 8.940 0.011 0.007
(c) 2 8 14.075 0.007 31.349 9.177 0.010 0.007
(d) 3 7 48.198 0.000 63.180 41.008 0.000 0.000
(f) 5 5 52.901 0.000 63.418 41.246 0.000 0.000
(e) 5 5 54.478 0.000 64.995 42.823 0.000 0.000
(g) 5 5 70.061 0.000 80.578 58.406 0.000 0.000 
bt <- best(mo)
plot(bt, algorithm = 'mds') + theme(legend.position= 'none')
A visualization of the best supported causal model, and the standardized path coefficients.

A visualization of the best supported causal model, and the standardized path coefficients. 

##png("figureS1.png", width = 7, height = 7, units = "in", res = 360)


barplot(NA, beside = TRUE, space = FALSE, horiz = TRUE, xlim = c(0, 1), axes = FALSE)

grid()

par(new = TRUE)

## Adjust the col argument to make the "b" bar red

barplot(bar.pl, beside = TRUE, space = FALSE, horiz = TRUE, xlim = c(0, 1), yaxt = 'n', ylab = 'Candidate models', xlab='Weight of evidence (w)', col = c('gray', alpha('red', 0.3), rep('gray', 6)))

axis(2, at = (0:7)+0.5, labels = names(bar.pl), las = 1)
box()

text(x = round(bar.pl, 3)+0.05, y = (0:7)+0.5, labels = round(p.val$p, 3))
legend('topright', legend = 'within 2 CICc', fill = alpha('red', 0.3), bty = 'n')
Relative importance of the causal models describing the evolution of telomere length across vertebrates. Bar labels are p-values; significance indicates rejection.

Relative importance of the causal models describing the evolution of telomere length across vertebrates. Bar labels are p-values; significance indicates rejection. 

##dev.off()
knitr::kable(mo$d_sep$'(b)'[1:2, 1:3], digits = 3, caption = 'Conditional independence statements tested', align = 'l', row.names = FALSE, table.attr = "style='width:30%;'")
Conditional independence statements tested
d_sep p phylo_par
TL ~ BM + LS 0.302 0.538
TL ~ LS + BM + M 0.286 0.529



Fitting models of evolution based on the conditional independence statements tested in the most likely path model 



pagel.model0 <- gls(log(Average_Telomere_Length_kb) ~ log.mass + log.lifespan + Endo_ectotherm, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model0)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ log.mass + log.lifespan + Endo_ectotherm 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   239.4886 256.3128 -113.7443
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.5294452 
> 
> Coefficients:
>                        Value Std.Error   t-value p-value
> (Intercept)         3.524075 0.3908061  9.017450  0.0000
> log.mass           -0.036372 0.0265466 -1.370136  0.1732
> log.lifespan       -0.130401 0.1185554 -1.099913  0.2736
> Endo_ectothermendo -0.396064 0.3697523 -1.071160  0.2863
> 
>  Correlation: 
>                    (Intr) lg.mss lg.lfs
> log.mass            0.125              
> log.lifespan       -0.557 -0.698       
> Endo_ectothermendo -0.276 -0.033  0.062
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -2.8434026 -0.3795892  0.2053350  0.7166696  2.1165328 
> 
> Residual standard error: 0.7751007 
> Degrees of freedom: 122 total; 118 residual
pagel.model1 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.lifespan + log.mass, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model1)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.lifespan +      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   236.9611 256.5893 -111.4806
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.3972721 
> 
> Coefficients:
>                                      Value Std.Error   t-value p-value
> (Intercept)                      2.5930652 0.5014725  5.170902  0.0000
> Endo_ectothermendo               0.9749552 0.6485376  1.503313  0.1355
> log.lifespan                     0.1503982 0.1632280  0.921400  0.3587
> log.mass                        -0.0325013 0.0257028 -1.264504  0.2086
> Endo_ectothermendo:log.lifespan -0.4017335 0.1745440 -2.301617  0.0231
> 
>  Correlation: 
>                                 (Intr) End_ct lg.lfs lg.mss
> Endo_ectothermendo              -0.728                     
> log.lifespan                    -0.802  0.646              
> log.mass                         0.070  0.024 -0.461       
> Endo_ectothermendo:log.lifespan  0.717 -0.879 -0.717 -0.041
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -3.2633469 -0.4613341  0.1622935  0.7012716  2.2145824 
> 
> Residual standard error: 0.7022405 
> Degrees of freedom: 122 total; 117 residual
pagel.model2 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass + log.lifespan, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model2)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass +      log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   235.6954 255.3236 -110.8477
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.4719298 
> 
> Coefficients:
>                                  Value Std.Error   t-value p-value
> (Intercept)                  2.9464763 0.4285141  6.876032  0.0000
> Endo_ectothermendo           0.3990165 0.4646744  0.858701  0.3923
> log.mass                     0.0491984 0.0435613  1.129404  0.2610
> log.lifespan                -0.1561544 0.1155990 -1.350828  0.1794
> Endo_ectothermendo:log.mass -0.1052738 0.0435059 -2.419761  0.0171
> 
>  Correlation: 
>                             (Intr) End_ct lg.mss lg.lfs
> Endo_ectothermendo          -0.531                     
> log.mass                    -0.348  0.540              
> log.lifespan                -0.431 -0.033 -0.500       
> Endo_ectothermendo:log.mass  0.516 -0.687 -0.805  0.108
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -3.0213728 -0.4001941  0.2822635  0.7533736  2.3837050 
> 
> Residual standard error: 0.729365 
> Degrees of freedom: 122 total; 117 residual
pagel.model3 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.mass * log.lifespan, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model3)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.mass *      log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   240.8611 260.4892 -113.4305
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.5130464 
> 
> Coefficients:
>                           Value Std.Error   t-value p-value
> (Intercept)            3.227847 0.5298438  6.092073  0.0000
> Endo_ectothermendo    -0.401744 0.3632469 -1.105980  0.2710
> log.mass               0.007697 0.0618546  0.124444  0.9012
> log.lifespan          -0.037529 0.1662055 -0.225801  0.8218
> log.mass:log.lifespan -0.013075 0.0166540 -0.785086  0.4340
> 
>  Correlation: 
>                       (Intr) End_ct lg.mss lg.lfs
> Endo_ectothermendo    -0.174                     
> log.mass              -0.578 -0.051              
> log.lifespan          -0.772  0.014  0.421       
> log.mass:log.lifespan  0.684  0.041 -0.903 -0.702
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -2.8960657 -0.3762918  0.2127450  0.7700036  2.1258145 
> 
> Residual standard error: 0.7646861 
> Degrees of freedom: 122 total; 117 residual
pagel.model4 <- gls(log(Average_Telomere_Length_kb) ~ log.lifespan * Endo_ectotherm * log.mass, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model4)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ log.lifespan * Endo_ectotherm *      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   240.9733 269.0135 -110.4866
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.4414292 
> 
> Coefficients:
>                                               Value Std.Error   t-value p-value
> (Intercept)                               2.9544921 0.7343734  4.023147  0.0001
> log.lifespan                             -0.1493016 0.3398880 -0.439267  0.6613
> Endo_ectothermendo                        0.7884744 1.0325390  0.763627  0.4467
> log.mass                                 -0.0117205 0.1157297 -0.101275  0.9195
> log.lifespan:Endo_ectothermendo          -0.1270265 0.4050161 -0.313633  0.7544
> log.lifespan:log.mass                     0.0169346 0.0417340  0.405775  0.6857
> Endo_ectothermendo:log.mass              -0.0752221 0.1429326 -0.526276  0.5997
> log.lifespan:Endo_ectothermendo:log.mass -0.0065637 0.0472436 -0.138934  0.8897
> 
>  Correlation: 
>                                          (Intr) lg.lfs End_ct lg.mss lg.:E_ lg.l:. End_:.
> log.lifespan                             -0.850                                          
> Endo_ectothermendo                       -0.682  0.634                                   
> log.mass                                 -0.563  0.430  0.361                            
> log.lifespan:Endo_ectothermendo           0.714 -0.838 -0.869 -0.363                     
> log.lifespan:log.mass                     0.704 -0.761 -0.483 -0.873  0.640              
> Endo_ectothermendo:log.mass               0.456 -0.348 -0.591 -0.809  0.484  0.707       
> log.lifespan:Endo_ectothermendo:log.mass -0.622  0.672  0.694  0.772 -0.766 -0.884 -0.883
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -3.1311098 -0.4398032  0.2639437  0.6540715  2.3785288 
> 
> Residual standard error: 0.714069 
> Degrees of freedom: 122 total; 114 residual
pagel.model5 <- gls(log(Average_Telomere_Length_kb) ~ log.lifespan + log.mass, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model5)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ log.lifespan + log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   238.6662 252.6863 -114.3331
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.5375615 
> 
> Coefficients:
>                  Value Std.Error   t-value p-value
> (Intercept)   3.414242 0.3783179  9.024795  0.0000
> log.lifespan -0.123003 0.1185862 -1.037249  0.3017
> log.mass     -0.037411 0.0265794 -1.407516  0.1619
> 
>  Correlation: 
>              (Intr) lg.lfs
> log.lifespan -0.560       
> log.mass      0.119 -0.697
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -3.1999393 -0.6749240 -0.1149099  0.4002848  1.7134961 
> 
> Residual standard error: 0.7831801 
> Degrees of freedom: 122 total; 119 residual
pagel.model6 <- gls(log(Average_Telomere_Length_kb) ~ log.lifespan * log.mass, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model6)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ log.lifespan * log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   240.1265 256.9507 -114.0633
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.5217597 
> 
> Coefficients:
>                            Value Std.Error   t-value p-value
> (Intercept)            3.1343979 0.5245486  5.975419  0.0000
> log.lifespan          -0.0363369 0.1666007 -0.218108  0.8277
> log.mass               0.0036685 0.0619088  0.059256  0.9528
> log.lifespan:log.mass -0.0121864 0.0166738 -0.730868  0.4663
> 
>  Correlation: 
>                       (Intr) lg.lfs lg.mss
> log.lifespan          -0.779              
> log.mass              -0.595  0.422       
> log.lifespan:log.mass  0.701 -0.703 -0.903
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -3.2569026 -0.6993456 -0.1267147  0.4314327  1.7152585 
> 
> Residual standard error: 0.7731187 
> Degrees of freedom: 122 total; 118 residual
output.models <- model.sel(pagel.model0, pagel.model1, pagel.model2, pagel.model3, pagel.model4, pagel.model5, pagel.model6)

sel.table <- round(as.data.frame(output.models)[9:13], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]


knitr::kable(sel.table, digits = 3, caption = 'Model selection procedure based on AIC criterion', row.names = FALSE, align = 'l')
Model selection procedure based on AIC criterion
Model K logLik AICc delta weight
pagel.model2 7 -110.848 236.678 0.000 0.448
pagel.model1 7 -111.481 237.944 1.266 0.238
pagel.model5 5 -114.333 239.183 2.506 0.128
pagel.model0 6 -113.744 240.219 3.541 0.076
pagel.model6 6 -114.063 240.857 4.179 0.055
pagel.model3 7 -113.431 241.844 5.166 0.034
pagel.model4 10 -110.487 242.955 6.277 0.019 
## Estimating lambda from PGLS models

pagel.model0.8 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass + log.lifespan, correlation = corPagel(0.8, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model0.8)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass +      log.lifespan 
>   Data: pruned_dat 
>       AIC      BIC   logLik
>   240.666 257.4901 -114.333
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> lambda 
>    0.8 
> 
> Coefficients:
>                                 Value Std.Error   t-value p-value
> (Intercept)                  3.201521 0.5778259  5.540632  0.0000
> Endo_ectothermendo           0.237381 0.6632568  0.357902  0.7211
> log.mass                     0.036292 0.0501916  0.723071  0.4711
> log.lifespan                -0.179061 0.1275018 -1.404384  0.1629
> Endo_ectothermendo:log.mass -0.093997 0.0506152 -1.857091  0.0658
> 
>  Correlation: 
>                             (Intr) End_ct lg.mss lg.lfs
> Endo_ectothermendo          -0.448                     
> log.mass                    -0.357  0.459              
> log.lifespan                -0.354 -0.007 -0.467       
> Endo_ectothermendo:log.mass  0.475 -0.573 -0.830  0.117
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -2.2119959 -0.3002200  0.1898719  0.5320637  1.7384591 
> 
> Residual standard error: 1.000879 
> Degrees of freedom: 122 total; 117 residual
pagel.model0 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass + log.lifespan, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model0)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass +      log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   260.5828 277.4069 -124.2914
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> lambda 
>      0 
> 
> Coefficients:
>                                  Value Std.Error   t-value p-value
> (Intercept)                  2.2828151 0.3338815  6.837201  0.0000
> Endo_ectothermendo           0.9867956 0.3490592  2.827015  0.0055
> log.mass                     0.0789860 0.0449376  1.757682  0.0814
> log.lifespan                -0.1441376 0.1110251 -1.298243  0.1968
> Endo_ectothermendo:log.mass -0.0989801 0.0444715 -2.225698  0.0280
> 
>  Correlation: 
>                             (Intr) End_ct lg.mss lg.lfs
> Endo_ectothermendo          -0.617                     
> log.mass                    -0.428  0.771              
> log.lifespan                -0.477 -0.256 -0.507       
> Endo_ectothermendo:log.mass  0.586 -0.893 -0.847  0.191
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -3.60527368 -0.60608579 -0.01119332  0.60835477  2.26235421 
> 
> Residual standard error: 0.6702151 
> Degrees of freedom: 122 total; 117 residual
intervals(pagel.model0, which = 'var-cov')
> Approximate 95% confidence intervals
> 
>  Residual standard error:
>     lower      est.     upper 
> 0.5956358 0.6702151 0.7663126
pagel.model0.5 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass + log.lifespan, correlation = corPagel(0.5, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model0.5)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass +      log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC   logLik
>   233.7381 250.5622 -110.869
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> lambda 
>    0.5 
> 
> Coefficients:
>                                  Value Std.Error   t-value p-value
> (Intercept)                  2.9683261 0.4366935  6.797275  0.0000
> Endo_ectothermendo           0.3837552 0.4755101  0.807039  0.4213
> log.mass                     0.0484022 0.0438711  1.103282  0.2722
> log.lifespan                -0.1582973 0.1162848 -1.361289  0.1760
> Endo_ectothermendo:log.mass -0.1045713 0.0438371 -2.385452  0.0187
> 
>  Correlation: 
>                             (Intr) End_ct lg.mss lg.lfs
> Endo_ectothermendo          -0.524                     
> log.mass                    -0.348  0.533              
> log.lifespan                -0.426 -0.030 -0.499       
> Endo_ectothermendo:log.mass  0.513 -0.678 -0.807  0.109
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -2.9670009 -0.3917590  0.2693331  0.7427991  2.3437865 
> 
> Residual standard error: 0.7424674 
> Degrees of freedom: 122 total; 117 residual
intervals(pagel.model0.5, which = 'var-cov')
> Approximate 95% confidence intervals
> 
>  Residual standard error:
>     lower      est.     upper 
> 0.6598481 0.7424674 0.8489247
pagel.model1 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass + log.lifespan, correlation = corPagel(1, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model1)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass +      log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   282.7728 299.5969 -135.3864
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> lambda 
>      1 
> 
> Coefficients:
>                                 Value Std.Error   t-value p-value
> (Intercept)                  3.477902 1.0698466  3.250842  0.0015
> Endo_ectothermendo           0.083476 1.3090875  0.063767  0.9493
> log.mass                     0.019098 0.0736050  0.259463  0.7957
> log.lifespan                -0.212424 0.1406667 -1.510120  0.1337
> Endo_ectothermendo:log.mass -0.077842 0.0773762 -1.006014  0.3165
> 
>  Correlation: 
>                             (Intr) End_ct lg.mss lg.lfs
> Endo_ectothermendo          -0.371                     
> log.mass                    -0.410  0.395              
> log.lifespan                -0.220  0.000 -0.329       
> Endo_ectothermendo:log.mass  0.436 -0.453 -0.881  0.099
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.10748607 -0.15843091  0.07926082  0.26006154  0.86623002 
> 
> Residual standard error: 2.008683 
> Degrees of freedom: 122 total; 117 residual
intervals(pagel.model1, which = 'var-cov')
> Approximate 95% confidence intervals
> 
>  Residual standard error:
>    lower     est.    upper 
> 1.785163 2.008683 2.296694
output.models <- model.sel(pagel.model0.8, pagel.model0, pagel.model0.5, pagel.model1)

sel.table <- round(as.data.frame(output.models)[7:11], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]


knitr::kable(sel.table, digits = 3, caption = 'Model selection procedure based on AIC criterion', row.names = FALSE, align = 'l')
Model selection procedure based on AIC criterion
Model K logLik AICc delta weight
pagel.model0.5 6 -110.869 234.469 0.000 0.97
pagel.model0.8 6 -114.333 241.396 6.928 0.03
pagel.model0 6 -124.291 261.313 26.845 0.00
pagel.model1 6 -135.386 283.503 49.035 0.00 
pruned_dat$Endo_ectotherm <- as.factor(pruned_dat$Endo_ectotherm)

mod.plot <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass, correlation = corPagel(0.5, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(mod.plot)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   233.6552 247.6753 -111.8276
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> lambda 
>    0.5 
> 
> Coefficients:
>                                  Value Std.Error   t-value p-value
> (Intercept)                  2.7150448 0.3964989  6.847547  0.0000
> Endo_ectothermendo           0.3642954 0.4770102  0.763706  0.4466
> log.mass                     0.0186163 0.0381623  0.487818  0.6266
> Endo_ectothermendo:log.mass -0.0980766 0.0437340 -2.242573  0.0268
> 
>  Correlation: 
>                             (Intr) End_ct lg.mss
> Endo_ectothermendo          -0.594              
> log.mass                    -0.715  0.598       
> Endo_ectothermendo:log.mass  0.622 -0.679 -0.873
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -3.0262099 -0.4803355  0.1497991  0.6394544  2.3543177 
> 
> Residual standard error: 0.7483241 
> Degrees of freedom: 122 total; 118 residual
anova(mod.plot, pagel.model0.5) ## Given that these models do not differ from each other, we are going to use mod.plot for visualization purposes
>                Model df      AIC      BIC    logLik   Test  L.Ratio p-value
> mod.plot           1  5 233.6552 247.6753 -111.8276                        
> pagel.model0.5     2  6 233.7381 250.5622 -110.8690 1 vs 2 1.917158  0.1662
pred.mod <- data.frame(log.mass = seq(min(pruned_dat$log.mass), max(pruned_dat$log.mass), len = 122), Endo_ectotherm = pruned_dat$Endo_ectotherm, row.names = rownames(pruned_dat))
pred.mod$tl.lg <- predict(mod.plot, pred.mod, type = 'response', se = TRUE)$fit
pred.mod$se <- predict(mod.plot, pred.mod, type = 'response', se = TRUE)$se.fit
str(pred.mod)
> 'data.frame': 122 obs. of  4 variables:
>  $ log.mass      : num  0.631 0.778 0.925 1.072 1.219 ...
>  $ Endo_ectotherm: Factor w/ 2 levels "ecto","endo": 2 2 2 2 2 2 2 1 2 2 ...
>  $ tl.lg         : num  3.03 3.02 3.01 2.99 2.98 ...
>  $ se            : num  0.395 0.394 0.393 0.392 0.391 ...
ecto <- pred.mod[pruned_dat$Endo_ectotherm == 'ecto', ]
endo <- pred.mod[pruned_dat$Endo_ectotherm == 'endo', ]
pred.mod
>                              log.mass Endo_ectotherm    tl.lg        se
> Hirundo_rustica             0.6312718           endo 3.029179 0.3950648
> Accipiter_gentilis          0.7782917           endo 3.017497 0.3939497
> Amazona_amazonica           0.9253116           endo 3.005815 0.3928565
> Buteo_jamaicensis           1.0723315           endo 2.994132 0.3917853
> Buteo_swainsoni             1.2193514           endo 2.982450 0.3907363
> Chrysolophus_pictus         1.3663713           endo 2.970768 0.3897096
> Colinus_virginianus         1.5133912           endo 2.959086 0.3887055
> Dicentrarchus_labrax        1.6604111           ecto 2.745955 0.3539540
> Diomedea_exulans            1.8074310           endo 2.935721 0.3867656
> Falco_sparverius            1.9544510           endo 2.924039 0.3858302
> Haliaeetus_leucocephalus    2.1014709           endo 2.912357 0.3849180
> Macronectes_giganteus       2.2484908           endo 2.900674 0.3840293
> Macronectes_halli           2.3955107           endo 2.888992 0.3831640
> Meleagris_gallopavo         2.5425306           endo 2.877310 0.3823225
> Rissa_tridactyla            2.6895505           endo 2.865628 0.3815049
> Strigops_habroptila         2.8365704           endo 2.853945 0.3807114
> Tachymarptis_melba          2.9835903           endo 2.842263 0.3799420
> Pseudocaranx_wrighti        3.1306102           ecto 2.773325 0.3220499
> Cyprinus_carpio             3.2776301           ecto 2.776062 0.3192267
> Lutjanus_argentimaculatus   3.4246500           ecto 2.778799 0.3164778
> Upeneichthys_vlamingii      3.5716700           ecto 2.781536 0.3138051
> Acanthopagrus_schlegelii    3.7186899           ecto 2.784273 0.3112106
> Macquaria_ambigua           3.8657098           ecto 2.787010 0.3086963
> Nothobranchius_furzeri      4.0127297           ecto 2.789747 0.3062642
> Nothobranchius_rachovii     4.1597496           ecto 2.792484 0.3039162
> Gadus_morhua                4.3067695           ecto 2.795221 0.3016542
> Platycephalus_bassensis     4.4537894           ecto 2.797958 0.2994804
> Oncorhynchus_mykiss         4.6008093           ecto 2.800695 0.2973965
> Neoceratodus_forsteri       4.7478292           ecto 2.803432 0.2954044
> Urolophus_paucimaculatus    4.8948491           ecto 2.806169 0.2935061
> Callorhinchus_milii         5.0418690           ecto 2.808906 0.2917034
> Meriones_unguiculatus       5.1888890           endo 2.667029 0.3713942
> Chinchilla_lanigera         5.3359089           endo 2.655347 0.3710291
> Ondatra_zibethicus          5.4829288           endo 2.643665 0.3706901
> Mesocricetus_auratus        5.6299487           endo 2.631983 0.3703773
> Myocastor_coypus            5.7769686           endo 2.620300 0.3700908
> Marmota_monax               5.9239885           endo 2.608618 0.3698307
> Tamias_striatus             6.0710084           endo 2.596936 0.3695969
> Zalophus_californianus      6.2180283           endo 2.585254 0.3693895
> Castor_canadensis           6.3650482           endo 2.573571 0.3692086
> Ateles_geoffroyi            6.5120681           endo 2.561889 0.3690542
> Crocuta_crocuta             6.6590880           endo 2.550207 0.3689264
> Aplodontia_rufa             6.8061080           endo 2.538525 0.3688252
> Balaena_mysticetus          6.9531279           endo 2.526842 0.3687505
> Homo_sapiens                7.1001478           endo 2.515160 0.3687025
> Peromyscus_maniculatus      7.2471677           endo 2.503478 0.3686811
> Pteropus_rodricensis        7.3941876           endo 2.491796 0.3686863
> Saimiri_sciureus            7.5412075           endo 2.480113 0.3687182
> Ceratotherium_simum         7.6882274           endo 2.468431 0.3687767
> Giraffa_camelopardalis      7.8352473           endo 2.456749 0.3688617
> Hydrochoerus_hydrochaeris   7.9822672           endo 2.445067 0.3689734
> Pan_paniscus                8.1292871           endo 2.433384 0.3691117
> Pongo_pygmaeus              8.2763070           endo 2.421702 0.3692765
> Pan_troglodytes             8.4233269           endo 2.410020 0.3694677
> Equus_grevyi                8.5703469           endo 2.398338 0.3696855
> Myrmecophaga_tridactyla     8.7173668           endo 2.386655 0.3699296
> Tapirus_indicus             8.8643867           endo 2.374973 0.3702002
> Ursus_maritimus             9.0114066           endo 2.363291 0.3704970
> Chaetophractus_vellerosus   9.1584265           endo 2.351609 0.3708200
> Eschrichtius_robustus       9.3054464           endo 2.339926 0.3711693
> Loxodonta_africana          9.4524663           endo 2.328244 0.3715446
> Muntiacus_muntjak           9.5994862           endo 2.316562 0.3719460
> Muntiacus_reevesi           9.7465061           endo 2.304880 0.3723733
> Elephas_maximus             9.8935260           endo 2.293197 0.3728264
> Procavia_capensis          10.0405459           endo 2.281515 0.3733053
> Tupaia_tana                10.1875659           endo 2.269833 0.3738098
> Heterocephalus_glaber      10.3345858           endo 2.258151 0.3743399
> Macaca_mulatta             10.4816057           endo 2.246468 0.3748954
> Choloepus_hoffmanni        10.6286256           endo 2.234786 0.3754763
> Tursiops_truncatus         10.7756455           endo 2.223104 0.3760824
> Lemur_catta                10.9226654           endo 2.211422 0.3767135
> Ochotona_princeps          11.0696853           endo 2.199739 0.3773696
> Ailurus_fulgens            11.2167052           endo 2.188057 0.3780506
> Tadarida_brasiliensis      11.3637251           endo 2.176375 0.3787562
> Myotis_lucifugus           11.5107450           endo 2.164693 0.3794864
> Didelphis_virginiana       11.6577649           endo 2.153010 0.3802411
> Galegeeska_rufescens       11.8047849           endo 2.141328 0.3810200
> Macroscelides_proboscideus 11.9518048           endo 2.129646 0.3818230
> Atelerix_albiventris       12.0988247           endo 2.117964 0.3826500
> Lepus_californicus         12.2458446           endo 2.106281 0.3835008
> Oryctolagus_cuniculus      12.3928645           endo 2.094599 0.3843753
> Sylvilagus_aquaticus       12.5398844           endo 2.082917 0.3852733
> Panthera_tigris            12.6869043           endo 2.071235 0.3861946
> Sciurus_carolinensis       12.8339242           endo 2.059552 0.3871391
> Setifer_setosus            12.9809441           endo 2.047870 0.3881066
> Alligator_sinensis         13.1279640           ecto 2.959438 0.3521694
> Zootoca_vivipara           13.2749839           ecto 2.962175 0.3556599
> Liasis_fuscus              13.4220039           ecto 2.964912 0.3592041
> Trachemys_scripta          13.5690238           ecto 2.967649 0.3628005
> Taeniopygia_guttata        13.7160437           endo 1.989459 0.3932830
> Aphelocoma_coerulescens    13.8630636           endo 1.977777 0.3943848
> Aphelocoma_ultramarina     14.0100835           endo 1.966094 0.3955084
> Branta_leucopsis           14.1571034           endo 1.954412 0.3966535
> Calidris_alpina            14.3041233           endo 1.942730 0.3978200
> Calonectris_diomedea       14.4511432           endo 1.931048 0.3990078
> Cepphus_grylle             14.5981631           endo 1.919365 0.4002165
> Fregata_minor              14.7451830           endo 1.907683 0.4014460
> Fulmarus_glacialis         14.8922029           endo 1.896001 0.4026962
> Haematopus_ostralegus      15.0392228           endo 1.884319 0.4039668
> Larus_crassirostris        15.1862428           endo 1.872636 0.4052577
> Lonchura_striata           15.3332627           endo 1.860954 0.4065686
> Oceanodroma_leucorhoa      15.4802826           endo 1.849272 0.4078994
> Parus_major                15.6273025           endo 1.837590 0.4092498
> Passerculus_sandwichensis  15.7743224           endo 1.825907 0.4106197
> Pygoscelis_adeliae         15.9213423           endo 1.814225 0.4120089
> Riparia_riparia            16.0683622           endo 1.802543 0.4134171
> Sterna_hirundo             16.2153821           endo 1.790861 0.4148443
> Sula_sula                  16.3624020           endo 1.779178 0.4162901
> Tachycineta_albilinea      16.5094219           endo 1.767496 0.4177544
> Tachycineta_bicolor        16.6564418           endo 1.755814 0.4192371
> Turdus_merula              16.8034618           endo 1.744132 0.4207378
> Uria_lomvia                16.9504817           endo 1.732449 0.4222565
> Vireo_olivaceus            17.0975016           endo 1.720767 0.4237929
> Lacerta_agilis             17.2445215           ecto 3.036073 0.4658332
> Alligator_mississippiensis 17.3915414           ecto 3.038810 0.4703553
> Xiphophorus_maculatus      17.5385613           ecto 3.041547 0.4749006
> Xiphophorus_hellerii       17.6855812           ecto 3.044284 0.4794685
> Danio_rerio                17.8326011           ecto 3.047021 0.4840582
> Macaca_nemestrina          17.9796210           endo 1.650674 0.4333728
> Macaca_fascicularis        18.1266409           endo 1.638991 0.4350279
> Gorilla_gorilla            18.2736608           endo 1.627309 0.4366993
> Emys_orbicularis           18.4206808           ecto 3.057969 0.5026246
cols <- c('purple', 'orange')


##png('figure2.png', width = 7, height = 7, units = 'in', res = 360)


plot(log(Average_Telomere_Length_kb) ~ log.mass, data = pruned_dat, pch = 21, bg = cols[pruned_dat$Endo_ectotherm], col = 'black', las = 1, xlab = "Mass log(gr)", ylab = "Telomere length log(kb)", type = 'n')

grid(nx = NULL, ny = NULL, col = "lightgray", lwd = 1)
par(new = TRUE)

plot(log(Average_Telomere_Length_kb) ~ log.mass, data = pruned_dat, pch = 21, bg = cols[pruned_dat$Endo_ectotherm], col = 'black', las = 1, xlab = "Mass log(gr)", ylab = "Telomere length log(kb)")

lines(ecto$tl.lg ~ ecto$log.mass, data = ecto, lwd = 2.5, col = 'purple')
polygon(c(ecto$log.mass, rev(ecto$log.mass)), c(ecto$tl.lg+ecto$se, rev(ecto$tl.lg-ecto$se)), border = FALSE, col = alpha('purple', 0.2))
lines(endo$tl.lg ~ endo$log.mass, data = endo, lwd = 2.5, col = 'orange')
polygon(c(endo$log.mass, rev(endo$log.mass)), c(endo$tl.lg+endo$se, rev(endo$tl.lg-endo$se)), border = FALSE, col = alpha('orange', 0.2))

legend('bottomleft', legend = c('Ectotherms', 'Endotherms'), pch = 16, col = c('purple', 'orange'), lwd = 1, bty = 'n', cex = 0.8)

##dev.off()
pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 0] <- "absent"
pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 1] <- "present"
pruned_dat$Telomerase_activity[is.na(pruned_dat$Telomerase_activity)] <- "N/A"

tel.act <- setNames(pruned_dat$Telomerase_activity, rownames(pruned_dat))
tel.act
>            Hirundo_rustica         Accipiter_gentilis          Amazona_amazonica 
>                      "N/A"                      "N/A"                      "N/A" 
>          Buteo_jamaicensis            Buteo_swainsoni        Chrysolophus_pictus 
>                      "N/A"                      "N/A"                      "N/A" 
>        Colinus_virginianus       Dicentrarchus_labrax           Diomedea_exulans 
>                      "N/A"                      "N/A"                      "N/A" 
>           Falco_sparverius   Haliaeetus_leucocephalus      Macronectes_giganteus 
>                      "N/A"                      "N/A"                      "N/A" 
>          Macronectes_halli        Meleagris_gallopavo           Rissa_tridactyla 
>                      "N/A"                      "N/A"                      "N/A" 
>        Strigops_habroptila         Tachymarptis_melba       Pseudocaranx_wrighti 
>                      "N/A"                      "N/A"                      "N/A" 
>            Cyprinus_carpio  Lutjanus_argentimaculatus     Upeneichthys_vlamingii 
>                      "N/A"                      "N/A"                      "N/A" 
>   Acanthopagrus_schlegelii          Macquaria_ambigua     Nothobranchius_furzeri 
>                      "N/A"                      "N/A"                  "present" 
>    Nothobranchius_rachovii               Gadus_morhua    Platycephalus_bassensis 
>                      "N/A"                  "present"                      "N/A" 
>        Oncorhynchus_mykiss      Neoceratodus_forsteri   Urolophus_paucimaculatus 
>                      "N/A"                      "N/A"                      "N/A" 
>        Callorhinchus_milii      Meriones_unguiculatus        Chinchilla_lanigera 
>                      "N/A"                  "present"                  "present" 
>         Ondatra_zibethicus       Mesocricetus_auratus           Myocastor_coypus 
>                  "present"                  "present"                  "present" 
>              Marmota_monax            Tamias_striatus     Zalophus_californianus 
>                  "present"                  "present"                   "absent" 
>          Castor_canadensis           Ateles_geoffroyi            Crocuta_crocuta 
>                   "absent"                   "absent"                   "absent" 
>            Aplodontia_rufa         Balaena_mysticetus               Homo_sapiens 
>                   "absent"                   "absent"                   "absent" 
>     Peromyscus_maniculatus       Pteropus_rodricensis           Saimiri_sciureus 
>                  "present"                   "absent"                   "absent" 
>        Ceratotherium_simum     Giraffa_camelopardalis  Hydrochoerus_hydrochaeris 
>                   "absent"                   "absent"                  "present" 
>               Pan_paniscus             Pongo_pygmaeus            Pan_troglodytes 
>                   "absent"                   "absent"                      "N/A" 
>               Equus_grevyi    Myrmecophaga_tridactyla            Tapirus_indicus 
>                   "absent"                   "absent"                   "absent" 
>            Ursus_maritimus  Chaetophractus_vellerosus      Eschrichtius_robustus 
>                   "absent"                   "absent"                   "absent" 
>         Loxodonta_africana          Muntiacus_muntjak          Muntiacus_reevesi 
>                   "absent"                   "absent"                   "absent" 
>            Elephas_maximus          Procavia_capensis                Tupaia_tana 
>                   "absent"                   "absent"                   "absent" 
>      Heterocephalus_glaber             Macaca_mulatta        Choloepus_hoffmanni 
>                  "present"                   "absent"                   "absent" 
>         Tursiops_truncatus                Lemur_catta          Ochotona_princeps 
>                   "absent"                   "absent"                  "present" 
>            Ailurus_fulgens      Tadarida_brasiliensis           Myotis_lucifugus 
>                   "absent"                  "present"                  "present" 
>       Didelphis_virginiana       Galegeeska_rufescens Macroscelides_proboscideus 
>                  "present"                  "present"                  "present" 
>       Atelerix_albiventris         Lepus_californicus      Oryctolagus_cuniculus 
>                  "present"                   "absent"                   "absent" 
>       Sylvilagus_aquaticus            Panthera_tigris       Sciurus_carolinensis 
>                   "absent"                  "present"                  "present" 
>            Setifer_setosus         Alligator_sinensis           Zootoca_vivipara 
>                  "present"                   "absent"                      "N/A" 
>              Liasis_fuscus          Trachemys_scripta        Taeniopygia_guttata 
>                      "N/A"                      "N/A"                      "N/A" 
>    Aphelocoma_coerulescens     Aphelocoma_ultramarina           Branta_leucopsis 
>                      "N/A"                      "N/A"                      "N/A" 
>            Calidris_alpina       Calonectris_diomedea             Cepphus_grylle 
>                      "N/A"                      "N/A"                      "N/A" 
>              Fregata_minor         Fulmarus_glacialis      Haematopus_ostralegus 
>                      "N/A"                      "N/A"                      "N/A" 
>        Larus_crassirostris           Lonchura_striata      Oceanodroma_leucorhoa 
>                      "N/A"                      "N/A"                      "N/A" 
>                Parus_major  Passerculus_sandwichensis         Pygoscelis_adeliae 
>                      "N/A"                      "N/A"                      "N/A" 
>            Riparia_riparia             Sterna_hirundo                  Sula_sula 
>                      "N/A"                      "N/A"                      "N/A" 
>      Tachycineta_albilinea        Tachycineta_bicolor              Turdus_merula 
>                      "N/A"                      "N/A"                      "N/A" 
>                Uria_lomvia            Vireo_olivaceus             Lacerta_agilis 
>                      "N/A"                      "N/A"                      "N/A" 
> Alligator_mississippiensis      Xiphophorus_maculatus       Xiphophorus_hellerii 
>                      "N/A"                      "N/A"                      "N/A" 
>                Danio_rerio          Macaca_nemestrina        Macaca_fascicularis 
>                  "present"                      "N/A"                      "N/A" 
>            Gorilla_gorilla           Emys_orbicularis 
>                      "N/A"                      "N/A"
TA <- to.matrix(tel.act, unique(tel.act))
TA <- TA[pruned_data_tree$tip.label, ]

life.span <- setNames(pruned_dat$log.lifespan, rownames(pruned_dat))
log_mass <- setNames(pruned_dat$log.mass, rownames(pruned_dat))
telo.length <- setNames(log(pruned_dat$Average_Telomere_Length_kb), rownames(pruned_dat))
plotTree(pruned_data_tree, ftype = "off", mar = c(3, 2, 2, 3))

tiplabels(pie = TA, piecol = c("white", "gray", "black"), cex = 0.22, offset = 4.3)

par(xpd = TRUE)

legend(x = 150, y = 0.2, legend = unique(tel.act), pch = 21, pt.bg = c("white", "gray", "black"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)


par(new = TRUE)
par(mar = c(3, 32, 2, 1.1))
barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1, space = 0,
        ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)

axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)

## Ancestral state reconstruction of telomere size

fit <- fastAnc(pruned_data_tree, telo.length, vars = TRUE, CI = TRUE)
fit$CI[1, ]
> [1] 1.182278 4.654196
obj <- contMap(pruned_data_tree, telo.length, plot = FALSE)



Figure 1 

##png("figure1.png", width = 7, height = 7, units = "in", res = 360)
##pdf("figure1.pdf")


plot(obj, ftype = "off", legend = FALSE, ylim = c(1-0.09*(Ntip(obj$tree)-1), Ntip(obj$tree)), mar = c(1, 0.1, 1, 4), lwd = 1.5)
add.color.bar(150, obj$cols,title = "Log telomere length (kb)", lims = obj$lims, digits = 3, prompt = FALSE, x = 0, y = 1-0.08*(Ntip(obj$tree)-1), lwd = 4, fsize = 0.6, subtitle = "")

par(xpd = TRUE)

text(300, 118, 'Testudines')
text(300, 100.5, 'Aves', col = 'red')
text(300, 82, 'Crocodilia')
text(240, 73, 'Squamata')
text(240, 39, 'Mammalia', col = 'red')
text(160, 15, 'Osteichthyes')
text(160, 4.5, 'Chondrichthyes')


#cladelabels(pruned_data_tree, node = 217, "Ectotherms", offset = 6)
#cladelabels(pruned_data_tree, node = 153, "Endotherms", offset = 6)
#segments(518, 0, 518, 18)
#segments(518, 0, 508, 0)
#segments(518, 18, 508, 18)
#text(528, 9, 'Ectotherms', srt = 90)

tiplabels(pie = TA, piecol = c("gray", "white", "black"), cex = 0.16, offset = 5.7)
legend(x = 200, y = -5, legend = unique(tel.act), pch = 21, pt.bg = c("gray", "white", "black"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)

par(new = TRUE)
par(mar = c(3.6, 30.8, 3, 2.2))

barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1.07, space = 0,
        ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)

axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Log lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)

##dev.off()
## Correlated evolution under the threshold model


## Removing NAs from the dataset

pruned_dat$Telomerase_activity
>   [1] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
>  [11] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
>  [21] "N/A"     "N/A"     "N/A"     "present" "N/A"     "present" "N/A"     "N/A"     "N/A"     "N/A"    
>  [31] "N/A"     "present" "present" "present" "present" "present" "present" "present" "absent"  "absent" 
>  [41] "absent"  "absent"  "absent"  "absent"  "absent"  "present" "absent"  "absent"  "absent"  "absent" 
>  [51] "present" "absent"  "absent"  "N/A"     "absent"  "absent"  "absent"  "absent"  "absent"  "absent" 
>  [61] "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "present" "absent"  "absent"  "absent" 
>  [71] "absent"  "present" "absent"  "present" "present" "present" "present" "present" "present" "absent" 
>  [81] "absent"  "absent"  "present" "present" "present" "absent"  "N/A"     "N/A"     "N/A"     "N/A"    
>  [91] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
> [101] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
> [111] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "present" "N/A"     "N/A"    
> [121] "N/A"     "N/A"
pruned_dat_not_NAs <- pruned_dat[!pruned_dat$Telomerase_activity == "N/A", ]
pruned_dat_not_NAs$Telomerase_activity
>  [1] "present" "present" "present" "present" "present" "present" "present" "present" "present" "absent" 
> [11] "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "present" "absent"  "absent"  "absent" 
> [21] "absent"  "present" "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "absent" 
> [31] "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "present" "absent"  "absent"  "absent" 
> [41] "absent"  "present" "absent"  "present" "present" "present" "present" "present" "present" "absent" 
> [51] "absent"  "absent"  "present" "present" "present" "absent"  "present"
check2 <- name.check(pruned_data_tree, pruned_dat_not_NAs)
rm_phy2 <- check2$tree_not_data
rm_dat2 <- check2$data_not_tree
pruned_data_tree_not_NAs <- drop.tip(pruned_data_tree, rm_phy2)
pruned_dat_not_NAs <- subset(pruned_dat_not_NAs, subset = pruned_dat_not_NAs$Scientific_name %in% pruned_data_tree_not_NAs$tip, select = names(pruned_dat_not_NAs))

pruned_dat_not_NAs$Telomerase_activity <- as.factor(pruned_dat_not_NAs$Telomerase_activity)
name.check(pruned_data_tree_not_NAs, pruned_dat_not_NAs)
> [1] "OK"
## Set the number of generations

ngen <- 5e6

## Run MCMC

pruned_dat_not_NAs$Average_Telomere_Length_kb <- log(pruned_dat_not_NAs$Average_Telomere_Length_kb)

mcmc.model <- threshBayes(pruned_data_tree_not_NAs, pruned_dat_not_NAs[ , c(14, 15)], type = c("cont", "disc"), ngen = ngen, plot = FALSE, control = list(print.interval = 5e+05))
> Starting MCMC....
> generation: 500000; mean acceptance rate: 0.23
> generation: 1000000; mean acceptance rate: 0.26
> generation: 1500000; mean acceptance rate: 0.23
> generation: 2000000; mean acceptance rate: 0.22
> generation: 2500000; mean acceptance rate: 0.2
> generation: 3000000; mean acceptance rate: 0.24
> generation: 3500000; mean acceptance rate: 0.24
> generation: 4000000; mean acceptance rate: 0.23
> generation: 4500000; mean acceptance rate: 0.22
> generation: 5000000; mean acceptance rate: 0.24
> Done MCMC.
## Pull out the post burn-in sample and compute HPD

r.mcmc <- tail(mcmc.model$par$r, 0.8*nrow(mcmc.model$par))
class(r.mcmc) <- "mcmc"

hpd.r <- HPDinterval(r.mcmc)
hpd.r
>          lower     upper
> var1 0.1833235 0.7869591
> attr(,"Probability")
> [1] 0.9500012



Figure S3 

## Profile plots from a Bayesian MCMC analysis of the threshold model

#png("figureS3.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS3.pdf")



layout(matrix(c(0, 0, 0, 0,
                1, 1, 1, 1,
                1, 1, 1, 1,
                2, 2, 2, 2,
                2, 2, 2, 2,
                0, 0, 0, 0), nrow = 6, ncol = 4, byrow = TRUE))

plot(mcmc.model$par$logL ~ mcmc.model$par$gen, col = 'gray', type = 'n', ylab = 'Log(L)', xlab = 'Generation', las = 1)
grid()
par(new = TRUE)
plot(mcmc.model$par$logL ~ mcmc.model$par$gen, col = 'gray', type = 'l', ylab = 'Log(L)', xlab = 'Generation', las = 1)
mtext('(a)', side = 2, line = 1.5, at = -160, las = 1)

plot(mcmc.model$par$r ~ mcmc.model$par$gen, col = 'purple', type = 'n', ylab = 'r', xlab = 'Generation', las = 1)
grid()
par(new = TRUE)
plot(mcmc.model$par$r ~ mcmc.model$par$gen, col = 'purple', type = 'l', ylab = 'r', xlab = 'Generation', las = 1)
mtext('(b)', side = 2, line = 1.5, at = 1.2, las = 1)

#dev.off()



Figure 4 

## Plot posterior density

#png("figure4.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure4.pdf")


par(las = 1)

plot(NA, xlim = c(-1, 1.5), ylim = c(0, 3), axes = FALSE, xlab = '', ylab = '')
grid()

par(new = TRUE)
plot(density(mcmc.model), xlim = c(-1, 1.5))
box()

## add whiskers to show HPD

h <- 0-par()$usr[3]
lines(x = hpd.r, y = rep(-h/2, 2))
lines(x = rep(hpd.r[1], 2), y = c(-0.3, -0.7)*h)
lines(x = rep(hpd.r[2], 2), y = c(-0.3, -0.7)*h)

#dev.off()