###### This is the complete code for all talbes, figures, and analyses presented in the paper ### "Sexual Signals Reflect Telomere Dynamics in a Wild Bird" ### by Conor C Taff, and Corey R Freeman-Gallant, Ecology & Evolution 2017. ### Code was last run in R v. 3.3.3 for Mac OS X 10.11.6. #General data description: These data were collected from a free living population of common yellowthroat #warblers breeding at two field sites near Saratoga Springs, NY, USA from May-July 2008-2012. Please #see the accompanying paper (Taff & Freeman-Gallant 2014, Hormones & Behavior) for more detail on the #study system, data collection, and sample processing. Additional information, data, and analyses of #this bird population can be found in the following papers: #(1) Mitchell et al. 2007. Animal Behaviour. 73: 165-170. #(2) Dunn et al. 2008. Journal of Avian Biology. 39: 66-72. #(3) Dunn et al. 2010. Functional Ecology. 24: 159-158. #(4) Freeman-Gallant et al. 2010. Evolution. 64: 1007-1017. #(5) Freeman-Gallant et al. 2011. Biology Letters. 7: 429-432. #(6) Taff et al. 2011. Animal Behaviour. 81: 619-625. #(7) Taff et al. 2012. Wilson Journal of Ornithology. 124: 370-374. #(8) Taff et al. 2012. Animal Behaviour. 84: 813-821. #(9) Taff et al. 2013. Journal of Evolutionary Biology. 26: 1392-1405. #(10) Freeman-Gallant et al. 2014. Journal of Evolutionary Biology. 27: 982-991. #(11) Taff et al. 2014. Hormones & Behavior. 66: 276-282. #(12) Taff et al. 2014. Proceedings of the Royal Society of London, B. 281: 20141974. #(13) Whittingham et al. 2015. Molecular Ecology. 7: 1584-1595. #(14) Taff & Freeman-Gallant. 2016. Ethology. 122: 319-328. #NOTE: In all cases, blank values in data indicate unknown or missing data. Why the data are missing varies by #case (e.g. some are from lack of observation, insufficient blood sample, or failed assay); #thus, blank values should not be taken as meaningful in most cases. ## Before beginning, set the working directory to the folder with the data files. ## Load in some packages that contain functions that will be used for analysis later on. ## Should all be available on CRAN; newer versions may necessitate code changes below. library(plotrix) library(rethinking) library(lme4) library(bbmle) library(LMERConvenienceFunctions) library(lmtest) library(colorspace) library(MASS) ## Load in data and make some new columns. ## 'DataBySeason' is organized with each male breeding season as one row. This data table ## includes many males that have only one year of data, or non-consecutive years. ## Those males are not included in most analyses, but are used in a few cases. data1<-read.csv("DataBySeason.csv") ## 'DataByYear' is used for cross year comparisons. Each row represents one pair of years with ## all relevant measures in both year n and n+1 plus some measures of changes between ## years. For some males a third consecutive year is also included in a row to make ## plotting easier, but that third observation also appears as year n+1 in a separate ## row (e.g. male observed in three years would have two rows, one has 2008-2009 + 2010 ## data, the other has 2009-2010 data only). data2<-read.csv("DataByYear.csv") ## Arcsine transform OXS data in 'data1' and 'data2' data2$as_AvgDNAnn1<-asin(sqrt(data2$AvgDNAnn1*.01))*100 data1$as_OXS<-asin(sqrt(data1$OXS*.01))*100 ## Set size of posterior sampling for confidence interval calculation posterior.size<-1000000 ## Plot into R device or save to pdf? plot.to.dev<-0 # Set as 0 (plot on screen) or 1 (plot to pdf) ## TABLE 1. This model includes average TAC, average oxidative damage, ## year n telo rq, and a random effect of male id. ## The data are limited to those in which all of the predictor ## variables were measured, otherwise the likelihood ratio tests for the full ## vs. reduced models do not work because the sample sizes and data included differ. ## This is a particular problem for TAC, as the reduced model that does not include ## TAC would have a larger sample size (and thus show TAC to be ns). data3<-subset(data2,data2$All_Data=="Yes") m1<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+n_Bib+n_Mask+Experience2+ (1|Male),data=data3) m1.telo<-lmer(TeloDStar~as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+n_Bib+n_Mask+Experience2+ (1|Male),data=data3) m1.DNA<-lmer(TeloDStar~n_ModTelo+AvgTACnn1+n_UVB+n_YelB+n_Ccar+n_Bib+n_Mask+Experience2+ (1|Male),data=data3) m1.TAC<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+n_UVB+n_YelB+n_Ccar+n_Bib+n_Mask+Experience2+ (1|Male),data=data3) m1.UVB<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_YelB+n_Ccar+n_Bib+n_Mask+Experience2+ (1|Male),data=data3) m1.YelB<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_Ccar+n_Bib+n_Mask+Experience2+ (1|Male),data=data3) m1.Ccar<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Bib+n_Mask+Experience2+ (1|Male),data=data3) m1.Bib<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+n_Mask+Experience2+ (1|Male),data=data3) m1.Mask<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+n_Bib+Experience2+ (1|Male),data=data3) m1.Age<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+n_Bib+n_Mask+ (1|Male),data=data3) anova(m1,m1.telo) anova(m1,m1.DNA) anova(m1,m1.TAC) anova(m1,m1.UVB) anova(m1,m1.YelB) anova(m1,m1.Ccar) anova(m1,m1.Bib) anova(m1,m1.Mask) anova(m1,m1.Age) ## Look for two way interactions m1.Age.TAC<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+ n_Bib+n_Mask+Experience2+Experience2*AvgTACnn1+(1|Male),data=data3) m1.Age.UVB<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+ n_Bib+n_Mask+Experience2+Experience2*n_UVB+(1|Male),data=data3) m1.Age.YelB<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+ n_Bib+n_Mask+Experience2+Experience2*n_YelB+(1|Male),data=data3) m1.Age.Ccar<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+ n_Bib+n_Mask+Experience2+Experience2*n_Ccar+(1|Male),data=data3) m1.Age.Bib<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+ n_Bib+n_Mask+Experience2+Experience2*n_Bib+(1|Male),data=data3) m1.Age.Mask<-lmer(TeloDStar~n_ModTelo+as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+ n_Bib+n_Mask+Experience2+Experience2*n_Mask+(1|Male),data=data3) anova(m1.Age.TAC,m1) anova(m1.Age.UVB,m1) anova(m1.Age.YelB,m1) anova(m1.Age.Ccar,m1) anova(m1.Age.Bib,m1) anova(m1.Age.Mask,m1) ## Proceed using m1.Age.TAC as the model for later predictions ## Figure 1a. Photo of common yellowthroat added in photoshop ## Figure 1b. Comparison of year n uv brightness to telomere length in years n, n+1, and n+2. ## This includes only a group of 47 males that were sampled as 1 year olds initially. Of those ## males, 27 were also sampled as two year olds, and 9 were sampled as three year olds. ## Conditionally turn on plotting to pdf file if(plot.to.dev==1){ pdf("Figure1.pdf",width=10,height=5.3) } ## Set plotting parameters par(mfrow=c(1,2)) data.a<-subset(data1,data1$Age=="FY") data.b<-subset(data2,data2$Age=="FY-SY") plot(data.a$MODTELRQ~data.a$syUVB,ylim=c(0.25,1.75), pch=18,col=col.alpha(rgb(230/255,159/255,0/255),alpha=0.7),ylab="Telomere RQ", xlab="Year N UV Brightness",xlim=c(-2.5,2.5)) mod.a<-lm(data.a$MODTELRQ~data.a$syUVB) abline(mod.a,col=rgb(230/255,159/255,0/255),lwd=2) cor.test(data.a$MODTELRQ,data.a$UVB) points(data.b$n_UVB,data.b$n1_ModTelo,pch=17,col="slateblue") mod.b<-lm(data.b$n1_ModTelo~data.b$n_UVB) abline(mod.b,col="slateblue",lwd=2) cor.test(data.b$n1_ModTelo,data.b$n_UVB) points(data.b$n_UVB,data.b$n2_ModTELO,pch=21,col="black") mod.c<-lm(data.b$n2_ModTELO~data.b$n_UVB) abline(mod.c,col="black",lwd=2) cor.test(data.b$n2_ModTELO,data.b$n_UVB) text(0,.45,"Year N, P = 0.06",pos=4,col=rgb(230/255,159/255,0)) text(0,.35,"Year N+1, P = 0.003",col="slateblue",pos=4) text(0,.25,"Year N+2, P = 0.31",col="black",pos=4) corner.label("B") ## Figure 1c. Comparison of year n coloration to change in telomere length between years. ## Delta telomere values are corrected for regression to the mean. ## This is equivalent to the relationship predicted in the ## multivariate model in Table 1. The plotted relationship is derived directly from the model ## and therefore also controls for the other predictors included in that model. ## All predictors other than year n color are held at their mean values. ## Panel 1. Year n uvb vs. change in telomeres. plot(data2$n_UVB,data2$TeloDStar,pch=16,col=col.alpha("black",alpha=0.5), xlab="Year N UV Brightness",ylab=expression(paste(Delta," TRQ Year N to N+1")) ,xlim=c(-2.5,2.5)) text(2,-0.46,"P = 0.002") corner.label("C") ### Add in predicted fit and CI from model 1 in Table 1. m1.Age.TAC2<-lmer(TeloDStar~as_AvgDNAnn1+AvgTACnn1+n_UVB+n_YelB+n_Ccar+ n_Bib+n_Mask+Experience2+Experience2*AvgTACnn1+n_ModTelo+(1|Male),data=data3) post.m1<-mvrnorm(n=posterior.size,mu=fixef(m1.Age.TAC2),Sigma=vcov(m1.Age.TAC2)) x.post<-c(mean(data3$n_ModTelo),mean(data3$as_AvgDNAnn1),mean(data3$AvgTACnn1),mean(data3$n_YelB),mean(data3$n_Ccar), mean(data3$n_Bib),mean(data3$n_Mask),mean(data3$Experience2)) r<-seq(from=-3.5,to=3.5,by=0.1) x.PC1<-sapply(r,function(z)mean(post.m1[,1]+ post.m1[,2]*x.post[2]+post.m1[,3]*x.post[3]+post.m1[,4]*z+post.m1[,5]*x.post[4]+ post.m1[,6]*x.post[5]+post.m1[,7]*x.post[6]+post.m1[,8]*x.post[7]+post.m1[,9]*x.post[8]+ post.m1[,10]*x.post[1]+post.m1[,11]*x.post[3]*x.post[8])) ci.PC1<-sapply(r,function(z)HPDI(post.m1[,1]+ post.m1[,2]*x.post[2]+post.m1[,3]*x.post[3]+post.m1[,4]*z+post.m1[,5]*x.post[4]+ post.m1[,6]*x.post[5]+post.m1[,7]*x.post[6]+post.m1[,8]*x.post[7]+post.m1[,9]*x.post[8]+ post.m1[,10]*x.post[1]+post.m1[,11]*x.post[3]*x.post[8],prob=0.95)) par(new=TRUE) lines(r,x.PC1,lwd=2,col="black") shade(ci.PC1,r,col=col.alpha("slateblue",0.4)) ## Conditionall turn off plotting device. if(plot.to.dev==1){dev.off()} ## Other analyses, statistics calculated and reported in text of manuscript, but not ## included in a figure or table. ## Correlations between size PC1 and telomere length in years N, N+1, and N+2 cor.test(data.a$MODTELRQ,data.a$Bib_Size) cor.test(data.b$n1_ModTelo,data.b$n_Bib) cor.test(data.b$n2_ModTELO,data.b$n_Bib) cor.test(data.a$MODTELRQ,data.a$Mask_Size) cor.test(data.b$n1_ModTelo,data.b$n_Mask) cor.test(data.b$n2_ModTELO,data.b$n_Mask) ## Tests for effect of TRQ on survival a few different ways ## Test for relationship between survival and telomere for all males pooled: ## RQ in year n: data3<-subset(data1,data1$Survive_Next=="Yes"|data1$Survive_Next=="No") m.s<-glmer(Survive_Next~MODTELRQ+(1|Male),family=binomial,data=data3) m.s.r<-glmer(Survive_Next~(1|Male),family=binomial,data=data3) anova(m.s,m.s.r) ## delta RQ to survival in year n+2 data4<-subset(data2,data2$SurviveNext=="Yes"|data2$SurviveNext=="No") m.s2<-glmer(SurviveNext~TeloDStar+(1|Male),family=binomial,data=data4) m.s2.r<-glmer(SurviveNext~(1|Male),family=binomial,data=data4) anova(m.s2,m.s2.r) ## Tests for shortening of telomeres between years pooled or within males ## Within male change in telo rq year n to n+1 t.test(data2$n_ModTelo,data2$n1_ModTelo,paired=TRUE) ## Differences in TRQ between years pooling age groups m<-lmer(MODTELRQ~Year_In_Site+(1|Male),data=data1) m.r<-lmer(MODTELRQ~(1|Male),data=data1) anova(m,m.r) ### Figure 2. TAC vs. telomere erosion for young vs. old males ## Conditionally turn on plotting to pdf file if(plot.to.dev==1){ pdf("Figure2.pdf",width=5.2,height=5.2) } ## Plot the main figure data3<-subset(data2,data2$All_Data=="Yes") data3$Experience2<-as.numeric(data3$Age) data3$PointColor<-ifelse(data3$Experience2==1,"gray75","slateblue") data3$PointShape<-ifelse(data3$Experience2==1,17,16) plot(data3$AvgTACnn1,data3$TeloDStar,pch=data3$PointShape,col=data3$PointColor, xlim=c(0.5,3),ylim=c(-.5,.5), ylab=expression(paste(Delta," TRQ Year N to N+1")), xlab="Average TAC") ### Add in predicted fit and CI from model 1 in Table 1. Note that the model here is ## modified slightly from Table 1 because the response variable includes a correction ## for year N telomere RQ rather than including year N RQ as a predictor. In the table ## year N RQ is retained as a predictor so that the effect can be seen more easily. ## First plot for males transitioning from year 1 to year 2 of life r<-seq(from=0,to=3.5,by=0.1) x.post<-c(mean(data3$as_AvgDNAnn1),mean(data3$n_UVB), mean(data3$n_YelB),mean(data3$n_Ccar),mean(data3$n_Bib),mean(data3$n_Mask), mean(data3$Experience2),mean(data3$n_ModTelo)) x.PC1<-sapply(r,function(z)mean(post.m1[,1]+post.m1[,2]*x.post[1]+ post.m1[,3]*z+post.m1[,4]*x.post[2]+post.m1[,5]*x.post[3]+ post.m1[,6]*x.post[4]+post.m1[,7]*x.post[5]+post.m1[,8]*x.post[6]+ post.m1[,9]*1+post.m1[,10]*x.post[8]+post.m1[,11]*z*1)) ci.PC1<-sapply(r,function(z)HPDI(post.m1[,1]+post.m1[,2]*x.post[1]+ post.m1[,3]*z+post.m1[,4]*x.post[2]+post.m1[,5]*x.post[3]+ post.m1[,6]*x.post[4]+post.m1[,7]*x.post[5]+post.m1[,8]*x.post[6]+ post.m1[,9]*1+post.m1[,10]*x.post[8]+post.m1[,11]*z*1,prob=0.95)) lines(r,x.PC1,lwd=2,col="black") shade(ci.PC1,r,col=col.alpha("gray80",0.4)) ## Now plot for males transitioning from later years x.PC1<-sapply(r,function(z)mean(post.m1[,1]+post.m1[,2]*x.post[1]+ post.m1[,3]*z+post.m1[,4]*x.post[2]+post.m1[,5]*x.post[3]+ post.m1[,6]*x.post[4]+post.m1[,7]*x.post[5]+post.m1[,8]*x.post[6]+ post.m1[,9]*2+post.m1[,10]*x.post[8]+post.m1[,11]*z*2)) ci.PC1<-sapply(r,function(z)HPDI(post.m1[,1]+post.m1[,2]*x.post[1]+ post.m1[,3]*z+post.m1[,4]*x.post[2]+post.m1[,5]*x.post[3]+ post.m1[,6]*x.post[4]+post.m1[,7]*x.post[5]+post.m1[,8]*x.post[6]+ post.m1[,9]*2+post.m1[,10]*x.post[8]+post.m1[,11]*z*2,prob=0.95)) lines(r,x.PC1,lwd=2,col="slateblue") shade(ci.PC1,r,col=col.alpha("slateblue",0.4)) ## Conditionall turn off plotting device. if(plot.to.dev==1){dev.off()}