Experimental design

  • Treatments were created in 24 tanks of 6 replicate treatents:
  • Low Light-Ambient pCO2 (LL-ACO2), Low light-High pCO2 (LL-HCO2), High light-Ambient pCO2 (HL-ACO2), and High light-High pCO2 (HL-HCO2).
  • 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 (pCO2 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 pCO2 among all treatment tanks and among replicate treatment tanks within a CO2 treatment (ambient or high).

Attach packages

## 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)

Attach the first dataset for seawater chemistry.

##################################################
## 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")

pH and CO2 treatment figure

Generate a figure of mean +/- SE pHT and pCO2 for each treatment tank.

# 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 (AT or ALK), CO2, and pHT.
  • First, testing for treatment effects, ie., are there two CO2 treatments and is AT stable across all tanks?
####################################################################
#  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

## Analysis of Variance Table
## 
## Response: ALK
##            Df Sum Sq Mean Sq F value Pr(>F)
## CO2_Treat   1   1136 1135.84  2.5833 0.1101
## Residuals 151  66392  439.68

pCO2 across all tanks

## Analysis of Variance Table
## 
## Response: pCO2
##            Df   Sum Sq  Mean Sq F value    Pr(>F)    
## CO2_Treat   1 10299281 10299281  338.83 < 2.2e-16 ***
## Residuals 151  4589914    30397                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pHT across all tanks

## Analysis of Variance Table
## 
## Response: pHT
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## CO2_Treat   1 3.0899 3.08992  402.26 < 2.2e-16 ***
## Residuals 151 1.1599 0.00768                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • There are differences in treatments (ACO2 and HCO2) and total alkalinity does not differ among treatments.
  • NEXT, testing whether replicate tanks are different, based on CO2 treatments (Ambient CO2 vs. High CO2). To do this, divide dataframe into High and Ambient pCO2 treatments and examine replicate tanks.

Testing replicate tanks

FIRST, testing for differences among replicate Ambient pCO2 tanks…

######### 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 ACO2 tanks

## Analysis of Variance Table
## 
## Response: ALK
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Tank      11   365.4   33.22  0.0788      1
## Residuals 72 30338.4  421.37

pCO2 across ACO2 tanks

## Analysis of Variance Table
## 
## Response: pCO2
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Tank      11  98427  8947.9  1.8557 0.06011 .
## Residuals 72 347164  4821.7                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pHT across ACO2 tanks

## Analysis of Variance Table
## 
## Response: pHT
##           Df   Sum Sq   Mean Sq F value  Pr(>F)  
## Tank      11 0.068451 0.0062228  1.7293 0.08396 .
## Residuals 72 0.259094 0.0035985                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

…THEN, testing for differences among replicate High pCO2 tanks…

########## 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 HCO2 tanks

## Analysis of Variance Table
## 
## Response: ALK
##           Df Sum Sq Mean Sq F value Pr(>F)
## Tank       9   1800  199.99  0.3482 0.9544
## Residuals 59  33889  574.39

pCO2 across HCO2 tanks

## Analysis of Variance Table
## 
## Response: pCO2
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Tank       9  177530   19726  0.2934 0.9739
## Residuals 59 3966793   67234

pHT across HCO2 tanks

## Analysis of Variance Table
## 
## Response: pHT
##           Df  Sum Sq   Mean Sq F value Pr(>F)
## Tank       9 0.04425 0.0049163  0.3681 0.9458
## Residuals 59 0.78810 0.0133576

Physiological Responses

  • Response variables are normalized to area (cm2) 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 c2 is normalized to area.
  • Photopigments (chlorophyll a and c2) are also normalized to symbiont cell densities.

First, normalize data found in raw physiology data file (code details below)

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
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 <<<normal
Y16<-energy.gdw #energy/gdw <<<<normal

###############################
## transformation assistance ##
###############################

# boxcox(**specify Y**~Light+CO2, lambda=seq(-1, 2, length=10))

###############################
###### normality test #########
###############################
df.tests<-cbind(df[1:4], Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8, Y9, Y10, Y11, Y12, Y13, Y14, Y15, Y16)
colnames(df.tests)<-c("CO2", "Light", "Treatment", "Tank", "calcif.cm2.trans", "biomass.trans", "zoox.trans", "chla", "chlc2.trans", "chlacell.trans", "chlc2cell.trans", "prot.cm2.trans", "carbs.cm2.trans", "lipids.cm2.trans", "energy.cm2.trans", "calcif.gdw.trans",  "prot.gdw", "carbs.gdw.trans", "lipids.gdw", "energy.gdw")

names(df.tests)
##  [1] "CO2"              "Light"            "Treatment"       
##  [4] "Tank"             "calcif.cm2.trans" "biomass.trans"   
##  [7] "zoox.trans"       "chla"             "chlc2.trans"     
## [10] "chlacell.trans"   "chlc2cell.trans"  "prot.cm2.trans"  
## [13] "carbs.cm2.trans"  "lipids.cm2.trans" "energy.cm2.trans"
## [16] "calcif.gdw.trans" "prot.gdw"         "carbs.gdw.trans" 
## [19] "lipids.gdw"       "energy.gdw"

Code chunk below shows a for-loop that will perform tests for assumptions of ANOVA (normal distribution of residuals, equal variance) and generates diagnostic plots

names(df.tests)
for(i in 5:20){
  Y<-df.tests[,i]
  full<-lmer(Y~CO2*Light+(1|Treatment/Tank)+(1|Colony), REML=FALSE, na.action=na.exclude)
  R <- resid(full) #save glm residuals
  Y.shapiro <- shapiro.test(R) #runs a normality test on residuals
  print(Y.shapiro) # null = normally distrubuted (P<0.05 = non-normal
  
  Levene <- leveneTest(Y~CO2*Light); print(Levene)
  
  op<-par(mfrow = c(2,2), mar=c(5,4,1,2), pty="sq")
  plot(full, add.smooth = FALSE, which=1)
  QQ <- qqnorm(R, main = colnames(df.tests)[i]) 
  QQline <- qqline(R)
  hist(R, xlab="Residuals", main = colnames(df.tests)[i])
  plot(Treatment, R, xlab="CO2-Light Treatment", ylab="Residuals")
  plot(Tank, R, xlab="Tank", ylab="Residuals")
}

Statistical Models

In each case, dependent variables were tested in a complete mixed effect model testing CO2*Light effects with Tank as a nested factor nested within CO2-Light Treatment (1|Treatment/Tank) and coral colony (i.e., genet) as a random factor (1|Colony). Models were all fit using restricted maximum likelihood when comparing models with different random effects. The final model is reported with maximum likelihood with anova tables generated from Type II Sum of Squares

  • mod<-lmer(Y1CO2~*Light+(1|Treatment:Tank)+(1|Colony), REML=FALSE, na.action=na.exclude)

The full model was tested against a model with Colony dropped and a separate model with Tank dropped. The best fit model was determined from log-likelihood tests of the three models. All final models contained at least 1 random effect.

Principal Component Analysis

The data matrix consists of corals where all biomass properties were measured (proteins, lipids, carbs, total biomass, and biomass energy content). Two dataset dataframes were generated: one for area-normalized responses, and the other for biomass normalized responses. In each dataframe, total biomass (AFDW cm-2) was included, as this is the only way to measure the totality of biomass.

Area normalized PCA

These data are all normalized to surface area. First a dataframe is generated with specific varaibles of interest. Two outliers are removed. The table generated shows importance of components with standard deviation, variance proportion captured, and the proportion of cumulative variance.

####################################
## PCA on OA x Light data: area normalized data ##
## using calcif, biomass, kJ, lipids, carbs, protein ##
####################################
library(devtools)
install_github("ggbiplot", "vqv")
library(ggbiplot)
require(graphics)

# names(df) # normalized dataframe

##############################
####### area PCA #############
df.2<-df[, c(1:6,9,11,17,19,21,23)] ## area normalized data
df.area<-na.omit(df.2) # NAs removed
df.area.outlier.removed<-df.area[c(-16, -69), ] #outliers

pca.area.df<-(df.area.outlier.removed[,c(3,7:12)]) # all data with only "Treatment" as category

#only 83 rows remain
# nrow(pca.area.df)

colnames(pca.area.df)<-c("Treatment", "calcification", "biomass", "protein","carbohydrates", "lipids", "energy")


############
############
# PCA area 
# apply PCA - scale. = TRUE for center as advised
PC.area <- prcomp(pca.area.df[,-1], center = TRUE, scale= TRUE)  
# correlatin matrix helps to standardize the data and account for different scales

# print(PC.area) <<< to show loadings
PC.area.summary<-(summary(PC.area)) #signficance of PCs: proportion of variance captured, cumulative variance

#In general, we are interested in keeping only those principal components whose eigenvalues are greater than 1. 

# the SDEV^2 is the eignevalue
ev.area<-PC.area$sdev^2 # PC1-2 have eignevalues >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.
######################
# 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.
######################
# 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 CO2, 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.

# 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

##################
# 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)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##            Sum Sq Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       1.00371 1.00371     1 21.369 0.48768 0.4925
## Light     0.51029 0.51029     1 21.204 0.24794 0.6237
## CO2:Light 0.06967 0.06967     1 21.223 0.03385 0.8558

Results for PC2-area lmer: no effect of treatment on PC2-area

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)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##            Sum Sq Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       2.89892 2.89892     1 21.912 2.70394 0.1144
## Light     0.37208 0.37208     1 21.729 0.34705 0.5619
## CO2:Light 0.93994 0.93994     1 21.741 0.87672 0.3594

