##----------------------------------------------------------------------------- ##Larval biology of the octocoral species Dentomuricea aff. meteor ##under different temperature regimes ## Rakka M., Godinho A., Orejas C., Carreiro-Silva M. ## Code and statistical analysis ##Creation Date:7/4/2020 ##Last update:3/3/2021 ##----------------------------------------------------------------------------- #Libraries library(tidyverse) library(ggplot2) library(gridExtra) library(visreg) library(Cairo) library(readxl) library(ggpubr) ##---------------------------------------------------------------------------- ##Read Data ##---------------------------------------------------------------------------- path<-"MRakka_etal_LarvalBiologyDmeteor_Data.xlsx" devdat<-read_excel(path,sheet=1) str(devdat) mylevs<-c("unfertilized egg","fertilized egg","2-cell","4-cell","8-cell","16-cell","32-cell", "64-cell","blastula","gastrula","late gastrula","planula","anomalous") devdat<-devdat %>% filter(!is.na(Stage2))%>% mutate(mystage=as_factor(Stage2), mystage=fct_relevel(mystage,mylevs), mystage=recode(mystage,"gastrula"="early gastrula")) ##---------------------------------------------------------------------------- ##Developmental stages ##---------------------------------------------------------------------------- hourlevs<-c("1-2","3-6","7-10","16-24","40-48","64-72","96") totcount<-devdat %>% mutate(Temp=as_factor(Temperature), Temp=recode(Temp,"13"="13°C","15"="15°C"), HourFromSpawn2=case_when(HourFromSpawn%in%c(1:2)~"1-2", HourFromSpawn%in%c(3:6)~"3-6", HourFromSpawn%in%c(7:10)~"7-10", HourFromSpawn%in%c(16:32)~"16-24", HourFromSpawn%in%c(40:52)~"40-48", HourFromSpawn%in%c(64:76)~"64-72", TRUE~"96"), HourFromSpawn2=fct_relevel(HourFromSpawn2,hourlevs))%>% select(Batch,Temp,HourFromSpawn2,mystage) %>% filter(!mystage%in%c(NA,"anomalous","unfertilized egg")) %>% count(Temp,HourFromSpawn2,mystage) %>% group_by(Temp,HourFromSpawn2) %>% mutate(totcount=sum(n))%>% mutate(perstage=n/totcount) %>% unique() fig1<-ggplot(totcount,aes(y=perstage,x=as.factor(HourFromSpawn2)))+ geom_bar(aes(fill=mystage),stat="identity")+ facet_wrap(~Temp,ncol=1,scales="free")+ scale_fill_viridis_d(direction=1)+ theme_bw()+ theme(legend.title=element_blank(), strip.background=element_rect(colour="black",fill="white"), axis.text.x = element_text(angle=45,vjust=0.9,hjust=0.8))+ ylab("Percentage(%)")+xlab("Time from spawning (h)")+ geom_label(aes(label=totcount,y=0),stat='identity', size=2.5) pdf(file="Graphs/Figure2.pdf", width=4.5,height=5) print(fig1) dev.off() ##---------------------------------------------------------------------------- ##Embryo and larval size ##---------------------------------------------------------------------------- devdat<-devdat%>% rename("Diam_small"="Diameter(small)","Diam_big"="Diameter(big)")%>% mutate(vol=4*3.14*(Diam_big/2)^2*(Diam_small/2)/3, lwratio=Diam_big/Diam_small) avtab<-devdat%>% group_by(Temperature,mystage)%>% summarize(volmean=mean(vol,na.rm=T), volsd=sd(vol,na.rm=T), lwmean=mean(lwratio,na.rm=T), lwsd=sd(lwratio,na.rm=T)) devdat%>%filter(mystage%in%levels(mystage)[1:6])%>% summarize(avsize=mean(vol,na.rm=TRUE), sdsize=sd(vol,na.rm=TRUE)) devdat%>%filter(mystage=="planula")%>% group_by(HourFromSpawn)%>% summarize(avsize=mean(vol,na.rm=TRUE), sdsize=sd(vol,na.rm=TRUE)) mylevs<-c("unfertilized egg","fertilized egg","2-cell","4-cell","8-cell","16-cell","32-cell", "64-cell","blastula","early gastrula","late gastrula","planula","planula_96","planula_336","anomalous") sizvol<-devdat%>% mutate(mystage2=as_factor(ifelse(HourFromSpawn==96|HourFromSpawn==336, paste(mystage,HourFromSpawn,sep="_"), paste(mystage))), Temperature=as_factor(Temperature), mystage2=fct_relevel(mystage2,mylevs), mystage2=recode(mystage2,"planula_96"="planula 4","planula_336"="planula 14")) thelmfun<-function(y,mydat){ mod0<-lm(as.formula(paste(y,"~1")),data=mydat) mod1<-lm(as.formula(paste(y,"~","mystage2")),data=mydat) mod2<-lm(as.formula(paste(y,"~","mystage2","+Temperature")),data=mydat) mod3<-lm(as.formula(paste(y,"~","mystage2","*Temperature")),data=mydat) myan=anova(mod0,mod1,mod2,mod3,test="Chisq") bestmod<-if(myan$`Pr(>Chi)`[2]<0.05){mod1}else{mod0} aicvals<-AIC(mod0,mod1,mod2,mod3) bestsum=summary(mod1) bestplots<-plot(bestmod) return(list(myan,aicvals,bestsum,bestplots)) } thelmfun(y="vol",mydat=sizvol) thelmfun(y="lwratio",mydat=sizvol) # Model Results graphically modsize<-lm(vol~mystage2,data=sizvol) vissize<-visreg(modsize,plot=FALSE) volplot<-ggplot() + geom_jitter(data=vissize$res,aes(mystage2, visregRes))+ geom_boxplot(data=vissize$fit,aes(y = visregFit,x=as.factor(mystage2)))+ geom_errorbar(data=vissize$fit,aes(ymin=visregLwr,ymax=visregUpr,x=as.factor(mystage2)))+ theme_light()+ theme(legend.title=element_blank(),legend.position=c(0.9,0.4), legend.background = element_blank(), axis.text.x=element_text(angle=90))+ ylab(bquote('Volume ('*'m'~m^-1*')'))+xlab("") modvol<-lm(lwratio~mystage2*Temperature,data=sizvol) vislw<-visreg(modvol,xvar="mystage2",by="Temperature",plot=F) lwplot<-ggplot() + geom_jitter(data=vislw$res,aes(mystage2, visregRes))+ geom_boxplot(data=vislw$fit,aes(y = visregFit,x=as.factor(mystage2)))+ geom_errorbar(data=vislw$fit,aes(ymin=visregLwr,ymax=visregUpr,x=as.factor(mystage2)))+ facet_wrap(~Temperature,ncol=1)+ theme_light()+ theme(legend.title=element_blank(),legend.position=c(0.9,0.4), legend.background = element_blank(), axis.text.x=element_text(angle=90), strip.background = element_rect(fill="white",colour="gray"), strip.text=element_text(colour="black"))+ ylab("Length/width ratio")+xlab("") pdf(file="Graphs/FigureS1.pdf", width=5.5,height=4) ggarrange(volplot,lwplot,ncol=2,labels=c("A","B")) dev.off() ##---------------------------------------------------------------------------- ##Survival analysis ##---------------------------------------------------------------------------- library(survival) library(survminer) survdat<-read_excel(path,sheet=2)%>% mutate(dead=ifelse(Trunc==0,1,0), Batch=as_factor(Batch)) survdat$Temperature <- factor(survdat$Temperature, levels = c("13", "15"), labels = c("13°C", "15°C")) survdat$dead<-ifelse(survdat$Trunc==0,1,0) surv_object <- Surv(time = survdat$LastDayAlive,event = survdat$dead) fit1 <- survfit(surv_object ~ Temperature, data = survdat) summary(fit1) ggsurvplot(fit1, data = survdat) survdiff(Surv(LastDayAlive,dead)~ Temperature, data = survdat) pooled<-ggsurvplot(fit1, conf.int=TRUE, pval=FALSE, legend.labs=c("13°C", "15°C"),legend.title="Temperature", palette=c("gray32", "gray60"), risk.table.height=.25,xlab="Time (days)",risk.table = FALSE) pdf(file="Graphs/Figure3.pdf", width=5,height=4.5) pooled dev.off() survdat13<-survdat[survdat$Temperature=="13°C",] surv_13 <- Surv(time = survdat13$LastDayAlive,event = survdat13$dead) fit2 <- survfit(surv_13 ~ Batch, data = survdat13) summary(fit2) ggsurvplot(fit2, data = survdat13, pval = TRUE) survdat15<-survdat[survdat$Temperature=="15°C",] surv_15 <- Surv(time = survdat15$LastDayAlive,event = survdat15$dead) fit3 <- survfit(surv_15 ~ Batch, data = survdat15) summary(fit3) ggsurvplot(fit3, data = survdat15, pval = TRUE) #Combine fits surv_object13 <- Surv(time = survdat13$LastDayAlive,event = survdat13$dead) fit13 <- survfit(surv_object13 ~ 1, data = survdat13) surv_object15 <- Surv(time = survdat15$LastDayAlive,event = survdat15$dead) fit15<-survfit(surv_object15 ~ 1, data = survdat15) fitall13 <- list(batches = fit2, pooled = fit13) temp13<-ggsurvplot_combine(fitall13, survdat,palette=c("gray60","gray61","gray62","gray64", "gray65","black"),legend="none",xlab="Days",title="13°C") fitall15<-list(batches = fit3, pooled = fit15) temp15<- ggsurvplot_combine(fitall15, survdat,palette=c("gray60","gray61","gray62","gray64", "gray65","black"),legend="none",xlab="Days",title="15°C",ylab="") tog<-list(temp13,temp15) figs1<-arrange_ggsurvplots(tog,nrow=1,print="TRUE") pdf(file="Graphs/FigureS2.pdf", width=5.5,height=3) print(figs1) dev.off() ##---------------------------------------------------------------------------- ##Larval stages ##---------------------------------------------------------------------------- library(reshape2) library(tidyverse) library(VIM) library(tidyr) library(lmtest) library(mgcv) mordat<-read_excel(path,sheet=3) mordat4<-mordat%>% select(Batch,Temperature,DaysFromSpawn,Count,Planula:Deformed) allmet<-pivot_longer(mordat4, cols=c("Planula","Metamorphosed","Settled","Deformed"), names_to="stage")%>% mutate(value=replace_na(value,0))%>% droplevels()%>% complete(nesting(Temperature,Batch),DaysFromSpawn=1:39,stage)%>% group_by(stage,Batch)%>% nest()%>% mutate(newdat=map(data,~approx(.x$DaysFromSpawn,.x$value,xout=.x$DaysFromSpawn)))%>% unnest_wider(newdat)%>% unnest(cols=c(data,x,y))%>% ungroup()%>% mutate(y=ifelse(y<1,0,y))%>% group_by(Batch,Temperature,DaysFromSpawn)%>% mutate(plann=sum(y, na.rm=TRUE),perc=y/plann, stage=fct_relevel(stage,"Planula","Metamorphosed","Settled","Deformed"))%>% ungroup()%>% mutate(Temperature=recode(Temperature,"13"="13°C","15"="15°C")) allmet%>% filter(DaysFromSpawn>=2)%>% group_by(Batch,Temperature)%>% mutate(inlar=ifelse(Temperature=="13°C",sum(y[DaysFromSpawn==2]),sum(y[DaysFromSpawn==3])))%>% mutate(totper=y/inlar)%>% group_by(Temperature,stage,DaysFromSpawn)%>% summarize(mean(totper,na.rm=TRUE))%>% filter(DaysFromSpawn==39) allmet%>% filter(stage=="Metamorphosed")%>% group_by(Temperature,DaysFromSpawn)%>% summarize(yall=sum(y,na.rm=T))%>% group_by(Temperature)%>% summarize(max(yall,na.rm=T)) figall<-allmet%>% filter(DaysFromSpawn>=2)%>% mutate()%>% group_by(DaysFromSpawn,Temperature)%>% summarize(survlar=sum(y,na.rm=TRUE))%>% ggplot(aes(y=survlar,x=DaysFromSpawn))+ geom_vline(xintercept = 36,linetype="dashed",colour="gray71")+ #geom_errorbar(aes(ymin=avsurv-sdsurv,ymax=avsurv+sdsurv),colour="gray")+ geom_point(aes(colour=as.factor(Temperature)))+ theme_bw()+ylab("Total larvae")+ scale_x_continuous(limits=c(3,40),breaks=seq(0, 100, by = 10))+ theme(legend.title=element_blank(), strip.background = element_rect(colour="black",fill="white"), legend.position=c(0.8,0.9), legend.background = element_blank(), legend.key.size = unit(0.2, 'cm'), axis.title.x = element_blank())+ scale_colour_manual(values=c("black", "gray60")) allnewmeans<-allmet%>% group_by(DaysFromSpawn,stage,Temperature)%>% summarize(mymean=mean(perc,na.rm=T),mysd=sd(perc,na.rm=T))%>% mutate(pery=mymean*100,persd=mysd*100) fig4<-ggplot(allnewmeans,aes(y=mymean,x=DaysFromSpawn))+ geom_vline(xintercept = 36,linetype="dashed",colour="gray71")+ geom_errorbar(aes(ymin=mymean-mysd,ymax=mymean+mysd),colour="gray")+ geom_point(aes(colour=as.factor(Temperature)))+ facet_wrap(~stage,ncol=1,scales="free")+ theme_bw()+xlab("Days")+ylab("Proportion of larvae (%)")+ scale_x_continuous(limits=c(3,40))+ theme(legend.title=element_blank(), strip.background = element_rect(colour="black",fill="white"), legend.position="none")+ scale_colour_manual(values=c("black", "gray60")) pdf(file="Graphs/Figure4.pdf", width=5.5,height=7) ggarrange(figall,fig4,heights = c(0.5, 2),align = "v",ncol=1,labels=c("A","B"), font.label = list(size = 11)) dev.off() ## Statistical analysis of different larval stages joinbatch<-allmet%>% filter(DaysFromSpawn>=2)%>% mutate(Temperature=as_factor(Temperature))%>% rename(props=y)%>% group_by(stage)%>% nest() bingam<-function(mydat){ m0<-gam(perc~1,data=mydat,family=binomial,method = 'REML',weights = plann) m1<-gam(perc~s(DaysFromSpawn,k=4),data=mydat,family=binomial,method = 'REML',weights = plann) m2<-gam(perc~s(DaysFromSpawn,k=4)+Temperature,data=mydat,family=binomial,method = 'REML',weights = plann) m3<-gam(perc~s(DaysFromSpawn, by=Temperature,k=4),data=mydat,family=binomial,method = 'REML',weights = plann) aictab<-AIC(m0,m1,m2,m3) anovtab<-lrtest(m0,m1,m2,m3) return(list(anovtab,aictab)) } map(joinbatch$data,~bingam(.)) models<-joinbatch%>%mutate(glms=map2(stage, data, ~ if(.x=="Planula"){gam(perc~s(DaysFromSpawn,k=4)+Temperature, data=.y,family=binomial,weights = plann)}else{ gam(perc~s(DaysFromSpawn,k=4,by=Temperature),data=.y, family=binomial,weights = plann)}), forms=paste(map(glms,"formula"))) map(models$glms,~summary(.)) #Model results makevis<-function(x,y,z){if(grepl("Temperature",x)){ visreg(y,xvar="DaysFromSpawn",by="Temperature", overlay=TRUE,plot=FALSE,data=z)}else{visreg(y,xvar="DaysFromSpawn", overlay=TRUE,plot=FALSE,data=z)}} visregs<-mapply(x=models$forms,y=models$glms,z=models$data,makevis,SIMPLIFY=FALSE) visplots<-lapply(visregs,function(x) if("Temperature"%in%names(x[["fit"]])){ ggplot(x[["fit"]],aes(DaysFromSpawn, visregFit, linetype=factor(Temperature), fill=factor(Temperature))) + geom_ribbon(aes(ymin=visregLwr, ymax=visregUpr), alpha=0.5, colour="grey50", linetype=1, size=0.2) + geom_line()+ scale_fill_grey(start=0.5, end=0.8)+ theme_light()+ theme(legend.title=element_blank(),legend.position=c(0.85,0.45), legend.background = element_blank())+ ylab("Proportion of larvae")+xlab("Days from spawning")}else{ ggplot(x[["fit"]],aes(DaysFromSpawn, visregFit)) + geom_ribbon(aes(ymin=visregLwr, ymax=visregUpr), alpha=0.5, colour="grey50", linetype=1, size=0.2) + geom_line()+ scale_fill_grey(start=0.5, end=0.8)+ theme_light()+ theme(legend.position="none")+ ylab("Proportion of larvae")+xlab("Days from spawning") }) pdf(file="Graphs/FigureS5.pdf", width=5.7,height=5) ggarrange(plotlist=visplots, labels = joinbatch$stage) dev.off() visplots[[1]] #Model diagnostics par(mfrow=c(4,4)) lapply(models$glms,function(x)gam.check(x) ) ##---------------------------------------------------------------------------- ##Swimming behaviour ##---------------------------------------------------------------------------- alltracks<-read_excel(path,sheet=4) alltracks<-as_tibble(alltracks) alltracks %>% filter(TRACK_DISPLACEMENT>2)%>% group_by(daysfspawn,Temp) %>% summarize(meansp=mean(TRACK_MEAN_SPEED,na.rm=T),sdt=sd(TRACK_MEAN_SPEED,na.rm=T)) ##Swimming speed-statistical analysis thelmfun2<-function(y,mydat){ mod0<-lm(as.formula(paste(y,"~1")),data=mydat) mod1<-lm(as.formula(paste(y,"~","daysfspawn")),data=mydat) mod2<-lm(as.formula(paste(y,"~","daysfspawn","+dir")),data=mydat) mod3<-lm(as.formula(paste(y,"~","daysfspawn","*dir")),data=mydat) myan=anova(mod0,mod1,mod2,mod3,test="Chisq") aicvals=AIC(mod0,mod1,mod2,mod3) bestmod<-if(myan$`Pr(>Chi)`[2]<0.05){mod1}else{mod0} bestsum=summary(mod1) bestplots<-plot(bestmod) return(list(myan,aicvals,bestsum,bestplots)) } nested_dat<-alltracks %>% mutate(daysfspawn=as_factor(daysfspawn))%>% arrange(Temp)%>% group_by(Temp)%>% nest() map(nested_dat$data,~thelmfun2(y="TRACK_MEAN_SPEED",mydat=.x)) nested_dat<-nested_dat%>% mutate(mods=map(data,~lm(TRACK_MEAN_SPEED~daysfspawn,data=.x))) #Plot model effects visswim<-lapply(nested_dat$mods,function(x) visreg(x,xvar="daysfspawn",plot=FALSE)) visswimplots<-lapply(visswim,function(x) ggplot() + geom_jitter(data=x$res,aes(daysfspawn, visregRes))+ geom_boxplot(data=x$fit,aes(y = visregFit,x=as.factor(daysfspawn)), colour="black")+ geom_errorbar(data=x$fit,aes(ymin=visregLwr,ymax=visregUpr,x=daysfspawn))+ theme_light()+ theme(legend.title=element_blank(),legend.position=c(0.9,0.4), legend.background = element_blank())+ ylab(expression(paste("Swimming velocity (mm/ ", s^1,")",sep="")))+ xlab("Days from spawning")) plotnames2<-as.list(nested_dat$Temp) pdf(file="Graphs/FigureS3.pdf", width=4.5,height=4.5) grid.arrange(grobs=mapply(X=visswimplots, Y=plotnames2,function(X,Y){arrangeGrob(X,top=Y)})) dev.off() ##Compare between temperatures alltracks<-alltracks%>% mutate(Temp=as_factor(Temp),daysfspawn=as_factor(daysfspawn)) thelmfun3<-function(y,mydat){ mod0<-lm(as.formula(paste(y,"~1")),data=mydat) mod1<-lm(as.formula(paste(y,"~","daysfspawn")),data=mydat) mod2<-lm(as.formula(paste(y,"~","daysfspawn","+Temp")),data=mydat) mod3<-lm(as.formula(paste(y,"~","daysfspawn","*Temp")),data=mydat) myan=anova(mod0,mod1,mod2,mod3,test="Chisq") aicvals=AIC(mod0,mod1,mod2,mod3) bestmod<-if(myan$`Pr(>Chi)`[2]<0.05){mod1}else{mod0} bestplots<-plot(bestmod) return(list(myan,aicvals,bestplots)) } thelmfun3(mydat=alltracks,y="TRACK_MEAN_SPEED") bestmodel<-lm(TRACK_MEAN_SPEED~ daysfspawn+Temp, data=alltracks) pdf(file="Graphs/FigureS4.pdf", width=4.5,height=4) visreg(bestmodel,xvar="daysfspawn",by="Temp",overlay=T,line=list(col=c("gray","black")), fill=list(col=c("gray","black"),a=0.2),points=list(col=c("gray","black")), ylab=expression(paste("Swimming velocity (mm/ ", s^1,")",sep="")), xlab="Days from spawning") dev.off() ##Percentage up/down-statistical analysis direc<-alltracks%>% group_by(daysfspawn,Batch,Temp,dir)%>% summarize(ndir=n())%>% mutate(per= prop.table(ndir) * 100) direc%>%group_by(daysfspawn,Temp,dir)%>% summarize(meansp=mean(per,na.rm=T),sdt=sd(per,na.rm=T)) direc%>%group_by(Temp,dir)%>% summarize(meansp=mean(per,na.rm=T),sdt=sd(per,na.rm=T)) nested_direc<-direc%>% group_by(Temp)%>% nest() map(nested_direc$data,~thelmfun2(y="per",mydat=.x))