--- title: "Elevated pCO~2~ affects tissue biomass composition, but not calcification, in a reef coral under two light regimes" author: "CB Wall, Hawai'i Institute of Marine Biology. Wall et al. (2017) Royal Society Open Science" date: "9/29/2017" output: html_document: code_folding: hide toc: yes toc_depth: 4 toc_float: yes --- ## Experimental design - Treatments were created in 24 tanks of 6 replicate treatents: - Low Light-Ambient pCO~2~ (LL-ACO~2~), Low light-High pCO~2~ (LL-HCO~2~), High light-Ambient pCO~2~ (HL-ACO~2~), and High light-High pCO~2~ (HL-HCO~2~). - two tanks malfunctioned and the the final tank replicates were 6 tanks for all treatments, but 4 replicate tanks for HL-HC treatment - Experimental factors (pCO~2~ and light) were treated as fixed effects. Tank replicates were treated as random effects nested within treatment. Coral colony was treated as a random effect. ## Seawater chemistry Exploring seawater chemistry data, generating data summaries to be used in tables and figures. Statistical tests performed to test for differences in pH, total alkalinity, and pCO~2~ among all treatment tanks and among replicate treatment tanks within a CO~2~ treatment (ambient or high). ```{r, eval=FALSE, include=FALSE} # Wherever you have the .Rmd file, R studio will force your working directory to be in this location, this can be changed (see 'setup' chunk in Rmarkdown help). # I made a file names 'markdown' that holds my .Rmd file. Within the 'markdown' folder there are two folders: 1 for data, another for figures. ``` Attach packages ```{r, message=FALSE, warning=FALSE, results="hide"} ## load packages library(lmerTest) library(lmtest) library(lme4) library(nlme) library(car) library(effects) library(gplots) library(plotrix) library(psych) library(ggplot2) library(grid) library(gridExtra) library(scales) library(MASS) library(lsmeans) library(pbkrtest) library(knitr) ``` ```{r, include=FALSE, results="hide"} ############################################################################################ #### Temperature data for Table 1--seawater chem all Treatments: mean, SE, N conditions #### ############################################################################################ rm(list=ls()) setwd('markdown/data') # Temperature data from treatment tanks temp.data<-read.csv("Exp_temp salinity_OAxLight.csv") temp.data$Treatment<-factor(temp.data$Treatment, levels=c("LL-AC", "LL-HC", "HL-AC", "HL-HC")) mean.trt<-aggregate(cbind(Temperature)~Treatment, temp.data, mean, NA.rm=TRUE) n.trt<-aggregate(cbind(Temperature)~Treatment, temp.data, length) se.trt<-aggregate(cbind(Temperature)~Treatment, temp.data, std.error) #TEMP df<-data.frame(mean.trt, se.trt[c(0, 2)], n.trt[c(0,2)]) colnames(df) <- c("Treatment", "Temperature", "Temp.SE", "Temp.n") kable(df, digits =2, "html") mean(df$Temp, na.rm=TRUE) std.error(df$Temp, na.rm=TRUE) sum(df$Temp.n) ################ #SALINITY mean.trt<-aggregate(cbind(Salinity)~Treatment, temp.data, mean, NA.rm=TRUE) n.trt<-aggregate(cbind(Salinity)~Treatment, temp.data, length) se.trt<-aggregate(cbind(Salinity)~Treatment, temp.data, std.error) df<-data.frame(mean.trt, se.trt[c(0, 2)], n.trt[c(0,2)]) colnames(df) <- c("Treatment", "Salinity", "Sal.SE", "Sal.n") kable(df, digits =2, "html") mean(df$Salinity, na.rm=TRUE) std.error(df$Salinity, na.rm=TRUE) sum(df$Sal.n) ``` Attach the first dataset for seawater chemistry. ```{r, message=FALSE, warning=FALSE} ################################################## ## EXPERIMENT PERIOD: 16 Dec 2014 - 14 Jan 2015 ## ################################################## rm(list=ls()) setwd('markdown/data') # attach the data with seawater chemistry information SWdata<-read.csv("SWchem_Treatments_OAxLight.csv", header=T, na.string=NA) SWdata$Treatment<-factor(SWdata$Treatment, levels=c("LL-AC", "LL-HC", "HL-AC", "HL-HC")) ###################################### ## Tank means ## ###################################### mean.summary<-aggregate(cbind(pCO2, pHT)~Tank, SWdata, mean) n.summary<-aggregate(cbind(pCO2, pHT)~Tank, SWdata, length) se.summary<-aggregate(cbind(pCO2, pHT)~Tank, SWdata, std.error) df<-data.frame(mean.summary, se.summary[c(0, 2,3)], n.summary[c(0,2)]) colnames(df) <- c("Tank", "pCO2", "pHT", "CO2-SE", "pH-SE", "n") ######################################################### #### across all tanks: mean, SE, N of CO2 conditions #### ######################################################### mean.summary<-aggregate(cbind(pCO2, pHT)~CO2_Treat, SWdata, mean) n.summary<-aggregate(cbind(pCO2, pHT)~CO2_Treat, SWdata, length) se.summary<-aggregate(cbind(pCO2, pHT)~CO2_Treat, SWdata, std.error) df<-data.frame(mean.summary, se.summary[c(0, 2,3)], n.summary[c(0,2)]) colnames(df) <- c("CO2_Treat", "pCO2", "pHT", "CO2-SE", "pH-SE", "n") ``` ```{r, include=FALSE, results="hide"} df kable(df, digits=2, "html") ###################################################################################### #### Table data for Table 1--seawater chem all Treatments: mean, SE, N conditions #### ###################################################################################### # note Temperature here reflects seawater measure at time of collection, not independent monitoring of treatment conditions over period of study (> frequency) SWdata$Treatment<-factor(SWdata$Treatment, levels=c("LL-AC", "LL-HC", "HL-AC", "HL-HC")) mean.trt<-aggregate(cbind(pHT, ALK, pCO2, HCO3, CO3, OmegaAragonite)~Treatment, SWdata, mean) n.trt<-aggregate(cbind(pHT, ALK, pCO2, HCO3, CO3, OmegaAragonite)~Treatment, SWdata, length) se.trt<-aggregate(cbind(pHT, ALK, pCO2, HCO3, CO3, OmegaAragonite)~Treatment, SWdata, std.error) df<-data.frame(mean.trt, se.trt[c(0, 2:7)], n.trt[c(0,7)]) colnames(df) <- c("Treatment", "pHT", "ALK", "pCO2", "HCO3", "CO3", "Omega", "pHT-SE", "ALK-SE", "pCO2-SE", "HCO3-SE", "CO3-SE", "Omega-SE", "n") kable(df, digits =2, "html") ``` #### pH and CO~2~ treatment figure Generate a figure of mean +/- SE pHT and pCO~2~ for each treatment tank. ```{r, fig.width=7, fig.height=2.5, message=FALSE, warning=FALSE, results="hide"} # set working directory to 'figures' for figure export setwd('markdown/figures') ######################################### ## EXPERIMENT: Figures of CO2 and pHT ## ######################################### colors=c("lightskyblue", "lightskyblue", "lightskyblue", "lightcoral","lightcoral", "lightcoral","lightskyblue", "lightskyblue", "lightskyblue", "lightcoral","lightcoral", "lightcoral","lightskyblue", "lightskyblue", "lightskyblue", "lightcoral","lightskyblue", "lightskyblue", "lightskyblue", "lightcoral","lightcoral", "lightcoral") # plot for all tank pH and pCO2 ################# ###### pH ###### ################# mean<-tapply(SWdata$pHT, SWdata$Tank, mean); se<-tapply(SWdata$pHT, SWdata$Tank, std.error); pH<-data.frame(mean, se); pH["Tank"] <- c("01", "02", "03", "04", "05","06","07","08","09","10","11","12","13","14","15","17","19","20","21","22","23","24"); colnames(pH) <- c("mean", "se","Tank") Fig.pH <- ggplot(pH, aes(x=Tank, y=mean)) + geom_bar(position="dodge", stat="identity", fill=colors) + geom_errorbar(aes(ymin=pH$mean - pH$se, ymax=pH$mean + pH$se), width=0.1) + xlab(expression("Experimental Tank")) + ylab(expression(pH[T])) + scale_y_continuous(limits=c(7.0,8.4),oob = rescale_none) Fig.pH # gives summary figure of tank pH mean +/-SE across all tanks dev.copy(pdf, "pHbytank.pdf", height=3, width=8) dev.off() ################# ###### CO2 ###### ################# mean<-tapply(SWdata$pCO2, SWdata$Tank, mean); se<-tapply(SWdata$pCO2, SWdata$Tank, std.error); pCO2<-data.frame(mean, se); pCO2["Tank"] <- c("01", "02", "03", "04", "05","06","07","08","09","10","11","12","13","14","15","17","19","20","21","22","23","24"); colnames(pCO2) <- c("mean", "se","Tank") Fig.pCO2 <- ggplot(pCO2, aes(x=Tank, y=mean)) + geom_bar(position="dodge", stat="identity", fill=colors) + geom_errorbar(aes(ymin=pCO2$mean - pCO2$se, ymax=pCO2$mean + pCO2$se), width=0.1) + xlab(expression("Experimental Tank")) + ylab(expression(mu*atm~pCO[2])) + scale_y_continuous(limits=c(0,1200),oob = rescale_none) Fig.pCO2 # gives summary figure of tank pCO2 mean +/-SE across all tanks ##### save the figure and export to directory? #### dev.copy(pdf, "CO2bytank.pdf", height=3, width=8) dev.off() ``` #### Testing treatment conditions - Statistical tests of tank effects, examining Total Alkalinity (*A~T~* or ALK), CO~2~, and pH~T~. - First, testing for treatment effects, ie., are there two CO~2~ treatments and is *A~T~* stable across all tanks? ```{r, ECHO=FALSE} #################################################################### # EXPERIMENT PERIOD: test CO2 TREATMENTS, either Ambient or High ## #################################################################### TA<-lm(ALK~CO2_Treat, data=SWdata); TAmod<-anova(TA) # no difference in TA across all CO2-TRT P=0.1101 CO2<-lm(pCO2~CO2_Treat, data=SWdata); CO2mod<-anova(CO2) # CO2 treatments differ, as expected pHT<-lm(pHT~CO2_Treat, data=SWdata); pHTmod<-anova(pHT) # pH differs, as expected ``` **Total Alkalinity** across all tanks ```{r, echo=FALSE} TAmod ``` **pCO~2~** across all tanks ```{r, echo=FALSE} CO2mod ``` **pH~T~** across all tanks ```{r, echo=FALSE} pHTmod ``` - There are differences in treatments (ACO~2~ and HCO~2~) and total alkalinity does not differ among treatments. - NEXT, testing whether replicate tanks are different, based on CO~2~ treatments (Ambient CO~2~ vs. High CO~2~). To do this, divide dataframe into High and Ambient pCO~2~ treatments and examine replicate tanks. #### Testing replicate tanks *FIRST, testing for differences among replicate Ambient pCO~2~ tanks...* ```{r} ######### examine differences among replicate pCO2-treatments ACO2<-SWdata[(SWdata$CO2_Treat=="ACO2"),] # only Ambient CO2 tanks HCO2<-SWdata[(SWdata$CO2_Treat=="HCO2"),] # only High CO2 tanks ########## ACO2 tanks ############# ################################### TA<-lm(ALK~Tank, data=ACO2); TAmod<-anova(TA) # no difference in TA across all CO2-TRT P=0.1101 CO2<-lm(pCO2~Tank, data=ACO2); CO2mod<-anova(CO2) posthoc<-lsmeans(CO2, pairwise~Tank, adjust="Tukey") contrasts<-cld(posthoc, Letters=letters) # CO2 do not differ from each other, but still some variance (P=0.0601) pHT<-lm(pHT~Tank, data=ACO2); pHTmod<-anova(pHT) posthoc<-lsmeans(pHT, pairwise~Tank, adjust="Tukey") contrasts<-cld(posthoc, Letters=letters) # pH do not differ from each other, but still some variance (P=0.0840) ``` **Total Alkalinity** across ACO~2~ tanks ```{r, echo=FALSE} TAmod ``` **pCO~2~** across ACO~2~ tanks ```{r, echo=FALSE} CO2mod ``` **pHT** across ACO~2~ tanks ```{r, echo=FALSE} pHTmod ``` *...THEN, testing for differences among replicate High pCO~2~ tanks...* ```{r} ########## HCO2 tanks ############# ################################### TA<-lm(ALK~Tank, data=HCO2); TAmod<-anova(TA) # no difference in TA across all CO2-TRT (P=0.1101) CO2<-lm(pCO2~Tank, data=HCO2); CO2mod<-anova(CO2) # CO2 do not differ from each other (P = 0.9739) pHT<-lm(pHT~Tank, data=HCO2); pHTmod<-anova(pHT) # pH do not differ from each other (P = 0.9458) ######## note: greatest variance in pH and pCO2 in ACO2 not HCO2 tanks ``` **Total Alkalinity** across HCO~2~ tanks ```{r, echo=FALSE} TAmod ``` **pCO~2~** across HCO~2~ tanks ```{r, echo=FALSE} CO2mod ``` **pHT** across HCO~2~ tanks ```{r, echo=FALSE} pHTmod ``` ## Physiological Responses - Response variables are normalized to area (cm^2^) or grams of ash-free dry weight biomass (gdw), or both. - Calcification, protein, carbohydrates, lipids, and tissue energy content are normalized to area and biomass. - Total biomass, *Symbiodinium* density, chlorophyll *a* and *c~2~* is normalized to area. - Photopigments (chlorophyll *a* and *c~2~*) are also normalized to symbiont cell densities. First, normalize data found in raw physiology data file **(code details below)** ```{r, message=FALSE, warning=FALSE} rm(list=ls()) setwd('markdown/data') data<-read.csv("Wall et al_OAxLight_R.csv", header=T, na.string=NA) ################################################ # calculate, normalized dependent variables ################################################ # helpful shorthand SA<-data$surface.area.cm2 # surface area in cm2 blastate<-data$total.blastate.ml # tissue slurry blastate in ml # CaCO3 area: calcif.mg.CaCO3..cm2..d == convert change in CaCO3 g to mg, dive by SA and # days between measurements = 23 d data$calcif.cm2<- (data$change.in.weight.g*1000)/SA/23 # CaCO3 biomass: mg.CaCO3..gdw..d == convert change in CaCO3 g to mg, dive by biomass and # days between measurements = 23 d data$calcif.gdw<-((data$change.in.weight.g*1000)/(data$AFDW.g.ml*blastate))/23 # AFDW.mg..cm2 == convert AFDW g to mg, mutiply by blastate volume, divide by cm2 data$biomass<- (data$AFDW.g.ml*1000*blastate)/SA # Symbiodinium.cells..cm2 == cell.ml * blastate / SA data$zoox<- (data$cells.ml*blastate)/SA # ug.chlorophyll.a..cm2 == ug.chl.a.ml * blastate / SA data$chla<- (data$ug.chl.a.ml*blastate)/SA # ug.chlorophyll.c2..cm2 == ug.chl.c2.ml * blastate / SA data$chlc2<- (data$ug.chl.c2.ml*blastate)/SA # pg.chlorophyll.a..cell == ug.chl.a.ml * 10^6 / cells.ml data$chlacell<- (data$ug.chl.a.ml*10^6)/data$cells.ml # pg.chlorophyll.c2..cell == ug.chl.c2.ml * 10^6 / cells.ml data$chlc2cell<- (data$ug.chl.c2.ml*10^6)/data$cells.ml # mg.protein..cm2 == mg.protein.ml * blastate / SA data$prot.cm2<- (data$mg.protein.ml*blastate)/SA # g.protein..gdw == mg.protein.ml/1000 * blastate / SA data$prot.gdw<- (data$mg.protein.ml/1000)/data$AFDW.g.ml # mg.carbs..cm2 == mg.carbs.ml * blastate / SA data$carb.cm2<- (data$mg.carb.ml*blastate)/SA # g.carbs..gdw == mg.carbs.ml/1000 * blastate / SA data$carb.gdw<- (data$mg.carb.ml/1000)/data$AFDW.g.ml # mg.lipids..cm2 == mg.lipids.ml * blastate / SA data$lipid.cm2<- (data$mg.lipid.ml*blastate)/SA # g.lipids..gdw == in this case, this is the data directly from lipid extractions of 3 ml blastate data$lipid.gdw<- data$g.lipid.gdw # kJ.energy.content..cm2 == first, convert g/gdw to mg/gdw, then multiply coefficents for kJ/g macromolecule and sum together. Divide by 1000 to get kJ # -23.9 kJ/g protein, -17.5 kJ/g carbs, -39.5 kJ/g lipids kJ.gdw<-((data$prot.gdw)*23.9 + (data$carb.gdw)*17.5 + (data$lipid.gdw)*39.5) data$energy.cm2<-(kJ.gdw)*(data$AFDW.g.ml*blastate)/SA # kJ/gdw * (biomass in blastate)/SA data$energy.gdw<-kJ.gdw df<-data[, c(1:8, 18:33)] ########################## #### Treatment levels #### ########################## CO2<- df$CO2 Light<- df$Light Colony<- df$Colony Treatment<- df$Treatment Tank<- df$Tank df$Tank<-as.factor(df$Tank) ##################################### ####### all standardized data ####### ##################################### ###### response variables ###### ####### data sobriquets ######## df$calcif.cm2[31]<-NA # remove outlier calcif.cm2<-df$calcif.cm2 calcif.gdw<-df$calcif.gdw biomass<-df$biomass zoox<-df$zoox chla<-df$chla chlc2<-df$chlc2 chlacell<-df$chlacell chlc2cell<-df$chlc2cell df$prot.cm2[10]<- NA #remove outlier prot.cm2<-df$prot.cm2 prot.gdw<-df$prot.gdw carb.cm2<-df$carb.cm2 carb.gdw<-df$carb.gdw lipid.cm2<-df$lipid.cm2 lipid.gdw<-df$lipid.gdw energy.cm2<-df$energy.cm2 energy.gdw<-df$energy.gdw ``` ## Assumptions of ANOVA: Checking for normality and homoscedasticity, using a combination of Shapiro-Wilks Test, graphical inspection of residuals and Q-Q plots. If data did not meet assuptions of ANOVA, a Box-Cox test was used to determine appropriate transformations. The code here transforms variables and tests for normality and homoscedasticity of these variables in a dataframe called "df.tests". - Variables with "trans" in header name have been transformed to meet assumptions of ANOVA ```{r, message=FALSE, warning=FALSE} Y1<-sqrt(calcif.cm2) # calcif.cm2 Y2<-biomass^-0.5 # AFDW cm2 Y3<-sqrt(zoox) #zoox cm2 Y4<-chla # chla cm2 Y5<-(chlc2)^0.25 # chlc2 cm2 Y6<-(chlacell)^-0.25 #chla per cell Y7<-(chlc2cell)^-0.25 #chlc2 per cell Y8<-prot.cm2^0.25 # protein cm2 Y9<-(carb.cm2)^0.25 #carbs cm2 Y10<-sqrt(lipid.cm2) #lipids cm2 Y11<-energy.cm2^0.25 # Joules/cm2 ################################## ####### biomass normalized ####### ################################# Y12<-sqrt(calcif.gdw) #calcification gdw Y13<-prot.gdw #protein gdw Y14<-sqrt(carb.gdw) #carbs gdw Y15<-lipid.gdw #lipids gdw <<1 # ev.area/sum(ev.area) #gives proportional eignevalues newdat.area<-PC.area$x[,1:4] # 2 PCAs explain 59% of variance; usable in regression SVM plot(PC.area, type="lines", main="PC.area eigenvalues") #plots PCA, equivalent to ev.area ``` ### Correlation tests for PC-area This code performs correlation tests on the relationship between dependent variables (as their PC loadings) and the PC of interest with eignevalue >1. The tests are performed for PC1 and PC2. For loop will perform statistical tests of significant correlation between PC "x" and variable "y" and generate scatter plot with trend line to show +/- correlations. - FIRST code chunk will run **PC1-Area and response variables** correlation tests. ```{r, results="hide", fig.show="hide"} ###################### # PC1 correlation tests ###################### # print(cor.test(PC.area$x[,1], Y)) ## 'print' to see results of all correlation tests with associated figure # names(pca.area.df) for(i in 2:7){ Y<-pca.area.df[,i] print(cor.test(PC.area$x[,1], Y)) ## correlation test for variables and PC1 plot(PC.area$x[,1], Y, ylab = colnames(pca.area.df)[i], xlab="PC1.area", main="PC1 correlation tests") ## shows plot of PC1 against variable of interest (Y) modPC<-lm(Y~PC.area$x[,1]) abline(modPC, col="red") } ``` - SECOND, code chunk will run **PC2-Area and response variables** correlation tests. ```{r, results="hide",fig.show="hide"} ###################### # PC2 correlation tests ###################### for(i in 2:7){ Y<-pca.area.df[,i] print(cor.test(PC.area$x[,2], Y)) ## correlation test for variables and PC2 plot(PC.area$x[,2], Y, ylab = colnames(pca.area.df)[i], xlab="PC2.area", main="PC2 correlation tests") ## shows plot of PC2 against variable of interest (Y) modPC<-lm(Y~PC.area$x[,2]) abline(modPC, col="red") } ``` ### Linear models for PC-area The component loadings were imported back to a Tank-Treatment dataframe to allign with CO~2~, Light, Tanks. The principal components (PC1 and PC2) were examined for normality, homoscedasticity, and run through an identical mixed effect linear model as is later used for univariate tests. ```{r, fig.show="hide", results="hide"} # head(pca.area.df) # area normalized data # new PCAdata dataframe--with treatment conditions and PCA PCdata<-data.frame(newdat.area) PCdata$CO2<-df.area.outlier.removed$CO2 PCdata$Light<-df.area.outlier.removed$Light PCdata$Treatment<-df.area.outlier.removed$Treatment PCdata$Tank<-df.area.outlier.removed$Tank PCdata$Colony<-df.area.outlier.removed$Colony ### test of random effects and models for normality mod1.lme<-lmer(PC1~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=PCdata, REML=FALSE) rand1<-rand(mod1.lme) mod2.lme<-lmer(PC2~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=PCdata, REML=FALSE) rand2<-rand(mod2.lme) R1 <- resid(mod1.lme); R2 <- resid(mod2.lme) R1.norm <- shapiro.test(R1) R2.norm <- shapiro.test(R2) # specify models and residuals op<-par(mfrow = c(2,2), mar=c(5,4,1,2)) fig<-plot(mod1.lme, add.smooth = FALSE, which=1) QQ <- qqnorm(R1, main = "") QQline <- qqline(R1) hist(R1, xlab="Residuals", main ="") ``` **Results for PC1-area lmer:** no effect of treatments on PC1-area ```{r} ################## # final models ################## mod1.lme<-lmer(PC1~CO2*Light+(1|Treatment:Tank), data=PCdata, REML=FALSE) rand1<-rand(mod1.lme) PC1.area<-anova(mod1.lme, type=2) print(PC1.area) ``` **Results for PC2-area lmer:** no effect of treatment on PC2-area ```{r} mod2.lme<-lmer(PC2~CO2*Light+(1|Treatment:Tank), data=PCdata, REML=FALSE) rand2<-rand(mod2.lme) PC2.area<-anova(mod2.lme, type=2) print(PC2.area) ``` ### PCA biomass normalization This code will run the identical analysis used on PCA-area, but now with repsonses normalized to biomass (gdw). ```{r, fig.height=3, fig.width=6} # names(df) df.gdw<-df[, c(1:6,10,11,18,20,22,24)] ## gdw normalized data df.gdw<-na.omit(df.gdw) # NAs removed # df.gdw.outlier.removed<-df.gdw[c(-16, -69), ] #outliers pca.gdw.df<-(df.gdw[,c(3,7:12)]) # all data with only "Treatment" as category colnames(pca.gdw.df)<-c("Treatment", "calcification", "biomass", "protein","carbohydrates", "lipids", "energy") ############################## # PCA biomass (treatment) PC.gdw <- prcomp(pca.gdw.df[,-1], center = TRUE, scale= TRUE) # print(PC.gdw) PCA.gdw.summary<-summary(PC.gdw) # the SDEV^2 is the eignevalue ev.gdw<-PC.gdw$sdev^2 # PC1-2 ave eignevalues >1 plot(PC.gdw, type="lines", main="PC.gdw eigenvalues") #plots PCA # ev.gdw/sum(ev.gdw) #gives proportional eignevalues newdat.gdw<-PC.gdw$x[,1:3] # 2 PCAs explain 73% of variance; usable in regression SVM ``` #### Correlation tests for PC-biomass This code performs correlation tests on the relationship between dependent variables (as their PC loadings) and the PC of interest with eignevalue >1. The tests are performed for PC1 and PC2. For loop will perform statistical tests of significant correlation between PC "x" and variable "y" and generate scatter plot with trend line to show +/- correlations. - FIRST, code chunk will run **PC1-Biomass and response variables** correlation tests ```{r, results="hide", fig.show="hide"} ###################### # correlation tests ###################### # head(pca.gdw.df) #### PC1 #### for(i in 2:7){ Y<-pca.gdw.df[,i] print(cor.test(PC.gdw$x[,1], Y)) ## correlation test for variables and PC2 plot(PC.gdw$x[,1],Y, ylab = colnames(pca.gdw.df)[i], xlab="PC1.gdw", main="PC1 correlation tests") ## shows plot of PC1 against variable of interest (Y) modPC<-lm(Y~PC.gdw$x[,1]) abline(modPC, col="red") } ``` - SECOND, code chunk will run **PC2-Biomass and response variables** correlation tests ```{r, results="hide", fig.show="hide"} for(i in 2:7){ Y<-pca.gdw.df[,i] print(cor.test(PC.gdw$x[,2], Y)) ## correlation test for variables and PC2 plot(PC.gdw$x[,2],Y, ylab = colnames(pca.gdw.df)[i], xlab="PC2.gdw", main="PC2 correlation tests") ## shows plot of PC2 against variable of interest (Y) modPC<-lm(Y~PC.gdw$x[,2]) abline(modPC, col="red") } ``` #### Linear models for PC-biomass Mixed effect linear models for biomass-normalized PC1 and PC2 ```{r, results="hide", fig.show="hide"} # biomass normalized data # head(pca.gdw.df) # head(df.gdw) #new PCAdata dataframe--with treatment conditions and PCA PCdata<-data.frame(newdat.gdw) PCdata$CO2<-df.gdw$CO2 PCdata$Light<-df.gdw$Light PCdata$Treatment<-df.gdw$Treatment PCdata$Tank<-df.gdw$Tank PCdata$Colony<-df.gdw$Colony ## testing model assumptions mod1.lme<-lmer(PC1~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=PCdata, REML=FALSE); rand(mod1.lme) mod2.lme<-lmer(PC2~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=PCdata, REML=FALSE); rand(mod2.lme) R1 <- resid(mod1.lme);R2 <- resid(mod2.lme) R1.norm <- shapiro.test(R1) R2.norm <- shapiro.test(R2) op<-par(mfrow = c(2,2), mar=c(5,4,1,2)) fig<-plot(mod1.lme, add.smooth = FALSE, which=1) QQ <- qqnorm(R1, main = "") QQline <- qqline(R1) hist(R1, xlab="Residuals", main ="") ``` **PC1-biomass model results:** no effect of treatments on PC1-biomass ```{r} ################## #### final models # PC1 mod1<-lmer(PC1~CO2*Light+(1|Treatment:Tank), data=PCdata, REML=FALSE) PC1.gdw<-anova(mod1, type=2) PC1.gdw ``` **PC2-biomass model results:** effect of pCO~2~ on PC2-biomass ```{r} #PC2 mod2<-lmer(PC2~CO2*Light+(1|Treatment:Tank), data=PCdata, REML=FALSE) PC2.gdw<-anova(mod2, type=2) # CO2 effect # capture output for all the final models (by area, and by biomass) capture.output(PC1.area, PC2.area, PC1.gdw, PC2.gdw, file="PCA-anova.results.txt") PC2.gdw ``` ## Figures: PCA These figures show principle components 1 and 2 (PC1 and PC2) for area-normalized and biomass-normalized data matrix. See manuscript for details. Together PCA-area and PCA-biomass figures can be found in **"Figure 1a and 1b"**. ### PCA-area graph of PC1 and PC2 for area-normalized responses. Together explaining ~60% of variance ```{r, message=FALSE, warning=FALSE, results="hide"} ################## #### graph it #### ################## setwd('markdown/figures') ## PC1 and PC2 PC1.PC2.area.fig <- ggbiplot(PC.area, choices = 1:2, obs.scale = 1, var.scale = 1, groups= pca.area.df[,1], ellipse = TRUE, ellipse.prob = 0.90, circle = FALSE) + scale_color_discrete(name = '') + theme_bw() + scale_x_continuous(breaks=pretty_breaks(n=5))+ scale_y_continuous(limits=c(-4, 4), breaks=pretty_breaks(n=5))+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) + theme(legend.text=element_text(size=15)) + theme(panel.background = element_rect(colour = "black", size=1))+ theme(legend.key = element_blank())+ theme(legend.direction = 'horizontal', legend.position = 'top') +theme(aspect.ratio=0.7) print(PC1.PC2.area.fig) dev.copy(pdf, "PC1.PC2.area.fig.pdf") dev.off() ``` ### PCA-biomass graph of PC1 and PC2 for biomass-normalized responses. Together explaining ~72% of variance ```{r, message=FALSE, warning=FALSE, results="hide"} ################## #### graph it #### ################## ## PC1 and PC2 setwd('markdown/figures') PC1.PC2.gdw.fig<- ggbiplot(PC.gdw, choices=1:2, obs.scale = 1, var.scale = 1, groups= pca.gdw.df[,1], ellipse = TRUE, ellipse.prob = 0.90, circle = FALSE) + scale_color_discrete(name = '') + theme_bw() + scale_x_continuous(limits=c(-5, 4), breaks=pretty_breaks(n=5))+ scale_y_continuous(limits=c(-5, 4), breaks=pretty_breaks(n=5))+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) + theme(legend.text=element_text(size=15)) + theme(panel.background = element_rect(colour = "black", size=1))+ theme(legend.key = element_blank())+ theme(legend.direction = 'horizontal', legend.position = 'top') +theme(aspect.ratio=0.7) print(PC1.PC2.gdw.fig) dev.copy(pdf, "PC1.PC2.gdw.fig.pdf") dev.off() ``` ## Univariate Models: lmer - all variables here are transformed if necessary and run through 3 separate models. A full model [(1|Colony and (1|Treatment:Tank)], a colony dropped model, and a nested Treatment:Tank dropped model. The best model is chosen from AIC tests and loglikihood ratio tests (LRT). ```{r, include=FALSE} #Y1 = calcif cm2 #Y2 = AFDW cm2 #Y3 = zoox cm2 #Y4 = chla cm2 #Y5 = chlc2 cm2 #Y6 = chla per cell #Y7 = chlc2 per cell #Y8 = protein cm2 #Y9 = carbs cm2 #Y10 = lipids cm2 #Y11 = kJoules cm2 #Y12 = calcification gdw #Y13 = protein gdw #Y14 = carbs gdw #Y15 = lipids gdw #Y16 = kJoules gdw ``` ### Area-normalized models **calcification cm^-2^** ```{r} ## ***********************************## ## AREA NORMALIZED ##################### ## ***********************************## ############################################ ** no colony effect ## Calcification ########################### ** effect of tank ############################################ ** no treatment effect #full model mod1<-lmer(Y1~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) #colony dropped mod2<-lmer(Y1~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) #tanks dropped mod3<-lmer(Y1~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand1<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is with colony dropped (=mod2) # tanks signficant, no effect of colony, no effect fixed effect final.mod<-lmer(Y1~CO2*Light+(1|Treatment:Tank), data=data, REML=FALSE, na.action=na.exclude) R1.table<-anova(final.mod, type=2) #ANOVA table R1.mod1<-summary(mod1) R1.table ``` **total biomass cm^-2^** ```{r} ############################################ ** no colony effect ## BIOMASS ################################# ** tank effect ############################################ ** no treatment effect mod1<-lmer(Y2~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y2~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y2~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand2<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is full model (=mod1) # tanks signficant, colony effect 0.1 (retain), no effect fixed effect final.mod<-lmer(Y2~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, REML=FALSE, na.action=na.exclude) R2.table<-anova(final.mod, type=2) #ANOVA table R2.mod1<-summary(mod1) R2.table ``` **Symbiodinium cells cm^-2^** ```{r} ############################################ ** no colony effect ## SYMBIODINIUM CELLS ###################### ** tank effects ############################################ ** no treatment effect ** mod1<-lmer(Y3~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y3~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y3~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand3<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is colony dropped, (=mod2) # tanks signficant, NO colony effect, NO effect fixed effect final.mod<-lmer(Y3~CO2*Light+(1|Treatment:Tank), REML=FALSE, data=data, na.action=na.exclude) R3.table<-anova(final.mod, type=2) #ANOVA table R3.mod1<-summary(mod1) R3.table ``` **chlorophyll *a* cm^-2^** ```{r} ############################################ ** effect of colony ## CHLOROPHYLL A ########################## ** effect of tank ############################################ ** TREATMENT EFFECT = LIGHT ** mod1<-lmer(Y4~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y4~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y4~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand4<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is full model, (=mod1) # tanks signficant, colony significant, fixed effect significant final.mod<-lmer(Y4~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, REML=FALSE, na.action=na.exclude) R4.table<-anova(final.mod, type=2) #ANOVA table R4.mod1<-summary(mod1) R4.table ``` **chlorophyll *c~2~* cm^-2^** ```{r} ############################################ ** effect of colony ## CHLOROPHYLL C2 ######################### ** effect of tank ############################################ ** TREATMENT EFFECT = LIGHT ** mod1<-lmer(Y5~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, REML=FALSE, na.action=na.exclude) mod2<-lmer(Y5~CO2*Light+(1|Treatment:Tank), data=data, REML=FALSE, na.action=na.exclude) mod3<-lmer(Y5~CO2*Light+(1|Colony), data=data, REML=FALSE, na.action=na.exclude) rand5<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is full model, (=mod1) # tanks signficant, colony significant, fixed effect significant final.mod<-lmer(Y5~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, REML=FALSE, na.action=na.exclude) R5.table<-anova(final.mod, type=2) #ANOVA table R5.mod1<-summary(mod1) R5.table ``` **chlorophyll *a* cell^-1^** ```{r} ############################################ ** no colony effect ## CHLOROPHYLL A /cell ##################### ** no tank effects ############################################ ** no treatment effect mod1<-lmer(Y6~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y6~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y6~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand6<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is full model, (=mod1) # tanks signficant, slight colony effect, NO fixed effects final.mod<-lmer(Y6~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, REML=FALSE, na.action=na.exclude) R6.table<-anova(final.mod, type=2) #ANOVA table R6.mod1<-summary(mod1) R6.table ``` **chlorophyll *c~2~* cell^-1^** ```{r} ############################################ ** colony effect ## CHLOROPHYLL C2/cell ##################### ** tank effect ############################################ ** no treatment effect mod1<-lmer(Y7~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y7~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y7~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand7<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is full model, (=mod1) # colony significant, tanks signficant, NO fixed effects final.mod<-lmer(Y7~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, REML=FALSE, na.action=na.exclude) R7.table<-anova(final.mod, type=2) #ANOVA table R7.mod1<-summary(mod1) R7.table ``` **protein cm^-2^** ```{r} ############################################ ** slight colony effect (p=0.10) ## PROTEIN ################################# ** slight tank effect (p=0.07) ############################################ ** TRTMNT EFFECTS: LIGHT, CO2 X LIGHT mod1<-lmer(Y8~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y8~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y8~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand8<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1,mod2, mod3) # best 'full model' (=mod1) ## final model is full model final.mod<-lmer(Y8~CO2*Light+(1|Treatment:Tank)+(1|Colony), REML=FALSE, data=data, na.action=na.exclude) R8.table<-anova(final.mod, type=2) R8.mod1<-summary(mod1) ### posthoc posthoc<-lsmeans(final.mod, pairwise~Light*CO2, adjust="tukey") # compact letter display for plotting contrasts<-cld(posthoc, Letters=letters); contrasts R8.table ``` **carbohydrate cm^-2^** ```{r} ############################################ ** no colony effect ## CARBOHYDRATE ############################ ** tank effect ############################################ ** no treatment effect mod1<-lmer(Y9~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y9~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y9~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand9<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is colony dropped model, (=mod2) # tanks signficant, colony significant, NO fixed effects final.mod<-lmer(Y9~CO2*Light+(1|Treatment:Tank), REML=FALSE, data=data, na.action=na.exclude) R9.table<-anova(final.mod, type=2) #ANOVA table R9.mod1<-summary(mod1) R9.table ``` **lipids cm^-2^** ```{r} ############################################ ** no colony effect ## LIPIDS ################################## ** slight tank effect (p=0.10) ############################################ ** no effect of treatment mod1<-lmer(Y10~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y10~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y10~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand10<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is colony dropped model, (=mod2) # tanks signficant, NO colony effect, NO fixed effects final.mod<-lmer(Y10~CO2*Light+(1|Treatment:Tank), data=data, REML=FALSE, na.action=na.exclude) R10.table<-anova(final.mod, type=2) #ANOVA table R10.mod1<-summary(mod1) R10.table ``` **energy content cm^-2^** ```{r} ############################################ ** no effect of colony ## ENERGY CONTENT ########################## ** slight effect of tank (0.1) ############################################ ** no treatment effect mod1<-lmer(Y11~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y11~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y11~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand11<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is colony dropped model, (=mod2) ## final model: tanks signficant, NO colony effect, NO fixed effects final.mod<-lmer(Y11~CO2*Light+(1|Treatment:Tank), REML=FALSE, data=data, na.action=na.exclude) R11.table<-anova(final.mod, type=2) R11.mod1<-summary(mod1) R11.table ``` ### Biomass-normalized models **calcification gdw^-1^** ```{r} ## **************************************## ## BIOMASS NORMALIZED ##################### ## **************************************## ############################################ ** no colony effect ## CALCIFICATION ########################### ** tank effect ############################################ ** no treatment effect mod1<-lmer(Y12~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y12~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y12~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand12<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is colony dropped model, (=mod2) # tank effect, NO colony effect, NO fixed effects final.mod<-lmer(Y12~CO2*Light+(1|Treatment:Tank), data=data, REML=FALSE, na.action=na.exclude) R12.table<-anova(final.mod, type=2) #ANOVA table R12.mod1<-summary(mod1) R12.table ``` **protein gdw^-1^** ```{r} ############################################ ** colony effect ## PROTEIN ################################ ** tank effect ############################################ ** no treatment effect mod1<-lmer(Y13~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y13~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y13~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand13<-rand(mod1) # p values for random effects mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is full model, (=mod1) ## final model final.mod<-lmer(Y13~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, REML=FALSE, na.action=na.exclude) R13.table<-anova(final.mod, type=2) #ANOVA table R13.mod1<-summary(mod1) R13.table ``` **carbohydrates gdw^-1^** ```{r} ############################################ ** colony effect ## CARBOHYDRATES ########################## ** tank effect ############################################ ** TREATMENT EFFECT = LIGHT mod1<-lmer(Y14~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y14~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y14~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand14<-rand(mod1) mod.AIC<-AIC(mod1, mod2, mod3) # best is full colony model, (=mod1) # tank effect, NO colony effect, NO fixed effects final.mod<-lmer(Y14~CO2*Light+(1|Treatment:Tank)+(1|Colony), REML=FALSE, data=data, na.action=na.exclude) R14.table<-anova(final.mod, type=2) #ANOVA table R14.mod1<-summary(mod1) R14.table ``` **lipids gdw^-1^** ```{r} ############################################ ** no colony effect ## LIPIDS ################################# ** no tank effect ############################################ ** TREATMENT EFFECT = CO2 mod1<-lmer(Y15~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y15~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y15~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand15<-rand(mod1) mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is drop colony model, (=mod2) #final model final.mod<-lmer(Y15~CO2*Light+(1|Treatment:Tank), data=data, REML=FALSE, na.action=na.exclude) R15.table<-anova(final.mod, type=2) R15.mod1<-summary(mod1) R15.table ``` **energy content gdw^-1^** ```{r} ############################################ ** no colony effect ## ENERGY CONTENT ######################### ** no tank effect ############################################ ** TREATMENT EFFECT = CO2 mod1<-lmer(Y16~CO2*Light+(1|Treatment:Tank)+(1|Colony), data=data, na.action=na.exclude) mod2<-lmer(Y16~CO2*Light+(1|Treatment:Tank), data=data, na.action=na.exclude) mod3<-lmer(Y16~CO2*Light+(1|Colony), data=data, na.action=na.exclude) rand16<-rand(mod1) mod.AIC<-AIC(mod1, mod2, mod3) # best 'mixed model' is tank dropped model, (=mod2) # NO tank effect, NO colony effect, NO fixed effects final.mod<-lmer(Y16~CO2*Light+(1|Treatment:Tank), data=data, REML=FALSE, na.action=na.exclude) R16.table<-anova(final.mod, type =2) #ANOVA table R16.mod1<-summary(mod1) R16.table ``` #### Figures To generate plots, summary tables of mean, SE, and n were generated. Plots below show response variable mean +/- SE for corals exposed to CO~2~-Light treatments. Responses are normalized to area (cm^2^) or biomass (gdw). ```{r, include=FALSE, results="hide"} ################################## ## means for area normalized df ## ################################## calcif.m.cm2<-aggregate(calcif.cm2~CO2+Light, df, mean); calcif.m.cm2 biomass.m.cm2<-aggregate(biomass~CO2+Light, df, mean); biomass.m.cm2 zoox.m.cm2<-aggregate(zoox~CO2+Light, df, mean); zoox.m.cm2 chla.m.cm2<-aggregate(chla~CO2+Light, df, mean); chla.m.cm2 chlc2.m.cm2<-aggregate(chlc2~CO2+Light, df, mean); chlc2.m.cm2 chlacell.m.cm2<-aggregate(chlacell~CO2+Light, df, mean); chlacell.m.cm2 chlc2cell.m.cm2<-aggregate(chlc2cell~CO2+Light, df, mean); chlc2cell.m.cm2 prot.m.cm2<-aggregate(prot.cm2~CO2+Light, df, mean); prot.m.cm2 carb.m.cm2<-aggregate(carb.cm2~CO2+Light, df, mean); carb.m.cm2 lipid.m.cm2<-aggregate(lipid.cm2~CO2+Light, df, mean); lipid.m.cm2 energy.m.cm2<-aggregate(energy.cm2~CO2+Light, df, mean); energy.m.cm2 df.m.cm2<-cbind(calcif.m.cm2, biomass.m.cm2[c(3,0)], zoox.m.cm2[c(3,0)], chla.m.cm2[c(3,0)], chlc2.m.cm2[c(3,0)], chlacell.m.cm2[c(3,0)], chlc2cell.m.cm2[c(3,0)], prot.m.cm2[c(3,0)], carb.m.cm2[c(3,0)], lipid.m.cm2[c(3,0)], energy.m.cm2[c(3,0)]) df.m.cm2 ##################################### # SE ##################################### calcif.SE.cm2<-aggregate(calcif.cm2~CO2+Light, df, std.error); calcif.SE.cm2 biomass.SE.cm2<-aggregate(biomass~CO2+Light, df, std.error); biomass.SE.cm2 zoox.SE.cm2<-aggregate(zoox~CO2+Light, df, std.error); zoox.SE.cm2 chla.SE.cm2<-aggregate(chla~CO2+Light, df, std.error); chla.SE.cm2 chlc2.SE.cm2<-aggregate(chlc2~CO2+Light, df, std.error); chlc2.SE.cm2 chlacell.SE.cm2<-aggregate(chlacell~CO2+Light, df, std.error); chlacell.SE.cm2 chlc2cell.SE.cm2<-aggregate(chlc2cell~CO2+Light, df, std.error); chlc2cell.SE.cm2 prot.SE.cm2<-aggregate(prot.cm2~CO2+Light, df, std.error); prot.SE.cm2 carb.SE.cm2<-aggregate(carb.cm2~CO2+Light, df, std.error); carb.SE.cm2 lipid.SE.cm2<-aggregate(lipid.cm2~CO2+Light, df, std.error); lipid.SE.cm2 energy.SE.cm2<-aggregate(energy.cm2~CO2+Light, df, std.error); energy.SE.cm2 df.SE.cm2<-cbind(calcif.SE.cm2, biomass.SE.cm2[c(3,0)], zoox.SE.cm2[c(3,0)], chla.SE.cm2[c(3,0)], chlc2.SE.cm2[c(3,0)], chlacell.SE.cm2[c(3,0)], chlc2cell.SE.cm2[c(3,0)], prot.SE.cm2[c(3,0)], carb.SE.cm2[c(3,0)], lipid.SE.cm2[c(3,0)], energy.SE.cm2[c(3,0)]) df.SE.cm2 ##################################### # sample size ##################################### calcif.n.cm2<-aggregate(calcif.cm2~CO2+Light, df, length); calcif.n.cm2 biomass.n.cm2<-aggregate(biomass~CO2+Light, df, length); biomass.n.cm2 zoox.n.cm2<-aggregate(zoox~CO2+Light, df, length); zoox.n.cm2 chla.n.cm2<-aggregate(chla~CO2+Light, df, length); chla.n.cm2 chlc2.n.cm2<-aggregate(chlc2~CO2+Light, df, length); chlc2.n.cm2 chlacell.n.cm2<-aggregate(chlacell~CO2+Light, df, length); chlacell.n.cm2 chlc2cell.n.cm2<-aggregate(chlc2cell~CO2+Light, df, length); chlc2cell.n.cm2 prot.n.cm2<-aggregate(prot.cm2~CO2+Light, df, length); prot.n.cm2 carb.n.cm2<-aggregate(carb.cm2~CO2+Light, df, length); carb.n.cm2 lipid.n.cm2<-aggregate(lipid.cm2~CO2+Light, df, length); lipid.n.cm2 energy.n.cm2<-aggregate(energy.cm2~CO2+Light, df, length); energy.n.cm2 df.n.cm2<-cbind(calcif.n.cm2, biomass.n.cm2[c(3,0)], zoox.n.cm2[c(3,0)], chla.n.cm2[c(3,0)], chlc2.n.cm2[c(3,0)], chlacell.n.cm2[c(3,0)], chlc2cell.n.cm2[c(3,0)], prot.n.cm2[c(3,0)], carb.n.cm2[c(3,0)], lipid.n.cm2[c(3,0)], energy.n.cm2[c(3,0)]) df.n.cm2 ##### create dfframe with means and SE df.cm2<-cbind(df.m.cm2, df.SE.cm2[c(3:13,0)]); colnames(df.cm2) <-c("CO2", "Light", "calcif", "biomass", "zoox", "chla", "chlc2", "chlacell", "chlc2cell", "prot", "carb", "lipid", "energy", "calcifSE", "biomassSE", "zooxSE", "chlaSE", "chlc2SE","chlacellSE", "chlc2cellSE", "protSE", "carbSE", "lipidSE", "energySE") ###################################### ## means for biomass normalized df ## ###################################### calcif.m.gdw<-aggregate(calcif.gdw~CO2+Light, df, mean); calcif.m.gdw prot.m.gdw<-aggregate(prot.gdw~CO2+Light, df, mean); prot.m.gdw carb.m.gdw<-aggregate(carb.gdw~CO2+Light, df, mean); carb.m.gdw lipid.m.gdw<-aggregate(lipid.gdw~CO2+Light, df, mean); lipid.m.gdw energy.m.gdw<-aggregate(energy.gdw~CO2+Light, df, mean); energy.m.gdw df.m.gdw<-cbind(calcif.m.gdw, prot.m.gdw[c(3,0)], carb.m.gdw[c(3,0)], lipid.m.gdw[c(3,0)], energy.m.gdw[c(3,0)]) df.m.gdw ##################################### # SE ##################################### calcif.SE.gdw<-aggregate(calcif.gdw~CO2+Light, df, std.error); calcif.SE.gdw prot.SE.gdw<-aggregate(prot.gdw~CO2+Light, df, std.error); prot.SE.gdw carb.SE.gdw<-aggregate(carb.gdw~CO2+Light, df, std.error); carb.SE.gdw lipid.SE.gdw<-aggregate(lipid.gdw~CO2+Light, df, std.error); lipid.SE.gdw energy.SE.gdw<-aggregate(energy.gdw~CO2+Light, df, std.error); energy.SE.gdw df.SE.gdw<-cbind(calcif.SE.gdw, prot.SE.gdw[c(3,0)], carb.SE.gdw[c(3,0)], lipid.SE.gdw[c(3,0)], energy.SE.gdw[c(3,0)]) df.SE.gdw ##################################### # sample size ##################################### calcif.n.gdw<-aggregate(calcif.gdw~CO2+Light, df, length); calcif.n.gdw prot.n.gdw<-aggregate(prot.gdw~CO2+Light, df, length); prot.n.gdw carb.n.gdw<-aggregate(carb.gdw~CO2+Light, df, length); carb.n.gdw lipid.n.gdw<-aggregate(lipid.gdw~CO2+Light, df, length); lipid.n.gdw energy.n.gdw<-aggregate(energy.gdw~CO2+Light, df, length); energy.n.gdw df.n.gdw<-cbind(calcif.n.gdw, prot.n.gdw[c(3,0)], carb.n.gdw[c(3,0)], lipid.n.gdw[c(3,0)], energy.n.gdw[c(3,0)]) df.n.gdw ##################################### df.gdw<-cbind(df.m.gdw, df.SE.gdw[c(3:7,0)]); colnames(df.gdw) <-c("CO2", "Light", "calcif", "prot", "carb", "lipid","energy", "calcifSE", "protSE", "carbSE", "lipidSE", "energySE") df.gdw ``` Figures here show calcification (surface area (cm~2~) and biomass normalized (gdw)), symbiont densities (cells/cm^2^), tissue biomass (mg/cm^2^), chlorophyll *a* and *c~2~* (ug/cm^2^), and chlorophyll *a* and *c~2~* per symbiont cell (pg/cell). These figures show data found in **"Figure 2"**. ```{r, results="hide"} ##################################### pd <- position_dodge(0.1) #offset for error bars ##################################### setwd('markdown/figures') # calcification per area (mg CaCO3/cm2/d) Fig.calcif<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$calcif, group=df.cm2$CO2, fill=df.cm2$CO2)) + geom_errorbar(aes(ymin=df.cm2$calcif-df.cm2$calcifSE, ymax=df.cm2$calcif+df.cm2$calcifSE), size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Calcification"~mg~CaCO[3]~cm^-2~d^-1), sep="")) + scale_x_discrete(limits=c("LL", "HL")) + scale_y_continuous(limits=c(0.22, 0.4), breaks=c(0.22, 0.26, 0.30, 0.34, 0.38))+ theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15)) + theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(aspect.ratio=1)+ theme(panel.background = element_rect(colour = "black", size=1)) + theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.calcif.pdf", height=8, width=8, units="cm") #calcification standardized to biomass (mg CaCO3/gdw/d) Fig.calcif.gdw<-ggplot(data=df.gdw, aes(x=df.gdw$Light, y=df.gdw$calcif, group=df.gdw$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.gdw$calcif-df.gdw$calcifSE, ymax=df.gdw$calcif+df.gdw$calcifSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ coord_cartesian(ylim=c(10, 22.5)) + xlab("") + ylab(expression(paste("Calcification" ~(mg~CaCO[3]~gdw~d^-1), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(10, 22), breaks=c(10, 12.5, 15, 17.5, 20, 22.5)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.calcif.gdw.pdf", height=8, width=8, units="cm") ##################################### # AFDW (biomass) normalized to area (mgdw/cm2) Fig.AFDW<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$biomass, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$biomass-df.cm2$biomassSE, ymax=df.cm2$biomass+df.cm2$biomassSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Biomass" ~(mg~cm^-2), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(15, 26), breaks=c(16, 18, 20, 22, 24, 26)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.AFDW.pdf", height=8, width=8, units="cm") ##################################### # symbiont densities (symbionts/cm2) Fig.symb<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$zoox/10^6, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$zoox/10^6-df.cm2$zooxSE/10^6, ymax=df.cm2$zoox/10^6+df.cm2$zooxSE/10^6),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste(~italic("Symbiodinium")~(10^6~cells~cm^-2), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0.5, 1.5), breaks=c(0.5, 0.75, 1.0, 1.25, 1.5)) + theme_classic() + theme(legend.key = element_blank())+ theme(legend.justification=c(1,1), legend.position=c(1,1))+ theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.symb.pdf", height=8, width=8, units="cm") ##################################### ##### chlorophyll a (ug/cm2) Fig.chla<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$chla, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$chla-df.cm2$chlaSE, ymax=df.cm2$chla+df.cm2$chlaSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Chlorophyll", ~italic("a"), ~(mu*g~cm^-2), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0, 6), breaks=c(0, 1, 2, 3, 4, 5, 6)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.chla.pdf", height=8, width=8, units="cm") ##################################### # chlorophyll c2 concentration (ug chl c2/cm2) Fig.chlc<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$chlc2, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$chlc2-df.cm2$chlc2SE, ymax=df.cm2$chlc2+df.cm2$chlc2SE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=22) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Chlorophyll", ~italic("c"[2]), ~(mu*g~cm^-2), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0, 6), breaks=c(0, 1, 2, 3, 4, 5, 6)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.chlc.pdf", height=8, width=8, units="cm") ##################################### # chlorophyll a per symbiont cell (pg chla/cell) Fig.chlacell<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$chlacell, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$chlacell-df.cm2$chlacellSE, ymax=df.cm2$chlacell+df.cm2$chlacellSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Chlorophyll", ~italic("a"), ~(pg~cell^-1), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0, 6), breaks=c(0, 1, 2, 3, 4, 5, 6)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.chlacell.pdf", height=8, width=8, units="cm") ##################################### # chlorophyll c2 per symbiont cell (pg chl c2/cell) Fig.chlccell<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$chlc2cell, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$chlc2cell-df.cm2$chlc2cellSE, ymax=df.cm2$chlc2cell+df.cm2$chlc2cellSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=22) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Chlorophyll", ~italic("c"[2]), ~(pg~cell^-1), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0, 6), breaks=c(0, 1, 2, 3, 4, 5, 6)) + theme_classic() + theme(legend.key = element_blank())+ theme(legend.justification=c(1,1), legend.position=c(1,1))+ theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.chlccell.pdf", height=8, width=8, units="cm") ``` ```{r,fig.width=3, fig.height=3, echo=FALSE} Fig.calcif Fig.calcif.gdw Fig.AFDW Fig.symb Fig.chla Fig.chlc Fig.chlacell Fig.chlccell ``` Figures here show biomass-normalized energy reserves (g/gdw): protein (g/gdw), carbohydrates (g/gdw), lipids (g/gdw1), and tissue energy content (kJ/gdw) These figures show data found in **"Figure 3"**. ```{r, results="hide"} # protein concentration (g protein/gdw) setwd('markdown/figures') Fig.prot.gdw<-ggplot(data=df.gdw, aes(x=df.gdw$Light, y=df.gdw$prot, group=df.gdw$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.gdw$prot-df.gdw$protSE, ymax=df.gdw$prot+df.gdw$protSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste(Protein ~(g~gdw^-1), sep=""))) + scale_y_continuous(limits=c(0.028, 0.045), breaks=c(0.028, 0.032, 0.036, 0.040, 0.044)) + scale_x_discrete(limits=c("LL", "HL"))+ theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.prot.gdw.pdf", height=8, width=8, units="cm") # carbohydrates (g carbs/gdw) Fig.carb.gdw<-ggplot(data=df.gdw, aes(x=df.gdw$Light, y=df.gdw$carb, group=df.gdw$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.gdw$carb-df.gdw$carbSE, ymax=df.gdw$carb+df.gdw$carbSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste(Carbohydrates ~(g~gdw^-1), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0.007, 0.015), breaks=c(0.007, 0.009, 0.011, 0.013, 0.015)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.carb.gdw.pdf", height=8, width=8, units="cm") # lipid concentration (g lipids/gdw) Fig.lip.gdw<-ggplot(data=df.gdw, aes(x=df.gdw$Light, y=df.gdw$lipid, group=df.gdw$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.gdw$lipid-df.gdw$lipidSE, ymax=df.gdw$lipid+df.gdw$lipidSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste(Lipids ~(g~gdw^-1), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0.2, 0.45), breaks=c(0.2, 0.25, 0.3, 0.35, 0.4, 0.45)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.lip.gdw.pdf", height=8, width=8, units="cm") # energy content (J /gdw) Fig.kJ.gdw<-ggplot(data=df.gdw, aes(x=df.gdw$Light, y=df.gdw$energy, group=df.gdw$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.gdw$energy-df.gdw$energySE, ymax=df.gdw$energy+df.gdw$energySE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Energy"~(kJ~gdw^-1), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(10, 20), breaks=c(10, 12, 14, 16, 18, 20)) + theme_classic() + theme(legend.key = element_blank())+ theme(legend.justification=c(0,1), legend.position=c(1.2,1))+ theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.kJ.gdw.pdf", height=8, width=8, units="cm") ``` ```{r,fig.width=3, fig.height=3, echo=FALSE} Fig.prot.gdw Fig.carb.gdw Fig.lip.gdw Fig.kJ.gdw ``` Figures here show area-normalized energy reserves (mg/cm^2^): protein (mg/cm^2^), carbohydrates (mg/cm^2^), lipids (mg/cm^2^), and tissue energy content (kJ/cm2) These figures show data found in **"Supplemental Figure 1"**. ```{r} ##################################### setwd('markdown/figures') # protein concentrations (mg protein/cm2) Fig.prot<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$prot, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$prot-df.cm2$protSE, ymax=df.cm2$prot+df.cm2$protSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Protein" ~(mg~cm^-2), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0.5, 0.9), breaks=c(0.5, 0.6, 0.7, 0.8, 0.9)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.prot.pdf", height=8, width=8, units="cm") ##################################### # carbohydrate concentration (mg carbs/cm2) Fig.carb<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$carb, group=CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$carb-df.cm2$carbSE, ymax=df.cm2$carb+df.cm2$carbSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ coord_cartesian(ylim=c(0.15, 0.30)) + xlab("") + ylab(expression(paste("Carbohydrates" ~(mg~cm^-2), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0.15, 0.30), breaks=c(0.15, 0.18, 0.21, 0.24, 0.27, 0.3)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.carb.pdf", height=8, width=8, units="cm") ##################################### # lipid concentrations (mg lipids/cm2) Fig.lip<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$lipid, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$lipid-df.cm2$lipidSE, ymax=df.cm2$lipid+df.cm2$lipidSE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])), guide=FALSE)+ xlab("") + ylab(expression(paste("Lipids" ~(mg~cm^-2), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(1, 5), breaks=c(1, 2, 3, 4, 5)) + theme_classic() + theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.lip.pdf", height=8, width=8, units="cm") ##################################### #energy kJ/cm2 Fig.kJ<-ggplot(data=df.cm2, aes(x=df.cm2$Light, y=df.cm2$energy, group=df.cm2$CO2, fill=CO2)) + geom_errorbar(aes(ymin=df.cm2$energy-df.cm2$energySE, ymax=df.cm2$energy+df.cm2$energySE),size=.5, width=0, position=pd) + geom_line(position=pd, size=.5) + geom_point(aes(fill=CO2), position=pd, size=5, pch=21) + scale_fill_manual(values=c('black','white'), labels=rep(expression(ACO[2], HCO[2])))+ coord_cartesian(ylim=c(0.18, 0.38)) + xlab("") + ylab(expression(paste("Energy" ~(kJ~cm^-2), sep=""))) + scale_x_discrete(limits=c("LL", "HL"))+ scale_y_continuous(limits=c(0.18, 0.38), breaks=c(0.2, 0.24, 0.28, 0.32, 0.36)) + theme_classic() + theme(panel.border = element_rect(fill=NA, color = "black", size = .1)) + theme(legend.key = element_blank()) + theme(legend.justification=c(0,1), legend.position=c(1,1))+ theme(text=element_text(size=12))+ theme(axis.title.x = element_text(face="bold", size=15))+ theme(legend.text=element_text(size=15)) + theme(legend.title = element_blank())+ theme(panel.background = element_rect(colour = "black", size=1))+ theme(aspect.ratio=1)+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.y=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm")), axis.text.x=element_text(margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"))) ggsave(file="Fig.kJ.pdf", height=8, width=8, units="cm") ``` ```{r,fig.width=3, fig.height=3, echo=FALSE} Fig.prot Fig.carb Fig.lip Fig.kJ ```