PCA biomass normalization

This code will run the identical analysis used on PCA-area, but now with repsonses normalized to biomass (gdw).

# 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
######################
# 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
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

# 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

##################
#### final models
# PC1
mod1<-lmer(PC1~CO2*Light+(1|Treatment:Tank), data=PCdata, REML=FALSE)
PC1.gdw<-anova(mod1, type=2) 
PC1.gdw
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##            Sum Sq Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       0.27910 0.27910     1 22.514 0.16421 0.6891
## Light     2.18454 2.18454     1 22.381 1.28526 0.2689
## CO2:Light 0.13932 0.13932     1 22.411 0.08197 0.7773

PC2-biomass model results: effect of pCO2 on PC2-biomass

#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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##           Sum Sq Mean Sq NumDF  DenDF F.value  Pr(>F)  
## CO2       9.7262  9.7262     1 22.419  5.5019 0.02823 *
## Light     0.2359  0.2359     1 22.266  0.1334 0.71835  
## CO2:Light 0.0189  0.0189     1 22.287  0.0107 0.91863  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

##################
#### 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

##################
#### 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).

Area-normalized models

calcification cm-2

## ***********************************##
## 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##              Sum Sq   Mean Sq NumDF  DenDF  F.value Pr(>F)
## CO2       0.0041979 0.0041979     1 22.057 0.274883 0.6053
## Light     0.0010740 0.0010740     1 21.943 0.070325 0.7933
## CO2:Light 0.0004808 0.0004808     1 21.952 0.031482 0.8608

