# Loading species morphometric data spp_data<-read.table("spp_morphom.txt",header=T) tree<-read.tree("brachymeles_tree.txt") # Reducing data to only variables included in PCAs spp_data_ppca<-data.frame(spp_data[,2:3],spp_data[,14:17]) spp_data_ppca_limb<-data.frame(spp_data[,2:3],spp_data[,14:16],spp_data[,18:21]) # Transformation of data for PCA spp_data_ppca$HD<-spp_data_ppca$HD+2 spp_data_ppca[,2:6]<-spp_data_ppca[,2:6]+0.01 spp_data_ppca<-log(spp_data_ppca) spp_data_ppca_limb$HD<-spp_data_ppca_limb$HD+2 spp_data_ppca_limb[,2:9]<-spp_data_ppca_limb[,2:9]+0.01 spp_data_ppca_limb<-log(spp_data_ppca_limb) # Running the PCA, three different versions ppca_cor_lam<-phyl.pca(tree,spp_data_ppca,method="lambda",mode="corr") ppca_cov_lam<-phyl.pca(tree,spp_data_ppca,method="lambda",mode="cov") ppca_cor_lam_limb<-phyl.pca(tree,spp_data_ppca_limb,method="lambda",mode="corr") # Loading and transforming individual level morphometrics in same way as species level data indiv_data<-read.table("indiv_morphom.txt",header=T) ind_data_for_ppca<-data.frame(indiv_data[,6:7],indiv_data[,18:21]) ind_data_for_ppca$HD<-ind_data_for_ppca$HD+2 ind_data_for_ppca[,2:6]<-ind_data_for_ppca[2:6]+0.01 # Repeats pPCA on species data with correlation matrix pca<-phyl.pca(tree,spp_data_ppca,method="lambda",mode="corr") # Calculates factor scores for individuals obj<-phyl.vcv(as.matrix(spp_data_ppca)[tree$tip.label,],vcv(tree),pca$lambda) as.matrix(spp_data_ppca)[tree$tip.label,] as.matrix(spp_data_ppca) V<-obj$R C<-obj$C nrow(spp_data_ppca) ncol(spp_data_ppca) n<-nrow(spp_data_ppca) m<-ncol(spp_data_ppca) Y<-spp_data_ppca/matrix(rep(sqrt(diag(V)),n),n,m,byrow=TRUE) diag(phyl.vcv(as.matrix(Y)[tree$tip.label,],C,1)$R) invC<-solve(C)[rownames(Y),rownames(Y)] a<-matrix(colSums(invC%*%as.matrix(Y))/sum(invC),m,1) A<-matrix(rep(a,n),n,m,byrow=TRUE) Y<-as.matrix(Y)-A phyl.vcv(as.matrix(Y)[tree$tip.label,],C,1)$alpha[,1] Sm<-Y%*%pca$Evec plot(pca$S,Sm,pch=21,bg=make.transparent("blue",0.5),cex=1.5,xlab="PC scores from phyl.pca",ylab="PC scores from this exercise") n<-nrow(ind_data_for_ppca) Yi<-ind_data_for_ppca/matrix(rep(sqrt(diag(V)),n),n,m,byrow=TRUE) A<-matrix(rep(a,n),n,m,byrow=TRUE) Yi<-as.matrix(Yi)-A Si<-Yi%*%pca$Evec Si->ind_ppca_scores # Supplementary Figure S1 - phylomorphopace for PC1 and PC2 with individuals plotted on # spp_col and ind_col are variables in the species and individual datasets that set the color for each point in the figure plot(ind_ppca_scores[,1],ind_ppca_scores[,2],pch=21,col="black",bg=ind_col,cex=1,cex.axis=1.5,xlim=c(-30,40),ylim=c(-18,15)) phylomorphospace(tree,cbind(ppca_cor_lam$S[,1],ppca_cor_lam$S[,2]),label="off",node.size=c(0,0),lwd=2,add=T) points(ppca_cor_lam$S[,1],ppca_cor_lam$S[,2],pch=23,col="black",bg=spp_col,cex=2,cex.axis=1.5,xlim=c(-30,40),ylim=c(-18,15)) # Sample plot of bivariate phylomoprhospace for two morphometric variables plot(spp_bur_f_data$PC1,spp_bur_f_data$D_slope,pch=21,col="black",bg=colors,cex=2,cex.axis=1.5,xlim=c(-25,35),ylim=c(0.30,0.60)) phylomorphospace(tree,cbind(td_spp_bF$data[,2],td_spp_bF$data[,10]),label="off",node.size=c(0,0),lwd=1,xlim=c(-25,35),ylim=c(0.3,0.6),add=T) points(spp_bur_f_data$PC1,spp_bur_f_data$D_slope,pch=21,col="black",bg=colors,cex=2,xlim=c(-25,35),ylim=c(0.3,0.6)) abline(coef=c(0.475,-0.001),col="black",lwd=6,lty=2) # PGLS analysis of species data # Code prepares the data for analysis and conducts a PGLS regression while estimating lambda library(geiger) library(caper) cp_spp_bur_f<-comparative.data(tree,spp_bur_f_data,species,na.omit=F) pgls_Vavg_F10<-pgls(Fmax_10~Vavg_10,cp_spp_bur_f,lambda="ML",bounds=list(lambda=c(0,1))) pgls_Vavg_F60<-pgls(Fmax_60~Vavg_60,cp_spp_bur_f,lambda="ML",bounds=list(lambda=c(0,1))) # Preparation of objects and generation of a phylomorphospace plot on two substrates # Plots two phylomorphospaces with different symbols for coarse and fine substrates, with species color coded by number of digits # Phylomorphospaces are ploted in the plot area, then the symbols, then the regression lines when the slope is significant # 10 and 60 in the various variables represent coarse and fine substrates, respectively library(geiger) library(phytools) library(caper) smp_tree_10<-paintSubTree(td_smp_bur_Luse10$phy,Ntip(td_smp_bur_Luse10$phy)+1,"1") smp_tree_60<-paintSubTree(td_smp_bur_Luse60$phy,Ntip(td_smp_bur_Luse60$phy)+1,"1") td_spp_bur_Luse10<-treedata(tree,spp_bur_Luse10,sort=T) td_smp_bur_Luse10<-treedata(smp_tree,spp_bur_Luse10,sort=T) td_spp_bur_Luse60<-treedata(tree,spp_bur_Luse60,sort=T) td_smp_bur_Luse60<-treedata(smp_tree,spp_bur_Luse60,sort=T) colors_10<-c("#ca0020","#ca0020","#0571b0","#0571b0","#f4a582","#ca0020","#ca0020","#ca0020","black") colors_60<-c("#92c5de","#ca0020","#ca0020","#ca0020","#0571b0","#0571b0","#f4a582","#ca0020","#ca0020","#ca0020","black") plot(spp_bur_Luse60$PC1,spp_bur_Luse60$FLuse_60,pch=21,col="black",bg=colors_60,cex=2,cex.axis=1.5,xlim=c(-25,35),ylim=c(0,1)) phylomorphospace(td_spp_bur_Luse60$phy,cbind(td_spp_bur_Luse60$data[,3],td_spp_bur_Luse60$data[,18]),label="off",node.size=c(0,0),lwd=1,add=T) phylomorphospace(smp_tree_10,cbind(td_smp_bur_Luse10$data[,3],td_smp_bur_Luse10$data[,10]),label="off",node.size=c(0,0),lwd=1,add=T,colors=setNames("red","1")) points(spp_bur_Luse60$PC1,spp_bur_Luse60$FLuse_60,pch=21,col="black",bg=colors_60,cex=2) points(spp_bur_Luse10$PC1,spp_bur_Luse10$FLuse_10,pch=22,col="black",bg=colors_10,cex=2) abline(coef=c(0.6989161,-0.0171536),col="black",lwd=6,lty=2) abline(coef=c(0.6395325,-0.0189446),col="red",lwd=6,lty=2) # Fitting of trait evolution models to head length data (Brachymeles_data) using the 40 species tree # First tree is matched with data, regimes for the tip species are defined and a dataset assembled # Then simmap trees are produced to represent selective regimes with two and three optima library(geiger) library(phytools) library(OUwie) tree<-read.tree("brachymeles_tree_use.txt") tree<-ladderize(tree) td_brachy<-treedata(tree,Brachymeles_data,sort=T) td_brachy$data[,4]->HL_data reg3<-c(0,1,2,2,1,2,1,0,1,1,1,2,2,2,1,0,2,2,0,0,1,2,1,2,2,1,2,1,0,0,0,1,0,0,2,2,2,1,2,2) reg2<-c(0,0,1,1,0,1,0,0,0,0,0,1,1,1,0,0,1,1,0,0,0,1,0,1,1,0,1,0,0,0,0,0,0,0,1,1,1,0,1,1) B_dat_use<-data.frame(reg2,reg3,Brachymeles_data$HL) colnames(B_dat_use)[3]<-"HL" rownames(B_dat_use)<-rownames(Brachymeles_data) HL_dat2<-data.frame(B_dat_use[,1],B_dat_use[,3]) colnames(HL_dat2)<-c("reg2","HL") rownames(HL_dat2)<-rownames(B_dat_use) HL_dat3<-B_dat_use[,2:3] HL_dat2<-data.frame(rownames(HL_dat2),HL_dat2) HL_dat3<-data.frame(rownames(HL_dat3),HL_dat3) colnames(HL_dat2)[1]<-"sp" colnames(HL_dat3)[1]<-"sp" tree_regime<-paintSubTree(tree,41,0,stem=F) tree_regime<-paintSubTree(tree_regime,41,0,stem=F) plot(tree_regime) tree_OU2<-paintSubTree(tree_regime,62,"1",stem=T) plot(tree_OU2) tree_OU3<-paintSubTree(tree_regime,62,"2",stem=T) tree_OU3<-paintSubTree(tree_OU3,44,1,stem=T) tree_OU3<-paintSubTree(tree_OU3,50,0,stem=T) plot(tree_OU3) mBM1_HL<-OUwie(tree_OU2,HL_dat2,model="BM1",simmap.tree=T) mOU1_HL<-OUwie(tree_OU2,HL_dat2,model="OU1",simmap.tree=T) mOU2_HL<-OUwie(tree_OU2,HL_dat2,model="OUM",simmap.tree=T) mOU3_HL<-OUwie(tree_OU3,HL_dat3,model="OUM",simmap.tree=T)