total biomass cm-2

############################################ ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##               Sum Sq    Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       0.00000302 0.00000302     1 22.079 0.00409 0.9496
## Light     0.00123193 0.00123193     1 22.058 1.66677 0.2101
## CO2:Light 0.00013157 0.00013157     1 22.056 0.17801 0.6772

Symbiodinium cells cm-2

############################################  ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##           Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## CO2        17978   17978     1    22 0.96158 0.3375
## Light      47910   47910     1    22 2.56251 0.1237
## CO2:Light   9532    9532     1    22 0.50984 0.4827

chlorophyll a cm-2

############################################ ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##            Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## CO2        0.0001  0.0001     1 21.368  0.0001    0.9932    
## Light     25.0845 25.0845     1 21.331 31.0554 1.486e-05 ***
## CO2:Light  1.8595  1.8595     1 21.328  2.3021    0.1439    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

chlorophyll c2 cm-2

############################################ ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##             Sum Sq  Mean Sq NumDF  DenDF F.value    Pr(>F)    
## CO2       0.000013 0.000013     1 21.796  0.0024 0.9614146    
## Light     0.103535 0.103535     1 21.751 18.8935 0.0002645 ***
## CO2:Light 0.014857 0.014857     1 21.747  2.7112 0.1140191    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

chlorophyll a cell-1

############################################  ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##              Sum Sq   Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       0.0000755 0.0000755     1 21.788 0.02109 0.8859
## Light     0.0099692 0.0099692     1 21.788 2.78608 0.1094
## CO2:Light 0.0010899 0.0010899     1 21.855 0.30459 0.5866

chlorophyll c2 cell-1

############################################  ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##              Sum Sq   Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       0.0004362 0.0004362     1 21.677 0.09203 0.7645
## Light     0.0076670 0.0076670     1 21.677 1.61758 0.2169
## CO2:Light 0.0028152 0.0028152     1 21.794 0.59394 0.4492

protein cm-2

############################################  ** 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
##  Light CO2     lsmean         SE    df  lower.CL  upper.CL .group
##  HL    HCO2 0.8724780 0.01614245 28.68 0.8394469 0.9055092  a    
##  HL    ACO2 0.9149473 0.01357862 27.81 0.8871243 0.9427704  ab   
##  LL    ACO2 0.9248994 0.01372553 28.74 0.8968165 0.9529822  ab   
##  LL    HCO2 0.9353848 0.01365200 28.27 0.9074318 0.9633379   b   
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
R8.table
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##             Sum Sq  Mean Sq NumDF  DenDF F.value  Pr(>F)  
## CO2       0.003914 0.003914     1 22.107  1.2788 0.27025  
## Light     0.024301 0.024301     1 22.056  7.9396 0.01001 *
## CO2:Light 0.014844 0.014844     1 22.020  4.8496 0.03842 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

carbohydrate cm-2

############################################   ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##              Sum Sq   Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       0.0017231 0.0017231     1 22.136 0.90856 0.3508
## Light     0.0008696 0.0008696     1 22.105 0.45854 0.5053
## CO2:Light 0.0046523 0.0046523     1 22.101 2.45307 0.1315

lipids cm-2

############################################  ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##            Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## CO2       0.40162 0.40162     1    22 2.28340 0.1450
## Light     0.01821 0.01821     1    22 0.10353 0.7507
## CO2:Light 0.03008 0.03008     1    22 0.17102 0.6832

energy content cm-2

############################################  ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##              Sum Sq   Mean Sq NumDF DenDF F.value Pr(>F)
## CO2       0.0097418 0.0097418     1    22 1.74027 0.2007
## Light     0.0021385 0.0021385     1    22 0.38203 0.5429
## CO2:Light 0.0001072 0.0001072     1    22 0.01915 0.8912

Biomass-normalized models

calcification gdw-1

## **************************************##
## 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##            Sum Sq Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       0.23537 0.23537     1 22.314 0.30561 0.5859
## Light     0.95754 0.95754     1 22.226 1.24331 0.2767
## CO2:Light 0.01836 0.01836     1 22.214 0.02384 0.8787

protein gdw-1

############################################  ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##               Sum Sq    Mean Sq NumDF  DenDF F.value Pr(>F)
## CO2       6.7983e-05 6.7983e-05     1 21.977 0.69180 0.4145
## Light     1.4785e-05 1.4785e-05     1 21.948 0.15046 0.7018
## CO2:Light 4.7995e-05 4.7995e-05     1 21.946 0.48840 0.4920

carbohydrates gdw-1

############################################  ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##               Sum Sq    Mean Sq NumDF  DenDF F.value  Pr(>F)  
## CO2       0.00019933 0.00019933     1 21.991  0.9431 0.34204  
## Light     0.00100331 0.00100331     1 21.953  4.7469 0.04037 *
## CO2:Light 0.00024106 0.00024106     1 21.950  1.1405 0.29715  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lipids gdw-1

############################################  ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##             Sum Sq  Mean Sq NumDF DenDF F.value  Pr(>F)  
## CO2       0.051069 0.051069     1    22  4.7616 0.04007 *
## Light     0.006752 0.006752     1    22  0.6295 0.43600  
## CO2:Light 0.000120 0.000120     1    22  0.0112 0.91663  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

energy content gdw-1

############################################   ** 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
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##           Sum Sq Mean Sq NumDF DenDF F.value  Pr(>F)  
## CO2       80.984  80.984     1    22  4.7211 0.04084 *
## Light     10.385  10.385     1    22  0.6054 0.44482  
## CO2:Light  0.062   0.062     1    22  0.0036 0.95247  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Figures

To generate plots, summary tables of mean, SE, and n were generated. Plots below show response variable mean +/- SE for corals exposed to CO2-Light treatments. Responses are normalized to area (cm2) or biomass (gdw).

Figures here show calcification (surface area (cm2) and biomass normalized (gdw)), symbiont densities (cells/cm2), tissue biomass (mg/cm2), chlorophyll a and c2 (ug/cm2), and chlorophyll a and c2 per symbiont cell (pg/cell). These figures show data found in “Figure 2”.

#####################################
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")

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”.

# 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")

Figures here show area-normalized energy reserves (mg/cm2): protein (mg/cm2), carbohydrates (mg/cm2), lipids (mg/cm2), and tissue energy content (kJ/cm2) These figures show data found in “Supplemental Figure 1”.

#####################################
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")