Evaluation of laser scanning confocal microscopy as a method for characterizing reef-building coral tissue thickness and Symbiodiniaceae fluorescence
AS Huffmyer, SB Matsuda, AR Eggers, JD Lemus, RD Gates
2020
Setup
Set up workspace, set options, and load required packages.
rm(list=ls(all=TRUE)) knitr::opts_chunk$set(root.dir = "~/Huffmyer_Confocal_Data",warning=FALSE, message=FALSE)#load packages
library("plotrix")
library("FSA")
library("ggplot2")
library("reshape2")
library("dplyr")
library("plyr")
library("pscl")
library("gridExtra")
library("car")
library("lsmeans")
library("nlme")
library("effects")
library("multcomp")
library("multcompView")
library("vegan")
library("exactRankTests")
library("mvtnorm")
library("PMCMR")
library("lmtest")
library("MuMIn")
library("glmmTMB")
library("lme4")
library("lmerTest")
library("HH")
library("FSA")
library("dunn.test")
library("emmeans")
library("tidyr")
library("ggsci")
library("cowplot")
library("vegan")
library("factoextra")
library("ggfortify")Analysis Overview
- Validation of measurements to compariative methodologies
- Examination of datasets and calculations
- Analysis of tissue thickness
- Analysis of symbiont densities and symbiont fluorescence
- Correlation of physiological characteristics
- Comparison of measurements over different tissue states
- Examination of datasets and calculations
- Utility of LSCM in documenting coral bleaching
- LSCM data collection: tissue thickness and symbiont fluorescence
- Correlation with comparative methodologies: chlorophyll content, tissue thickness, pigmentation
- LSCM data collection: tissue thickness and symbiont fluorescence
Validation of LSCM Methodology
1. Data Manipulation
A) Calculations
Raw data sets are inputted and manipulated to calculate the following in adult corals:
(1) Average cell counts on hemocytometer
(2) Average tissue thickness on dissecting scope
(3) Average tissue thickness on confocal
(4) Average symbiont fluorescence on confocal
#load data set for cell counts in adult tissues by hemocytometer
cells<-read.csv("Data/Validations/AdultCoral_SymbiontCellCounts.csv", header=TRUE, sep=",", na.strings="NA")
#load data set for thickness in adult tissues by a dissecting microscope
thickness<-read.csv("Data/Validations/AdultCoral_DissectScopeThickness.csv", header=TRUE, sep=",", na.strings="NA")
#load data set for LSCM symbiont fluoresence in adult tissues
fluorescence<-read.csv("Data/Validations/AdultCoral_ConfocalFluorescence.csv", header=TRUE, sep=",", na.strings="NA")
#load data set for LSCM thickness in adult tissues
zstack<-read.csv("Data/Validations/AdultCoral_ConfocalThickness.csv", header=TRUE, sep=",", na.strings="NA")
#summarize cell counts for adult tissues
#calculate average cell counts for each individual
totalcells <- aggregate(cells$Symbionts.cm2, by=list(cells$Sample, cells$Species), FUN=mean)
totalcells<-rename(totalcells, c("Group.1"="Sample", "Group.2"="Species", "x"="Density"))
#generate column to specify tissue type as "colony"
totalcells$Tissue<-c("Colony")
#summarize tissue thickness (via dissecting microscope) for adults
#calculate average tissue thickness for each individual
averagethickness <- aggregate(thickness$Thickness, by=list(thickness$Tunic), FUN=mean)
averagethickness<-rename(averagethickness, c("Group.1"="Sample", "x"="Thickness"))
#generate column to specify tissue type as "colony"
averagethickness$Tissue<-c("Colony")
#summarize symbiont fluorescence by LSCM for adults
#calculate average tissue intensity values
meantissue <- aggregate(fluorescence$Tissue.Intensity, by=list(fluorescence$Sample, fluorescence$Species), FUN=mean)
meantissue<-rename(meantissue, c("Group.1"="Sample", "Group.2"="Species", "x"="Tissue.Intensity"))
#calculate average symbiont fluorescence for adults
meansymbionts <- aggregate(fluorescence$Symbiont.Intensity, by=list(fluorescence$Sample, fluorescence$Species), FUN=mean)
meansymbionts<-rename(meansymbionts, c("Group.1"="Sample", "Group.2"="Species", "x"="Symbiont.Intensity"))
#merge symbionts and tissue dataframes together (note that tissue fluorescence is not analyzed here)
glow<-merge(meantissue, meansymbionts)
#summarize LSCM thickness for adults
meanthick <- aggregate(zstack$Thickness, by=list(zstack$Sample, zstack$Species), FUN=mean)
meanthick<-rename(meanthick, c("Group.1"="Sample", "Group.2"="Species", "x"="Zstack"))
#merge all dataframes together
master<-merge(totalcells, averagethickness)
master<-merge(master, glow)
master<-merge(master, meanthick)
#pull area data from cells dataframe
master$Area<-cells$AreaCM[match(master$Sample, cells$Sample)]
#write csv of master dataframe
write.csv(master, file="Output/AdultCoral_Master.csv")From the manipulated data frame, load the curated master data file. The master data frames will be used for all subsequent analyses.
B) Load master dataframes
Adult master dataframe
master<-read.csv("Output/AdultCoral_Master.csv", header=TRUE, sep=",", na.strings="NA")C) Calibrate LSCM Fluorescence Values
LSCM collects symbiont cell density data as the intensity of “red” (chlorophyll) fluorescence as specified in LSCM settings. This value can be calibrated to a % Relative Intensity value, by using the Red InSpeck fluorescent microbead kits. The intensity value calibration curves used in this study (specific settings are located in manuscript) are as follows:
| Species | Calibration Equation |
|---|---|
| M. capitata | (Red Intensity Value + 268.36) / 109.01 |
| P. compressa | (Red Intensity Value + 231.16) / 113.73 |
| P. varians | (Red Intensity Value + 32.77) / 41.087 |
| P. acuta | (Red Intensity Value + 464.76) / 193.16 |
Symbiont intensity values are calibrated with this equation to generate a % RI value and merged into the master dataframes for adult corals.
#Use conditional calculations by subsetting each species, then matching the values back into the master data sheets
mcap_cal<-subset(master, Species=="MCAP", select=c(Sample, Species, Symbiont.Intensity))
mcap_cal$Red<-(mcap_cal$Symbiont.Intensity+268.36)/109.01
pcom_cal<-subset(master, Species=="PCOM", select=c(Sample, Species, Symbiont.Intensity))
pcom_cal$Red<-(pcom_cal$Symbiont.Intensity+231.16)/113.73
pvar_cal<-subset(master, Species=="PVAR", select=c(Sample, Species, Symbiont.Intensity))
pvar_cal$Red<-(pvar_cal$Symbiont.Intensity+32.77)/41.087
pacu_cal<-subset(master, Species=="PACU", select=c(Sample, Species, Symbiont.Intensity))
pacu_cal$Red<-(pacu_cal$Symbiont.Intensity+464.76)/193.16
#merge all Red values into master data frame
calibrations<-bind_rows(mcap_cal, pcom_cal, pvar_cal, pacu_cal)
#combine into master, adult, juvies data sheets
master$Red<-calibrations$Red[match(master$Sample, calibrations$Sample)]D) Examining Data Structure
Adult Tissue Thickness: Dissecting Scope
qqPlot(master$Thickness)
## [1] 17 4Adult Tissue Thickness: LSCM
qqPlot(master$Zstack)
## [1] 8 4Adult Cell Density: Hemocytometer
qqPlot(master$Density)
## [1] 6 5Adult Cell Density: LSCM
qqPlot(master$Red)
## [1] 19 14All data for adult corals are normally distributed. All data are then subsetted by species for further analysis.
mcap<-subset(master, Species=="MCAP", c(Sample, Species, Tissue, Density, Thickness, Tissue.Intensity, Symbiont.Intensity, Zstack, Area, Red))
pcom<-subset(master, Species=="PCOM", c(Sample, Species, Tissue, Density, Thickness, Tissue.Intensity, Symbiont.Intensity, Zstack, Area, Red))
pvar<-subset(master, Species=="PVAR", c(Sample, Species, Tissue, Density, Thickness, Tissue.Intensity, Symbiont.Intensity, Zstack, Area, Red))
pacu<-subset(master, Species=="PACU", c(Sample, Species, Tissue, Density, Thickness, Tissue.Intensity, Symbiont.Intensity, Zstack, Area, Red))2. Analysis of Tissue Thickness
Tissue thickness of adult corals by two methods: dissecting scope and LSCM. All tissue thickness values will be analyzed by ANOVA analyses.
A) Tissue Thickness
Dissecting Microscope Method
model1a<-aov(Thickness~Species, data=master)
model1a<-aov(Thickness~Species, data=master)
summary(model1a)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 5665519 1888506 14.76 2.66e-05 ***
## Residuals 20 2558433 127922
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bartlett.test(Thickness~Species, data=master) #not violated##
## Bartlett test of homogeneity of variances
##
## data: Thickness by Species
## Bartlett's K-squared = 7.1466, df = 3, p-value = 0.06737test1a<-emmeans(model1a, "Species")
cld(test1a, Letter="abcdefg")## Species emmean SE df lower.CL upper.CL .group
## PACU 491 146 20 186 796 a
## PVAR 1313 146 20 1008 1617 b
## MCAP 1678 146 20 1373 1982 b
## PCOM 1684 146 20 1379 1988 b
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05TukeyHSD(model1a)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Thickness ~ Species, data = master)
##
## $Species
## diff lwr upr p adj
## PACU-MCAP -1186.753646 -1764.7225 -608.7848 0.0000697
## PCOM-MCAP 5.939604 -572.0292 583.9084 0.9999911
## PVAR-MCAP -365.141229 -943.1100 212.8276 0.3169675
## PCOM-PACU 1192.693250 614.7244 1770.6621 0.0000654
## PVAR-PACU 821.612417 243.6436 1399.5812 0.0037935
## PVAR-PCOM -371.080833 -949.0496 206.8880 0.3037791ANOVA analysis indicates a significant effect of species on tissue thickness measured by a dissecting microscope. The Bartlett’s test for homogeneity of variance is not violated.
LSCM Method model1b<-aov(Zstack~Species, data=master)
model1b<-aov(Zstack~Species, data=master)
summary(model1b)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 5527032 1842344 28.21 2.19e-07 ***
## Residuals 20 1306025 65301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bartlett.test(Zstack~Species, data=master) #not violated##
## Bartlett test of homogeneity of variances
##
## data: Zstack by Species
## Bartlett's K-squared = 4.0597, df = 3, p-value = 0.2551test1b<-emmeans(model1b, "Species")
cld(test1b, Letter="abcdefg")## Species emmean SE df lower.CL upper.CL .group
## PACU 496 104 20 278 714 a
## PVAR 1101 104 20 884 1319 b
## PCOM 1532 104 20 1314 1749 c
## MCAP 1754 104 20 1537 1972 c
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05TukeyHSD(model1b)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Zstack ~ Species, data = master)
##
## $Species
## diff lwr upr p adj
## PACU-MCAP -1258.5531 -1671.4990 -845.60729 0.0000002
## PCOM-MCAP -222.7537 -635.6996 190.19215 0.4504508
## PVAR-MCAP -653.0178 -1065.9636 -240.07192 0.0013687
## PCOM-PACU 1035.7994 622.8536 1448.74529 0.0000046
## PVAR-PACU 605.5354 192.5895 1018.48122 0.0028519
## PVAR-PCOM -430.2641 -843.2099 -17.31822 0.0392446ANOVA analysis indicates a significant effect of species on tissue thickness measured by LSCM. The Bartlett’s test for homogeneity of variance is not violated.
Analyze effect of methodology on adult sample thickness
#melt data frame to plot dissecting microscope scope and LSCM methods together
dfm <- melt(master[,c('Species','Thickness','Zstack')],id.vars = 1)
dfm$Group<-paste(dfm$Species,dfm$variable)model1c<-aov(value~variable, data=dfm), where value = thickness and variable = method.
model1c<-aov(value~variable, data=dfm)
summary(model1c)## Df Sum Sq Mean Sq F value Pr(>F)
## variable 1 59435 59435 0.182 0.672
## Residuals 46 15057009 327326model1c<-aov(value~Group, data=dfm)
test1c<-emmeans(model1c, "Group")
cld(test1c, Letter="abcdefg")## Group emmean SE df lower.CL upper.CL .group
## PACU Thickness 491 127 40 235 747 a
## PACU Zstack 496 127 40 239 752 a
## PVAR Zstack 1101 127 40 845 1358 b
## PVAR Thickness 1313 127 40 1056 1569 bc
## PCOM Zstack 1532 127 40 1275 1788 bc
## MCAP Thickness 1678 127 40 1421 1934 c
## PCOM Thickness 1684 127 40 1427 1940 c
## MCAP Zstack 1754 127 40 1498 2011 c
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 8 estimates
## significance level used: alpha = 0.05TukeyHSD(model1c)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = value ~ Group, data = dfm)
##
## $Group
## diff lwr upr p adj
## MCAP Zstack-MCAP Thickness 76.715338 -496.91031 650.340986 0.9998583
## PACU Thickness-MCAP Thickness -1186.753646 -1760.37929 -613.127998 0.0000017
## PACU Zstack-MCAP Thickness -1181.837804 -1755.46345 -608.212156 0.0000019
## PCOM Thickness-MCAP Thickness 5.939604 -567.68604 579.565252 1.0000000
## PCOM Zstack-MCAP Thickness -146.038366 -719.66401 427.587282 0.9913215
## PVAR Thickness-MCAP Thickness -365.141229 -938.76688 208.484419 0.4726364
## PVAR Zstack-MCAP Thickness -576.302440 -1149.92809 -2.676792 0.0482167
## PACU Thickness-MCAP Zstack -1263.468984 -1837.09463 -689.843336 0.0000004
## PACU Zstack-MCAP Zstack -1258.553142 -1832.17879 -684.927494 0.0000005
## PCOM Thickness-MCAP Zstack -70.775734 -644.40138 502.849914 0.9999175
## PCOM Zstack-MCAP Zstack -222.753704 -796.37935 350.871945 0.9142564
## PVAR Thickness-MCAP Zstack -441.856567 -1015.48222 131.769081 0.2409514
## PVAR Zstack-MCAP Zstack -653.017778 -1226.64343 -79.392130 0.0160724
## PACU Zstack-PACU Thickness 4.915842 -568.70981 578.541490 1.0000000
## PCOM Thickness-PACU Thickness 1192.693250 619.06760 1766.318898 0.0000016
## PCOM Zstack-PACU Thickness 1040.715280 467.08963 1614.340928 0.0000235
## PVAR Thickness-PACU Thickness 821.612417 247.98677 1395.238065 0.0010750
## PVAR Zstack-PACU Thickness 610.451206 36.82556 1184.076854 0.0299600
## PCOM Thickness-PACU Zstack 1187.777408 614.15176 1761.403056 0.0000017
## PCOM Zstack-PACU Zstack 1035.799438 462.17379 1609.425086 0.0000256
## PVAR Thickness-PACU Zstack 816.696575 243.07093 1390.322223 0.0011679
## PVAR Zstack-PACU Zstack 605.535364 31.90972 1179.161012 0.0321284
## PCOM Zstack-PCOM Thickness -151.977970 -725.60362 421.647678 0.9890164
## PVAR Thickness-PCOM Thickness -371.080833 -944.70648 202.544815 0.4520151
## PVAR Zstack-PCOM Thickness -582.242044 -1155.86769 -8.616396 0.0444602
## PVAR Thickness-PCOM Zstack -219.102863 -792.72851 354.522785 0.9208055
## PVAR Zstack-PCOM Zstack -430.264074 -1003.88972 143.361574 0.2702562
## PVAR Zstack-PVAR Thickness -211.161211 -784.78686 362.464438 0.9339213ANOVA analysis indicates that there is no significant effect of method on tissue thickness values in adult tissues.
Coefficient of variation (CV): Dissecting microscope and LSCM, respectively.
#variability on scope
sd(master$Thickness)/mean(master$Thickness)## [1] 0.4630851#variability on LSCM
sd(master$Zstack)/mean(master$Zstack)## [1] 0.4464446Thick1 <- ddply(dfm, c("Species", "variable"), summarise,
N = length(value[!is.na(value)]),
mean = mean(value, na.rm=TRUE),
sd = sd(value, na.rm=TRUE),
se = sd / sqrt(N),
max = max(value, na.rm=TRUE)
)
knitr::kable(Thick1)| Species | variable | N | mean | sd | se | max |
|---|---|---|---|---|---|---|
| MCAP | Thickness | 6 | 1677.7543 | 389.3004 | 158.93122 | 2270.9841 |
| MCAP | Zstack | 6 | 1754.4696 | 311.0328 | 126.97862 | 2096.9122 |
| PACU | Thickness | 6 | 491.0006 | 132.5836 | 54.12701 | 635.1969 |
| PACU | Zstack | 6 | 495.9165 | 120.4565 | 49.17616 | 659.5892 |
| PCOM | Thickness | 6 | 1683.6939 | 510.4328 | 208.38333 | 2282.3106 |
| PCOM | Zstack | 6 | 1531.7159 | 301.1120 | 122.92845 | 1855.2033 |
| PVAR | Thickness | 6 | 1312.6131 | 286.3769 | 116.91290 | 1603.2980 |
| PVAR | Zstack | 6 | 1101.4519 | 243.4859 | 99.40271 | 1465.7067 |
Generate a figure for display.
par(mar = c(1, 1, 1, 1))
ThickSpecies<-ggplot(Thick1, aes(x=Species, y=mean, fill=variable)) +
geom_bar(position=position_dodge(), stat="identity", color="black") +
scale_fill_manual(name="Method",
values=c("lightgray", "darkgray"),
labels=c("Dissecting \nMicroscope", "LSCM"))+
scale_x_discrete(limits=c("PACU", "PVAR", "MCAP","PCOM"), labels=c( "P. acuta", "P. varians", "M. capitata", "P. compressa"))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
width=0.1,
position=position_dodge(.9))+
theme_classic()+
ylim(0,2050)+
geom_text(x=0.78, y=625, label="a", size=8) + #pacuta dissecting
geom_text(x=1.22, y=625, label="a", size=8) + #pacuta confocal
geom_text(x=1.78, y=1525, label="bc", size=8) + #pvar dissecting
geom_text(x=2.22, y=1275, label="b", size=8) + #pvar confocal
geom_text(x=2.78, y=1925, label="c", size=8) + #mcap dissecting
geom_text(x=3.22, y=1975, label="c", size=8) + #mcap confocal
geom_text(x=3.78, y=1975, label="c", size=8) + #pcom dissecting
geom_text(x=4.22, y=1750, label="bc", size=8) + #pcom confocal
theme(legend.key.height=unit(1, "cm"))+
theme(text = element_text(size = 18))+
theme(axis.text = element_text(size = 18, color="black"))+
theme(axis.title = element_text(size = 18, face="bold"))+
theme(legend.title = element_text(size = 18, face="bold"))+
theme(legend.text = element_text(size = 18, color="black"))+
theme(axis.text.x = element_text(angle = 0, size=18, face="italic"))+
ylab(expression(bold(paste("Tissue Thickness ("*mu*"m)")))); ThickSpecies
ggsave("Figures/Fig1a.pdf", plot=ThickSpecies, height=8, width=10, units = c("in"), dpi=300) #output figureSummary - Adult Tissue Thickness
(1) There was a significant effect of species on tissue thickness as measured by both the dissecting microscope and LSCM methods.
(2) LSCM detected differences between species that was not distinguished using the dissecting microscope.
(3) There was no difference in tissue thickness using the two methods, indicating the LSCM method is valid.
B) Correlation of Methodologies
Correlations between LSCM and dissecting microscope measurements in adult coral tissue
Conduct Pearson correlation test between LSCM and dissecting microscope measurements of tissue thickness in adult tissues.
cor.test(master$Thickness, master$Zstack, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: master$Thickness and master$Zstack
## t = 6.6603, df = 22, p-value = 1.073e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6180457 0.9181690
## sample estimates:
## cor
## 0.8176003ThicknessCorrelation<-ggplot(master, aes(x=Thickness, y=Zstack)) +
geom_point(aes(color=Species), size=3) +
scale_color_manual(name="Species",
values=c("royalblue4", "red4", "darkgreen", "darkorchid4"),
labels=c("M. capitata", "P. acuta", "P. compressa", "P. varians"))+
geom_smooth(method=lm, color="black", se=FALSE)+
theme_classic()+
geom_text(x=2000, y=500, label="r=0.82, p<0.01", size=6, fontface="plain") + #correlation results
theme(axis.text = element_text(size = 18, color="black"))+
theme(axis.title = element_text(size = 18, face="bold"))+
theme(legend.title = element_text(size = 18, face="bold"))+
theme(legend.text = element_text(size = 18, color="black", face="italic"))+
ylab(expression(bold(paste("Depth of Fluorescence ("*mu*"m)")))) +
xlab(expression(bold(paste("Tissue Thickness ("*mu*"m)"))))
ThicknessCorrelation
ggsave("Figures/Fig1b.pdf", plot=ThicknessCorrelation, height=8, width=10, units = c("in"), dpi=300) #output figure3. Analysis of Symbiont Densities
Analysis of symbiont densities (cells per cm^2) as compared to symbiont fluorescence (%RI) by LSCM. All data analyzed using ANOVA analyses.
A) Symbiont Cell Densities - Cell Counts via Hemocytometer
Adult tissue hemocytometer symbiont cell densiites
model3a<-aov(Density~Species, data=master)
model3a<-aov(Density~Species, data=master)
summary(model3a)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 1.052e+12 3.507e+11 6.717 0.00257 **
## Residuals 20 1.044e+12 5.221e+10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bartlett.test(Density~Species, data=master) #not violated##
## Bartlett test of homogeneity of variances
##
## data: Density by Species
## Bartlett's K-squared = 4.9297, df = 3, p-value = 0.177test3a<-emmeans(model3a, "Species")
cld(test3a, Letter="abcdefg")## Species emmean SE df lower.CL upper.CL .group
## PACU 477615 93284 20 283028 672202 a
## PCOM 638188 93284 20 443602 832775 ab
## MCAP 905119 93284 20 710532 1099706 b
## PVAR 1004476 93284 20 809890 1199063 b
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05TukeyHSD(model3a)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Density ~ Species, data = master)
##
## $Species
## diff lwr upr p adj
## PACU-MCAP -427503.48 -796748.643 -58258.31 0.0196794
## PCOM-MCAP -266930.27 -636175.432 102314.90 0.2128283
## PVAR-MCAP 99357.68 -269887.486 468602.85 0.8742532
## PCOM-PACU 160573.21 -208671.956 529818.38 0.6235342
## PVAR-PACU 526861.16 157615.990 896106.32 0.0036677
## PVAR-PCOM 366287.95 -2957.221 735533.11 0.0523387ANOVA analysis indicate that there is a significant difference in symbiont cell densities between species as measured by hemocytometer. Bartlett’s test is not violated.
Cells1 <- ddply(master, c("Species"), summarise,
N = length(Density[!is.na(Density)]),
mean = mean(Density, na.rm=TRUE),
sd = sd(Density, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Density, na.rm=TRUE)
)
knitr::kable(Cells1)| Species | N | mean | sd | se | max |
|---|---|---|---|---|---|
| MCAP | 6 | 905118.7 | 355696.9 | 145212.65 | 1411016.9 |
| PACU | 6 | 477615.3 | 161399.4 | 65891.03 | 706967.2 |
| PCOM | 6 | 638188.5 | 179566.1 | 73307.57 | 901408.5 |
| PVAR | 6 | 1004476.4 | 155020.5 | 63286.84 | 1166163.1 |
Generate figure.
#colors are "royalblue4", "red4", "darkgreen", "darkorchid4"
AdultDensities<-ggplot(Cells1, aes(x=Species, y=mean, fill=Species)) +
geom_bar(position=position_dodge(), stat="identity", colour="black") +
scale_fill_manual(name="Species",
values=c("lightgray", "lightgray", "lightgray", "lightgray"))+
scale_x_discrete(limits=c("PACU", "PVAR", "MCAP", "PCOM"), labels=c("P. acuta", "P. varians", "M. capitata", "P. compressa"))+
scale_y_continuous(labels = scales::comma, limit=c(0,1250000))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
width=.1,
position=position_dodge(.9))+
geom_text(x=1, y=580000, label="a", size=6) + #pacuta
geom_text(x=2, y=1110000, label="b", size=6) + #pvarians
geom_text(x=3, y=1090000, label="b", size=6) + #mcapitata
geom_text(x=4, y=750000, label="ab", size=6) + #pcompressa
theme_classic()+
theme(text = element_text(size = 18))+
theme(axis.text = element_text(size = 18, color="black"))+
theme(axis.text.x = element_text(size = 18, color="black", face="italic"))+
theme(axis.title = element_text(size = 18, face="bold"))+
theme(legend.position = "none")+
theme(legend.text = element_text(size = 18, color="black"))+
ylab(expression(bold(paste("Symbiodiniaceae Cells cm"^-2))))
AdultDensities
ggsave("Figures/Fig1d.pdf", plot=AdultDensities, height=8, width=8, units = c("in"), dpi=300) #output figureB) Symbiont Fluorescence - LSCM
Adult tissue LSCM symbiont fluorescence (%RI)
model3b<-aov(Red~Species, data=adults)
model3b<-aov(Red~Species, data=master)
summary(model3b)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 1819 606.3 5.979 0.00442 **
## Residuals 20 2028 101.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bartlett.test(Red~Species, data=master) #not violated##
## Bartlett test of homogeneity of variances
##
## data: Red by Species
## Bartlett's K-squared = 0.25907, df = 3, p-value = 0.9675test3b<-emmeans(model3b, "Species")
cld(test3b, Letter="abcdefg")## Species emmean SE df lower.CL upper.CL .group
## PCOM 27.9 4.11 20 19.3 36.5 a
## PACU 32.8 4.11 20 24.3 41.4 ab
## MCAP 45.6 4.11 20 37.0 54.1 b
## PVAR 49.0 4.11 20 40.4 57.6 b
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05TukeyHSD(model3b)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Red ~ Species, data = master)
##
## $Species
## diff lwr upr p adj
## PACU-MCAP -12.716954 -28.9889388 3.55503 0.1609121
## PCOM-MCAP -17.641685 -33.9136692 -1.36970 0.0306139
## PVAR-MCAP 3.415138 -12.8568464 19.68712 0.9347125
## PCOM-PACU -4.924730 -21.1967149 11.34725 0.8313971
## PVAR-PACU 16.132092 -0.1398921 32.40408 0.0525143
## PVAR-PCOM 21.056823 4.7848382 37.32881 0.0084836ANOVA analysis indicates that symbiont fluorescence (% RI) is significantly different between species. Barlett’s test is not violated.
Red1 <- ddply(master, c("Species"), summarise,
N = length(Red[!is.na(Red)]),
mean = mean(Red, na.rm=TRUE),
sd = sd(Red, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Red, na.rm=TRUE)
)
knitr::kable(Red1)| Species | N | mean | sd | se | max |
|---|---|---|---|---|---|
| MCAP | 6 | 45.56163 | 9.499923 | 3.878327 | 58.53713 |
| PACU | 6 | 32.84467 | 8.970417 | 3.662158 | 44.45967 |
| PCOM | 6 | 27.91994 | 11.086587 | 4.526080 | 41.72114 |
| PVAR | 6 | 48.97676 | 10.580664 | 4.319538 | 63.00588 |
Generate figure.
AdultRed<-ggplot(Red1, aes(x=Species, y=mean, fill=Species)) +
geom_bar(position=position_dodge(), stat="identity", colour="black") +
scale_fill_manual(name="Species",
values=c("darkgray", "darkgray", "darkgray", "darkgray"))+
scale_x_discrete(limits=c("PACU", "PVAR", "MCAP", "PCOM"), labels=c("P. acuta", "P. varians", "M. capitata", "P. compressa"))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
width=.1,
position=position_dodge(.9))+
ylim(0,65)+
geom_text(x=1, y=39, label="ab", size=6) + #pacuta
geom_text(x=2, y=55.5, label="b", size=6) + #pvarians
geom_text(x=3, y=52, label="b", size=6) + #mcapitata
geom_text(x=4, y=35, label="a", size=6) + #pcompressa
theme_classic()+
theme(text = element_text(size = 18))+
theme(axis.text = element_text(size = 18, color="black"))+
theme(axis.text.x = element_text(size = 18, color="black", face="italic"))+
theme(axis.title = element_text(size = 18, face="bold"))+
theme(legend.position = "none")+
theme(legend.text = element_text(size = 18, color="black"))+
ylab(expression(bold(paste("Symbiodiniaceae Fluorescence (% RI)"))))
AdultRed
ggsave("Figures/Fig1c.pdf", plot=AdultRed, height=8, width=8, units = c("in"), dpi=300) #output figureC) Correlation of methodologies
Coefficient of variation (CV) of adult tissue cell densities on hemocytometer, fluorescence on LSCM, and chlorophyll content respectively:
sd(master$Density)/mean(master$Density)## [1] 0.3991586sd(master$Red)/mean(master$Red)## [1] 0.3330869Correlation of Red %RI and symbiont cell densities in adult corals
cor.test(master$Density, master$Red, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: master$Density and master$Red
## t = 3.0314, df = 22, p-value = 0.006131
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1784844 0.7762301
## sample estimates:
## cor
## 0.5427995CellsCorrelation<-ggplot(master, aes(x=Density, y=Red)) +
geom_point(aes(color=Species), size=3) +
geom_smooth(color="black", method=lm, se=FALSE)+
scale_color_manual(name="Species",
values=c("royalblue4", "red4", "darkgreen", "darkorchid4"),
labels=c("M. capitata", "P. acuta", "P. compressa", "P. varians"))+
theme_classic()+
scale_x_continuous(labels = scales::comma)+
geom_text(x=1200000, y=20, label="r=0.54, p<0.01", size=6) + #results
theme(text = element_text(size = 18))+
theme(axis.text = element_text(size = 18, color="black"))+
theme(axis.title = element_text(size = 18, face="bold"))+
theme(legend.title = element_text(size = 18, face="bold"))+
theme(legend.text = element_text(size = 18, color="black", face="italic"))+
ylab("Symbiodiniaceae Fluorescence (% RI)") +
xlab(expression(bold(paste("Symbiodiniaceae Cells cm"^-2))))
CellsCorrelation
ggsave("Figures/Fig1e.pdf", plot=CellsCorrelation, height=8, width=10, units = c("in"), dpi=300) #output figureThere is a significant correlation between cell density and LSCM red values.
Summary - Symbiont Densities compared to LSCM
- There was a significant effect of species on symbiont densities as measured by both the dissecting microscope and LSCM methods.
- Trends between species were similar between hemocytometer and LSCM measurements and were significantly correlated in adult samples.
- Low r value in correlation likely reflects decrease in chlorophyll fluorescence following decalcification and fixation.
4. Coefficients of variation
Summarize variation of each metric for each species.
CV1 <- ddply(master, c("Species"), summarise,
meanRed = mean(Red, na.rm=TRUE),
sdRed = sd(Red, na.rm=TRUE),
CVRed = (sdRed/meanRed),
meanDensity = mean(Density, na.rm=TRUE),
sdDensity = sd(Density, na.rm=TRUE),
CVDensity = (sdDensity/meanDensity),
meanZstack = mean(Zstack, na.rm=TRUE),
sdZstack = sd(Zstack, na.rm=TRUE),
CVZstack = (sdZstack/meanZstack),
meanThickness = mean(Thickness, na.rm=TRUE),
sdThickness = sd(Thickness, na.rm=TRUE),
CVThickness = (sdThickness/meanThickness)
)
knitr::kable(CV1)| Species | meanRed | sdRed | CVRed | meanDensity | sdDensity | CVDensity | meanZstack | sdZstack | CVZstack | meanThickness | sdThickness | CVThickness |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| MCAP | 45.56163 | 9.499923 | 0.2085071 | 905118.7 | 355696.9 | 0.3929837 | 1754.4696 | 311.0328 | 0.1772803 | 1677.7543 | 389.3004 | 0.2320366 |
| PACU | 32.84467 | 8.970417 | 0.2731164 | 477615.3 | 161399.4 | 0.3379276 | 495.9165 | 120.4565 | 0.2428968 | 491.0006 | 132.5836 | 0.2700273 |
| PCOM | 27.91994 | 11.086587 | 0.3970849 | 638188.5 | 179566.1 | 0.2813685 | 1531.7159 | 301.1120 | 0.1965847 | 1683.6939 | 510.4328 | 0.3031625 |
| PVAR | 48.97676 | 10.580664 | 0.2160344 | 1004476.4 | 155020.5 | 0.1543296 | 1101.4519 | 243.4859 | 0.2210591 | 1312.6131 | 286.3769 | 0.2181732 |
Bleaching Experiment
Juvenile and adult corals of M. capitata and P. acuta were exposed to a thermal stress to test the utility of LSCM in documenting physiological characteristics across enviornmental conditions. Here, we measure tissue thickness and symbiont fluorescence using LSCM at beginning and end points and destructively sample corals for tissue thickness, symbiont cell density and symbiont fluorescence at the end point. We display tissue thickness and symbiont fluorescence over time using LSCM measurements and calculate these metrics in ambient and high temperature groups at the final time point.
The purpose of this experiment is to test whether LSCM is 1) valid in use with living corals across physiological conditions; 2) investigate the correlation of LSCM tissue thickness and tissue thickness on the dissecting microsope; and 3) investigate the correlation of LSCM symbiont fluorescence with chlorophyll content and coral pigmentation.
As confocal microscopy allows for measurement on living corals, data is available at initial and end points of bleaching experiment. Destructive measures (chlorophyll a content, tissue thickness on dissecting scope) as well as comparative measurement of coral pigmentation collected by photographs was collected at the endpoint of the experiment.
1. LSCM Tissue Thickness
Load dataset and calculate mean tissue thickness values for each individual coral at start and ending timepoints collected by LSCM.
#load dataset
bleach_thick<-read.csv("Data/Bleaching/Bleaching_Confocal_Thickness.csv", header=TRUE, sep=",", na.strings="NA")
#calculate mean for each sample
bleach1 <- ddply(bleach_thick, c("Timepoint", "Sample", "Treatment", "Species", "LifeStage", "Colony"), summarise,
Thickness = mean(Thickness, na.rm=TRUE)
)
#summarize by species at each timepoint and treatment
bleachSum <- ddply(bleach1, c("Timepoint", "Treatment", "Species", "LifeStage"), summarise,
N = length(Thickness[!is.na(Thickness)]),
mean = mean(Thickness, na.rm=TRUE),
sd = sd(Thickness, na.rm=TRUE),
se = sd / sqrt(N)
)
knitr::kable(bleachSum)
#set the intial timepoint as the reference level
bleachSum$Timepoint <- relevel(bleachSum$Timepoint, ref = "Initial")
#set the intial timepoint as the reference level
bleach1$Timepoint <- relevel(bleach1$Timepoint, ref = "Initial")
#subset by species
bleach1$Group<-paste(bleach1$LifeStage, bleach1$Treatment)
bleach_mcap<-subset(bleach1, Species=="Montipora capitata", c(Treatment, Timepoint, Species, LifeStage, Sample, Thickness, Group, Colony))
bleach_pacu<-subset(bleach1, Species=="Pocillopora acuta", c(Treatment, Timepoint, Species, LifeStage, Sample, Thickness, Group, Colony))Display a figure of mean tissue thickness over the course of bleaching with a separate plot for each species.
FigBleach1 <- ggplot(data=bleachSum, aes(x=Timepoint, y=mean, group=interaction(Treatment, LifeStage), colour=Treatment, linetype=LifeStage)) +
facet_wrap(~Species, scale="free")+
geom_line(position=position_dodge(0.2), size=1.3) +
scale_linetype_manual(name = "Life Stage", values=c("solid", "dotted"))+
scale_color_manual(name = "Temperature", values=c("blue", "red")) +
geom_point(size=3, position=position_dodge(0.2)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0, linetype="solid", position=position_dodge(0.2), size=1.3) +
xlab("Timepoint") +
ylab(expression(bold(paste("Depth of Fluorescence (", mu, "m)")))) +
ylim(0,2000) +
theme_bw()+
theme_classic() +
theme(panel.background = element_blank()) +
theme(strip.text.x = element_text(face="italic"))+
theme(text = element_text(size=20),
panel.border = element_blank(),
plot.background =element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_text(size=18, colour="black"),
legend.title=element_text(face="bold"),
legend.position="right",
plot.title=element_text(),
legend.text = element_text(size = 20),
axis.title.y=element_text(margin=margin(), face="bold"),
axis.title.x=element_text(margin=margin(), face="bold"))
FigBleach1
ggsave(filename="Figures/Fig2a.pdf", plot=FigBleach1, dpi=300, width=10, height=7, units="in")Analyze the differences in tissue thickness between timepoints between temperature and life stage with a separate repeated measures linear mixed effect model for each species.
Montipora capitata:
modelB1a<-lmer(Thickness~LifeStage * Treatment * Timepoint + (1|Colony) + (1|Sample), data=bleach_mcap)
summary(modelB1a)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Thickness ~ LifeStage * Treatment * Timepoint + (1 | Colony) +
## (1 | Sample)
## Data: bleach_mcap
##
## REML criterion at convergence: 724.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07230 -0.60093 -0.05217 0.61615 2.21473
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample (Intercept) 6.280e+03 7.924e+01
## Colony (Intercept) 3.954e-05 6.288e-03
## Residual 5.772e+04 2.402e+02
## Number of obs: 59, groups: Sample, 30; Colony, 16
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 1783.95 103.28 50.53
## LifeStageJuvenile -234.80 136.63 50.53
## TreatmentHigh -150.95 140.75 50.53
## TimepointFinal -352.17 138.71 25.46
## LifeStageJuvenile:TreatmentHigh 33.45 186.87 50.53
## LifeStageJuvenile:TimepointFinal 146.95 186.55 26.22
## TreatmentHigh:TimepointFinal 34.67 189.03 25.46
## LifeStageJuvenile:TreatmentHigh:TimepointFinal -34.82 253.22 25.87
## t value Pr(>|t|)
## (Intercept) 17.273 <2e-16 ***
## LifeStageJuvenile -1.719 0.0918 .
## TreatmentHigh -1.072 0.2886
## TimepointFinal -2.539 0.0176 *
## LifeStageJuvenile:TreatmentHigh 0.179 0.8586
## LifeStageJuvenile:TimepointFinal 0.788 0.4379
## TreatmentHigh:TimepointFinal 0.183 0.8559
## LifeStageJuvenile:TreatmentHigh:TimepointFinal -0.138 0.8917
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) LfStgJ TrtmnH TmpntF LfSJ:TH LSJ:TF TrH:TF
## LifeStgJvnl -0.756
## TreatmntHgh -0.734 0.555
## TimepontFnl -0.672 0.508 0.493
## LfStgJvn:TH 0.553 -0.731 -0.753 -0.371
## LfStgJvn:TF 0.499 -0.661 -0.366 -0.744 0.483
## TrtmntHg:TF 0.493 -0.372 -0.672 -0.734 0.506 0.546
## LfStJ:TH:TF -0.368 0.487 0.501 0.548 -0.666 -0.737 -0.746
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(modelB1a, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## LifeStage 282119 282119 1 25.907 4.8877 0.0360828
## Treatment 187113 187113 1 26.078 3.2417 0.0833625
## Timepoint 1008578 1008578 1 25.893 17.4737 0.0002938
## LifeStage:Treatment 791 791 1 25.972 0.0137 0.9077271
## LifeStage:Timepoint 59503 59503 1 25.810 1.0309 0.3193735
## Treatment:Timepoint 847 847 1 25.982 0.0147 0.9044969
## LifeStage:Treatment:Timepoint 1091 1091 1 25.873 0.0189 0.8916907
##
## LifeStage *
## Treatment .
## Timepoint ***
## LifeStage:Treatment
## LifeStage:Timepoint
## Treatment:Timepoint
## LifeStage:Treatment:Timepoint
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm = emmeans(modelB1a, ~ LifeStage, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## LifeStage emmean SE df lower.CL upper.CL .group
## Juvenile 1388 46.7 13.5 1287 1488 A
## Adult 1541 52.5 13.1 1428 1654 B
##
## Results are averaged over the levels of: Treatment, Timepoint
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Adult - Juvenile 153 70.2 13.3 2.183 0.0476
##
## Results are averaged over the levels of: Treatment, Timepoint
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB1a, ~ Treatment, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## Treatment emmean SE df lower.CL upper.CL .group
## High 1402 47.2 25.3 1304 1499 A
## Ambient 1527 52.0 26.1 1420 1634 A
##
## Results are averaged over the levels of: LifeStage, Timepoint
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Ambient - High 126 70.2 13.3 1.788 0.0966
##
## Results are averaged over the levels of: LifeStage, Timepoint
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB1a, ~ Timepoint, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## Timepoint emmean SE df lower.CL upper.CL .group
## Final 1329 47.7 32.9 1232 1427 A
## Initial 1599 46.9 32.1 1504 1695 B
##
## Results are averaged over the levels of: LifeStage, Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Initial - Final 270 63.4 25.6 4.260 0.0002
##
## Results are averaged over the levels of: LifeStage, Treatment
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB1a, ~ LifeStage * Treatment *Timepoint, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## LifeStage Treatment Timepoint emmean SE df lower.CL upper.CL .group
## Juvenile High Final 1226 84.3 50.5 1057 1396 A
## Adult High Final 1316 95.6 50.5 1123 1508 AB
## Juvenile Ambient Final 1344 96.9 50.7 1149 1539 ABC
## Juvenile High Initial 1432 84.3 50.5 1262 1601 ABC
## Adult Ambient Final 1432 104.0 50.5 1223 1641 ABC
## Juvenile Ambient Initial 1549 89.9 50.5 1369 1730 ABC
## Adult High Initial 1633 95.6 50.5 1441 1825 BC
## Adult Ambient Initial 1784 104.0 50.5 1575 1993 C
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 8 estimates
## significance level used: alpha = 0.05Assumptions of analysis are met.
qqPlot(residuals(modelB1a))
## 41 98
## 19 48bartlett.test(residuals(modelB1a)~bleach_mcap$Group) #not violated##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelB1a) by bleach_mcap$Group
## Bartlett's K-squared = 1.0933, df = 3, p-value = 0.7787There was a significant effect of life stage and timepoint in Montipora capitata, but no effect of temperature treatment.
Pocillopora acuta:
modelB1b<-lmer(Thickness~LifeStage * Treatment * Timepoint + (1|Sample) + (1|Colony), data=bleach_pacu)
summary(modelB1b)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Thickness ~ LifeStage * Treatment * Timepoint + (1 | Sample) +
## (1 | Colony)
## Data: bleach_pacu
##
## REML criterion at convergence: 615
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1157 -0.4727 0.0667 0.5520 1.6198
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample (Intercept) 4335.5 65.84
## Colony (Intercept) 448.4 21.17
## Residual 11468.2 107.09
## Number of obs: 56, groups: Sample, 28; Colony, 14
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 742.25 48.18 44.11
## LifeStageJuvenile 219.35 68.14 44.11
## TreatmentHigh 23.92 67.20 25.43
## TimepointFinal 42.48 57.24 24.00
## LifeStageJuvenile:TreatmentHigh -118.30 95.03 25.43
## LifeStageJuvenile:TimepointFinal -177.95 80.95 24.00
## TreatmentHigh:TimepointFinal -284.39 80.95 24.00
## LifeStageJuvenile:TreatmentHigh:TimepointFinal 256.38 114.48 24.00
## t value Pr(>|t|)
## (Intercept) 15.404 < 2e-16 ***
## LifeStageJuvenile 3.219 0.00241 **
## TreatmentHigh 0.356 0.72480
## TimepointFinal 0.742 0.46526
## LifeStageJuvenile:TreatmentHigh -1.245 0.22451
## LifeStageJuvenile:TimepointFinal -2.198 0.03782 *
## TreatmentHigh:TimepointFinal -3.513 0.00178 **
## LifeStageJuvenile:TreatmentHigh:TimepointFinal 2.239 0.03466 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) LfStgJ TrtmnH TmpntF LfSJ:TH LSJ:TF TrH:TF
## LifeStgJvnl -0.707
## TreatmntHgh -0.697 0.493
## TimepontFnl -0.594 0.420 0.426
## LfStgJvn:TH 0.493 -0.697 -0.707 -0.301
## LfStgJvn:TF 0.420 -0.594 -0.301 -0.707 0.426
## TrtmntHg:TF 0.420 -0.297 -0.602 -0.707 0.426 0.500
## LfStJ:TH:TF -0.297 0.420 0.426 0.500 -0.602 -0.707 -0.707anova(modelB1b, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## LifeStage 134047 134047 1 11.999 11.6886 0.005090 **
## Treatment 102396 102396 1 11.999 8.9287 0.011315 *
## Timepoint 217351 217351 1 23.999 18.9525 0.000215 ***
## LifeStage:Treatment 195 195 1 11.999 0.0170 0.898463
## LifeStage:Timepoint 8667 8667 1 23.999 0.7557 0.393277
## Treatment:Timepoint 85394 85394 1 23.999 7.4462 0.011707 *
## LifeStage:Treatment:Timepoint 57515 57515 1 23.999 5.0151 0.034658 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm = emmeans(modelB1b, ~ LifeStage, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## LifeStage emmean SE df lower.CL upper.CL .group
## Adult 704 28 12 643 765 A
## Juvenile 840 28 12 779 901 B
##
## Results are averaged over the levels of: Treatment, Timepoint
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Adult - Juvenile -135 39.6 12 -3.419 0.0051
##
## Results are averaged over the levels of: Treatment, Timepoint
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB1b, ~ Treatment, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## Treatment emmean SE df lower.CL upper.CL .group
## High 715 27.4 24 659 772 A
## Ambient 829 27.4 24 772 885 B
##
## Results are averaged over the levels of: LifeStage, Timepoint
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Ambient - High 113 37.9 12 2.988 0.0113
##
## Results are averaged over the levels of: LifeStage, Timepoint
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB1b, ~ Timepoint, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## Timepoint emmean SE df lower.CL upper.CL .group
## Final 710 24.4 24.5 659 760 A
## Initial 834 24.4 24.5 784 885 B
##
## Results are averaged over the levels of: LifeStage, Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Initial - Final 125 28.6 24 4.353 0.0002
##
## Results are averaged over the levels of: LifeStage, Treatment
## Degrees-of-freedom method: kenward-rogerAssumptions of analysis are met.
qqPlot(residuals(modelB1b))
## 29 57
## 12 28bartlett.test(residuals(modelB1b)~bleach_pacu$Group) #not violated##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelB1b) by bleach_pacu$Group
## Bartlett's K-squared = 1.0457, df = 3, p-value = 0.7902There was a significant effect of treatment, timepoint, interaction between treatment and timepiont, and an interaction between lifestage, treatment, and timepoint on Pocillopora acuta tissue thickness. This analysis indicates that there was a negative effect of high temperature by the final timpoint. The three way interaction suggests that life stages responded differently (+ effect of juvenile stage) to high temperature treatment at the final timepoint.
Coefficient of variation (CV):
sd(bleach1$Thickness)/mean(bleach1$Thickness)## [1] 0.3715388CVa <- ddply(bleach1, c("Treatment", "Species", "LifeStage"), summarise,
mean = mean(Thickness, na.rm=TRUE),
sd = sd(Thickness, na.rm=TRUE),
CV = (sd / mean)
)
knitr::kable(CVa)| Treatment | Species | LifeStage | mean | sd | CV |
|---|---|---|---|---|---|
| Ambient | Montipora capitata | Adult | 1607.8606 | 328.4657 | 0.2042874 |
| Ambient | Montipora capitata | Juvenile | 1452.3133 | 255.5347 | 0.1759501 |
| Ambient | Pocillopora acuta | Adult | 763.4843 | 120.2092 | 0.1574482 |
| Ambient | Pocillopora acuta | Juvenile | 893.8610 | 129.7337 | 0.1451386 |
| High | Montipora capitata | Adult | 1474.2514 | 306.4991 | 0.2079015 |
| High | Montipora capitata | Juvenile | 1328.9706 | 244.6491 | 0.1840892 |
| High | Pocillopora acuta | Adult | 645.2095 | 177.2757 | 0.2747567 |
| High | Pocillopora acuta | Juvenile | 785.4726 | 160.2471 | 0.2040136 |
2. LSCM Symbiont Fluorescence: Initial and endpoints
Load dataset and calculate mean symbiont fluorescence values for each individual coral at start and ending timepoints collected by LSCM.
#load dataset
bleach_red<-read.csv("Data/Bleaching/Bleaching_Confocal_Fluorescence.csv", header=TRUE, sep=",", na.strings="NA")
#calculate mean for each sample
bleach2 <- ddply(bleach_red, c("Timepoint", "Sample"), summarise,
Intensity = mean(Symbiont.Intensity, na.rm=TRUE)
)
#match identifying values into this data frame
bleach2$Species<-bleach_red$Species[match(bleach2$Sample, bleach_red$Sample)]
bleach2$LifeStage<-bleach_red$LifeStage[match(bleach2$Sample, bleach_red$Sample)]
bleach2$Treatment<-bleach_red$Treatment[match(bleach2$Sample, bleach_red$Sample)]
bleach2$Unique<-paste(bleach2$Sample, bleach2$Timepoint)
bleach2$Colony<-bleach_red$Colony[match(bleach2$Sample, bleach_red$Sample)]
#conduct calibrations to convert intensity values to Red %RI values
#Use conditional calculations by subsetting each species, then matching the values back into the master data sheets
mcap_cal<-subset(bleach2, Species=="Montipora capitata", select=c(Sample, Treatment, Intensity, Unique))
mcap_cal$Red<-(mcap_cal$Intensity+268.36)/109.01
pacu_cal<-subset(bleach2, Species=="Pocillopora acuta", select=c(Sample, Treatment, Intensity, Unique))
pacu_cal$Red<-(pacu_cal$Intensity+464.76)/193.16
#merge all Red values into master data frame
calibrations<-bind_rows(mcap_cal, pacu_cal)
#combine into data sheet
bleach2$Red<-calibrations$Red[match(bleach2$Unique, calibrations$Unique)]
#set the intial timepoint as the reference level
bleach_red$Timepoint <- relevel(bleach_red$Timepoint, ref = "Initial")
#set the intial timepoint as the reference level
bleach2$Timepoint <- relevel(bleach2$Timepoint, ref = "Initial")
#summarize by species at each timepoint and treatment
bleachRed <- ddply(bleach2, c("Timepoint", "Treatment", "Species", "LifeStage"), summarise,
N = length(Red[!is.na(Red)]),
mean = mean(Red, na.rm=TRUE),
sd = sd(Red, na.rm=TRUE),
se = sd / sqrt(N)
)
knitr::kable(bleachRed)
#subset by species and add a group name for analysis of homogeneity of variance
bleach2$Group<-paste(bleach2$LifeStage, bleach2$Treatment)
bleach2$Colony<-bleach1$Colony[match(bleach2$Sample, bleach1$Sample)]
bleachRed_mcap<-subset(bleach2, Species=="Montipora capitata", c(Treatment, Timepoint, Species, LifeStage, Red, Sample, Group, Colony))
bleachRed_pacu<-subset(bleach2, Species=="Pocillopora acuta", c(Treatment, Timepoint, Species, LifeStage, Red, Sample, Group, Colony))Display a figure of mean symbiont fluorescence over the course of bleaching with a separate plot for each species.
FigBleach2 <- ggplot(data=bleachRed, aes(x=Timepoint, y=mean, group=interaction(Treatment, LifeStage), colour=Treatment, linetype=LifeStage)) +
facet_wrap(~Species, scale="free")+
geom_line(position=position_dodge(0.09), size=1.3) +
scale_linetype_manual(name = "Life Stage", values=c("solid", "dotted"))+
scale_color_manual(name = "Temperature", values=c("blue", "red")) +
geom_point(size=3, position=position_dodge(0.09)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0, linetype="solid", position=position_dodge(0.09), size=1.3) +
xlab("Timepoint") +
ylab(expression(bold(paste("Symbiodiniaceae Fluorescence (% RI)"))))+
ylim(0,75) +
theme_bw()+
theme_classic() +
theme(panel.background = element_blank()) +
theme(strip.text.x=element_text(face="italic")) +
theme(text = element_text(size=20),
panel.border = element_blank(),
plot.background =element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_text(size=18, colour="black"),
legend.title=element_text(face="bold"),
legend.position="right",
plot.title=element_text(),
legend.text = element_text(size = 20),
axis.title.y=element_text(margin=margin(), face="bold"),
axis.title.x=element_text(margin=margin(), face="bold"))
FigBleach2
ggsave(filename="Figures/Fig2b.pdf", plot=FigBleach2, dpi=300, width=10, height=7, units="in")Analyze the differences in symbiont fluorescence over time between temperature and life stage with a separate model for each species.
Montipora capitata:
modelB2a<-lmer(Red~LifeStage * Treatment * Timepoint + (1|Sample) + (1|Colony), data=bleachRed_mcap)
summary(modelB2a)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Red ~ LifeStage * Treatment * Timepoint + (1 | Sample) + (1 |
## Colony)
## Data: bleachRed_mcap
##
## REML criterion at convergence: 385.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94242 -0.52737 -0.00445 0.64938 2.29855
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample (Intercept) 5.873e+00 2.423e+00
## Colony (Intercept) 3.919e-13 6.260e-07
## Residual 7.682e+01 8.765e+00
## Number of obs: 59, groups: Sample, 30; Colony, 16
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 48.361 3.712 50.745
## LifeStageJuvenile 1.540 4.911 50.745
## TreatmentHigh -7.672 5.059 50.745
## TimepointFinal -8.488 5.060 25.031
## LifeStageJuvenile:TreatmentHigh 1.747 6.717 50.745
## LifeStageJuvenile:TimepointFinal -5.796 6.803 25.793
## TreatmentHigh:TimepointFinal -5.592 6.896 25.031
## LifeStageJuvenile:TreatmentHigh:TimepointFinal 1.909 9.236 25.442
## t value Pr(>|t|)
## (Intercept) 13.027 <2e-16 ***
## LifeStageJuvenile 0.314 0.755
## TreatmentHigh -1.516 0.136
## TimepointFinal -1.677 0.106
## LifeStageJuvenile:TreatmentHigh 0.260 0.796
## LifeStageJuvenile:TimepointFinal -0.852 0.402
## TreatmentHigh:TimepointFinal -0.811 0.425
## LifeStageJuvenile:TreatmentHigh:TimepointFinal 0.207 0.838
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) LfStgJ TrtmnH TmpntF LfSJ:TH LSJ:TF TrH:TF
## LifeStgJvnl -0.756
## TreatmntHgh -0.734 0.555
## TimepontFnl -0.682 0.515 0.500
## LfStgJvn:TH 0.553 -0.731 -0.753 -0.377
## LfStgJvn:TF 0.507 -0.671 -0.372 -0.744 0.490
## TrtmntHg:TF 0.500 -0.378 -0.682 -0.734 0.513 0.546
## LfStJ:TH:TF -0.373 0.494 0.509 0.548 -0.676 -0.737 -0.747
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(modelB2a, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## LifeStage 0.26 0.26 1 25.448 0.0034 0.953662
## Treatment 1011.98 1011.98 1 25.621 13.1735 0.001239
## Timepoint 2959.20 2959.20 1 25.462 38.5217 1.586e-06
## LifeStage:Treatment 22.60 22.60 1 25.513 0.2942 0.592271
## LifeStage:Timepoint 82.26 82.26 1 25.378 1.0708 0.310528
## Treatment:Timepoint 74.81 74.81 1 25.551 0.9739 0.332961
## LifeStage:Treatment:Timepoint 3.28 3.28 1 25.442 0.0427 0.837893
##
## LifeStage
## Treatment **
## Timepoint ***
## LifeStage:Treatment
## LifeStage:Timepoint
## Treatment:Timepoint
## LifeStage:Treatment:Timepoint
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm = emmeans(modelB2a, ~ LifeStage, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## LifeStage emmean SE df lower.CL upper.CL .group
## Juvenile 38.9 1.66 13.5 35.3 42.4 A
## Adult 38.9 1.86 13.1 34.9 42.9 A
##
## Results are averaged over the levels of: Treatment, Timepoint
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Adult - Juvenile 0.00761 2.49 13.3 0.003 0.9976
##
## Results are averaged over the levels of: Treatment, Timepoint
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB2a, ~ Treatment, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## Treatment emmean SE df lower.CL upper.CL .group
## High 34.3 1.68 25.3 30.9 37.8 A
## Ambient 43.4 1.85 26.1 39.6 47.2 B
##
## Results are averaged over the levels of: LifeStage, Timepoint
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Ambient - High 9.12 2.49 13.3 3.656 0.0028
##
## Results are averaged over the levels of: LifeStage, Timepoint
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB2a, ~ Timepoint, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) #letter display## Timepoint emmean SE df lower.CL upper.CL .group
## Final 32.0 1.72 33.6 28.5 35.5 A
## Initial 45.7 1.68 32.8 42.3 49.2 B
##
## Results are averaged over the levels of: LifeStage, Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Initial - Final 13.7 2.31 25.6 5.928 <.0001
##
## Results are averaged over the levels of: LifeStage, Treatment
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB2a, ~ Treatment * LifeStage * Timepoint, adjust="tukey") #can choose other adjustment methods
cld(emm, Letters=c(LETTERS)) ## Treatment LifeStage Timepoint emmean SE df lower.CL upper.CL .group
## High Juvenile Final 26.0 3.03 50.8 19.9 32.1 A
## High Adult Final 26.6 3.44 50.8 19.7 33.5 AB
## Ambient Juvenile Final 35.6 3.49 50.9 28.6 42.6 ABC
## Ambient Adult Final 39.9 3.74 50.8 32.4 47.4 ABC
## High Adult Initial 40.7 3.44 50.8 33.8 47.6 BC
## High Juvenile Initial 44.0 3.03 50.8 37.9 50.1 C
## Ambient Adult Initial 48.4 3.74 50.8 40.9 55.9 C
## Ambient Juvenile Initial 49.9 3.23 50.8 43.4 56.4 C
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 8 estimates
## significance level used: alpha = 0.05Assumptions of analysis are met.
qqPlot(residuals(modelB2a))
## 41 45
## 19 21bartlett.test(residuals(modelB2a)~bleachRed_mcap$Group) #not violated##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelB2a) by bleachRed_mcap$Group
## Bartlett's K-squared = 1.985, df = 3, p-value = 0.5755There is a significant negative effect of high temperature treatment and a decrease in values across groups at the final timepoint in Montipora capitata.
Pocillopora acuta:
modelB2b<-lmer(Red~LifeStage * Treatment * Timepoint + (1|Sample) + (1|Colony), data=bleachRed_pacu)
summary(modelB2b)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Red ~ LifeStage * Treatment * Timepoint + (1 | Sample) + (1 |
## Colony)
## Data: bleachRed_pacu
##
## REML criterion at convergence: 381.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.55939 -0.58201 0.00185 0.48205 1.87322
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample (Intercept) 35.3330 5.9442
## Colony (Intercept) 0.8215 0.9063
## Residual 90.0305 9.4884
## Number of obs: 56, groups: Sample, 28; Colony, 14
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 42.731 4.246 44.356
## LifeStageJuvenile 8.259 6.004 44.356
## TreatmentHigh 1.762 5.985 25.148
## TimepointFinal -5.420 5.072 23.999
## LifeStageJuvenile:TreatmentHigh -0.508 8.464 25.148
## LifeStageJuvenile:TimepointFinal -3.978 7.173 23.999
## TreatmentHigh:TimepointFinal -23.779 7.173 23.999
## LifeStageJuvenile:TreatmentHigh:TimepointFinal 5.998 10.144 23.999
## t value Pr(>|t|)
## (Intercept) 10.065 5e-13 ***
## LifeStageJuvenile 1.376 0.1759
## TreatmentHigh 0.294 0.7709
## TimepointFinal -1.069 0.2959
## LifeStageJuvenile:TreatmentHigh -0.060 0.9526
## LifeStageJuvenile:TimepointFinal -0.555 0.5843
## TreatmentHigh:TimepointFinal -3.315 0.0029 **
## LifeStageJuvenile:TreatmentHigh:TimepointFinal 0.591 0.5598
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) LfStgJ TrtmnH TmpntF LfSJ:TH LSJ:TF TrH:TF
## LifeStgJvnl -0.707
## TreatmntHgh -0.705 0.498
## TimepontFnl -0.597 0.422 0.424
## LfStgJvn:TH 0.498 -0.705 -0.707 -0.300
## LfStgJvn:TF 0.422 -0.597 -0.300 -0.707 0.424
## TrtmntHg:TF 0.422 -0.299 -0.599 -0.707 0.424 0.500
## LfStJ:TH:TF -0.299 0.422 0.424 0.500 -0.599 -0.707 -0.707anova(modelB2b, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## LifeStage 434.2 434.2 1 11.949 4.8229 0.048564 *
## Treatment 618.8 618.8 1 11.943 6.8734 0.022389 *
## Timepoint 4435.2 4435.2 1 23.999 49.2635 2.949e-07 ***
## LifeStage:Treatment 12.2 12.2 1 11.943 0.1352 0.719577
## LifeStage:Timepoint 3.4 3.4 1 23.999 0.0373 0.848568
## Treatment:Timepoint 1511.4 1511.4 1 23.999 16.7872 0.000412 ***
## LifeStage:Treatment:Timepoint 31.5 31.5 1 23.999 0.3497 0.559837
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm = emmeans(modelB2b, ~ LifeStage, adjust="tukey")
cld(emm, Letters=c(LETTERS)) #letter display## LifeStage emmean SE df lower.CL upper.CL .group
## Adult 35.0 2.42 12 29.7 40.2 A
## Juvenile 42.5 2.42 12 37.2 47.7 B
##
## Results are averaged over the levels of: Treatment, Timepoint
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Adult - Juvenile -7.52 3.42 12 -2.196 0.0485
##
## Results are averaged over the levels of: Treatment, Timepoint
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB2b, ~ Treatment, adjust="tukey")
cld(emm, Letters=c(LETTERS)) #letter display## Treatment emmean SE df lower.CL upper.CL .group
## High 34.3 2.41 24 29.3 39.2 A
## Ambient 43.2 2.41 24 38.2 48.1 B
##
## Results are averaged over the levels of: LifeStage, Timepoint
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Ambient - High 8.88 3.39 12 2.622 0.0223
##
## Results are averaged over the levels of: LifeStage, Timepoint
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB2b, ~ Timepoint, adjust="tukey")
cld(emm, Letters=c(LETTERS)) #letter display## Timepoint emmean SE df lower.CL upper.CL .group
## Final 29.8 2.13 25 25.4 34.2 A
## Initial 47.6 2.13 25 43.2 52.0 B
##
## Results are averaged over the levels of: LifeStage, Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05pairs(emm)## contrast estimate SE df t.ratio p.value
## Initial - Final 17.8 2.54 24 7.019 <.0001
##
## Results are averaged over the levels of: LifeStage, Treatment
## Degrees-of-freedom method: kenward-rogeremm = emmeans(modelB2b, ~ Timepoint * Treatment * LifeStage, adjust="tukey")
cld(emm, Letters=c(LETTERS)) #letter display## Timepoint Treatment LifeStage emmean SE df lower.CL upper.CL .group
## Final High Adult 15.3 4.25 44.4 6.74 23.8 A
## Final High Juvenile 25.1 4.25 44.4 16.51 33.6 AB
## Final Ambient Adult 37.3 4.25 44.4 28.76 45.9 BC
## Final Ambient Juvenile 41.6 4.25 44.4 33.04 50.1 BC
## Initial Ambient Adult 42.7 4.25 44.4 34.18 51.3 BC
## Initial High Adult 44.5 4.25 44.4 35.94 53.0 C
## Initial Ambient Juvenile 51.0 4.25 44.4 42.44 59.5 C
## Initial High Juvenile 52.2 4.25 44.4 43.69 60.8 C
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 8 estimates
## significance level used: alpha = 0.05Assumptions of analysis are met.
qqPlot(residuals(modelB2b))
## 63 89
## 34 42bartlett.test(residuals(modelB2b)~bleachRed_pacu$Group) #not violated##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelB2b) by bleachRed_pacu$Group
## Bartlett's K-squared = 3.7958, df = 3, p-value = 0.2844There was a significant effect of lifestage, treatment, timepoint, and the interaction of treatment and timepoint on Pocillopora acuta symbiont fluorescence. This indicates there was a decrease in symbiont fluorescence by the final timepoint at high temperature and an increase in fluorescence in juvenile corals as compared to adults.
Coefficient of variation (CV):
sd(bleach2$Red)/mean(bleach2$Red)## [1] 0.3639065CVb <- ddply(bleach2, c("Treatment", "Species", "LifeStage"), summarise,
mean = mean(Red, na.rm=TRUE),
sd = sd(Red, na.rm=TRUE),
CV = (sd / mean)
)
knitr::kable(CVb)| Treatment | Species | LifeStage | mean | sd | CV |
|---|---|---|---|---|---|
| Ambient | Montipora capitata | Adult | 44.11683 | 10.516427 | 0.2383768 |
| Ambient | Montipora capitata | Juvenile | 43.18784 | 10.147932 | 0.2349720 |
| Ambient | Pocillopora acuta | Adult | 40.02160 | 8.856768 | 0.2212997 |
| Ambient | Pocillopora acuta | Juvenile | 46.29210 | 12.926205 | 0.2792313 |
| High | Montipora capitata | Adult | 33.64837 | 12.154375 | 0.3612174 |
| High | Montipora capitata | Juvenile | 34.99137 | 12.773443 | 0.3650456 |
| High | Pocillopora acuta | Adult | 29.89379 | 18.328324 | 0.6131148 |
| High | Pocillopora acuta | Juvenile | 38.65536 | 18.562573 | 0.4802069 |
3. Comparative metrics: Chlorophyll, Thickness by Dissecting Microscope, and Pigmentation/Color at the endpoint of bleaching exposure
Calculate mean chlorophyll a content, color, and tissue thickness (by dissecting scope method) values in Ambient and High treatments and analyze these values with ANOVA tests and linear mixed effect model analysis. Repeated measures will be included in models using the paired colony ID. Montipora capitata juvenile colonies did not have a pair and are therefore randomly paired with another individual.
A) Chlorophyll Content
First, assemble data frame for chlorophyll a content averaged for each coral at High and Ambient temperature treatments.
#chlorophyll
compareChl<-read.csv("Data/Bleaching/Bleaching_Chlorophyll_Content.csv", header=TRUE, sep=",", na.strings="NA")
compareChl$Colony<-as.factor(compareChl$Colony)
compareChl$Treatment<-bleach2$Treatment[match(compareChl$Sample, bleach2$Sample)]
compareChl$Species<-bleach2$Species[match(compareChl$Sample, bleach2$Sample)]
compareChl$LifeStage<-bleach2$LifeStage[match(compareChl$Sample, bleach2$Sample)]
compareChl<-compareChl[c("Species", "LifeStage", "Colony", "Chla.SA", "Treatment", "Sample")]
compareChl_mcap <- compareChl[which(compareChl$Species == "Montipora capitata"),]
compareChl_mcap <- compareChl_mcap[which(compareChl_mcap$LifeStage == "Adult"),]
compareChl_pacu <- compareChl[which(compareChl$Species == "Pocillopora acuta"),]Summarize and display a chart of mean chlorophyll content.
compare3<- ddply(compareChl, c("Species", "LifeStage", "Treatment"), summarise,
N = length(Chla.SA[!is.na(Chla.SA)]),
mean = mean(Chla.SA, na.rm=TRUE),
sd = sd(Chla.SA, na.rm=TRUE),
se = sd / sqrt(N)
)
compare3## Species LifeStage Treatment N mean sd se
## 1 Montipora capitata Adult Ambient 3 0.15600685 0.10204437 0.058915342
## 2 Montipora capitata Adult High 4 0.08002048 0.02627553 0.013137764
## 3 Montipora capitata Juvenile Ambient 0 NaN NA NA
## 4 Montipora capitata Juvenile High 0 NaN NA NA
## 5 Pocillopora acuta Adult Ambient 4 0.11600707 0.07824723 0.039123613
## 6 Pocillopora acuta Adult High 4 0.02736657 0.01003481 0.005017407
## 7 Pocillopora acuta Juvenile Ambient 4 0.12500527 0.06831860 0.034159298
## 8 Pocillopora acuta Juvenile High 4 0.04144988 0.01100008 0.005500039CompareBleach3 <- ggplot(data=compare3, aes(x=Treatment, y=mean, group=LifeStage, shape=LifeStage, color=Treatment)) +
scale_shape_manual(name = "Life Stage", values = c(19,21))+
facet_wrap(~Species)+
scale_color_manual(values=c("blue", "red"))+
geom_point(aes(color=Treatment), size=3, position=position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black", width=0, linetype="solid", position=position_dodge(0.3), size=0.8) +
xlab("Treatment") +
ylab(expression(paste("Chlorophyll content (", mu, "g cm"^-2, ")"))) +
ylim(-0,0.3) +
theme_bw()+
theme_classic() +
theme(strip.text.x=element_text(face="italic")) +
theme(panel.background = element_blank()) +
theme(text = element_text(size=20),
panel.border = element_blank(),
plot.background =element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.y = element_text(size=18, colour="black"),
axis.text.x = element_text(size=18, colour="black", face="plain"),
legend.title=element_text(),
legend.position="none",
plot.title=element_text(),
legend.text = element_text(size = 20),
axis.title.y=element_text(margin=margin(t = 0, r = 5, b = 0, l = 0), face="bold"),
axis.title.x=element_text(margin=margin(t = 5, r = 0, b = 0, l = 0), face="bold"))
CompareBleach3
ggsave(filename="Figures/FigS1b.pdf", plot=CompareBleach3, dpi=300, width=6.5, height=8, units="in")Not that chlorophyll extraction from juvenile M. capitata was not successful and is not included in the dataset here.
Analyze effect of treatment and lifestage on chlorophyll content with a separate model for each species.
Montipora capitata:
modelB3a<-lmer(Chla.SA ~ Treatment + (1|Colony), data=compareChl_mcap)
summary(modelB3a)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Chla.SA ~ Treatment + (1 | Colony)
## Data: compareChl_mcap
##
## REML criterion at convergence: -10.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.18437 -0.35102 -0.00477 0.20307 1.48504
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.001035 0.03217
## Residual 0.003497 0.05914
## Number of obs: 7, groups: Colony, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.15598 0.03861 4.97765 4.039 0.010 *
## TreatmentHigh -0.07596 0.04590 2.93370 -1.655 0.199
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TreatmntHgh -0.695anova(modelB3a, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 0.0095785 0.0095785 1 2.9337 2.7388 0.1986Assumptions of analysis are not met.
qqPlot(residuals(modelB3a))
## 21 25
## 2 6bartlett.test(residuals(modelB3a)~compareChl_mcap$Treatment) #not violated##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelB3a) by compareChl_mcap$Treatment
## Bartlett's K-squared = 3.563, df = 1, p-value = 0.05908kruskal.test(compareChl_mcap$Chla.SA~compareChl_mcap$Treatment) #not significant##
## Kruskal-Wallis rank sum test
##
## data: compareChl_mcap$Chla.SA by compareChl_mcap$Treatment
## Kruskal-Wallis chi-squared = 1.125, df = 1, p-value = 0.2888#conduct non parametric post hoc tets of kruskal test
chl_ph = dunnTest(Chla.SA ~ Treatment,
data=compareChl_mcap,
method="bonferroni") # Can adjust p-values;
chl_ph## Comparison Z P.unadj P.adj
## 1 Ambient - High 1.06066 0.2888444 0.2888444There is no effect of treatment in M. capitata.
Pocillopora acuta:
modelB3b<-lmer(Chla.SA~LifeStage * Treatment + (1|Colony), data=compareChl_pacu)
summary(modelB3b)
anova(modelB3b, type=2)Assumptions of analysis are met.
qqPlot(residuals(modelB3b))
## 29 15
## 15 7compareChl_pacu$Group<-paste(compareChl_pacu$Treatment, compareChl_pacu$LifeStage)
bartlett.test(residuals(modelB3b)~compareChl_pacu$Group) #violated due to large variability in ambient treatment as compared to high##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelB3b) by compareChl_pacu$Group
## Bartlett's K-squared = 14.495, df = 3, p-value = 0.002304As assumption of homogeneity of variance is not met, we proceed with analyzing effect of treatment and lifestage using a Kruskal Wallis Test.
kruskal.test(compareChl_pacu$Chla.SA~compareChl_pacu$Treatment) #significant##
## Kruskal-Wallis rank sum test
##
## data: compareChl_pacu$Chla.SA by compareChl_pacu$Treatment
## Kruskal-Wallis chi-squared = 8.0404, df = 1, p-value = 0.004574kruskal.test(compareChl_pacu$Chla.SA~compareChl_pacu$LifeStage) #not significant##
## Kruskal-Wallis rank sum test
##
## data: compareChl_pacu$Chla.SA by compareChl_pacu$LifeStage
## Kruskal-Wallis chi-squared = 1.3346, df = 1, p-value = 0.248chl_ph2 = dunnTest(Chla.SA ~ Treatment,
data=compareChl_pacu,
method="bonferroni")
chl_ph2## Comparison Z P.unadj P.adj
## 1 Ambient - High 2.835567 0.004574439 0.004574439chl_ph3 = dunnTest(Chla.SA ~ LifeStage,
data=compareChl_pacu,
method="bonferroni")
chl_ph3## Comparison Z P.unadj P.adj
## 1 Adult - Juvenile -1.155231 0.2479958 0.2479958There is a significant effect of treatment on chlorophyll content. There is no effect of lifestage.
Coefficient of variation (CV):
sd(compareChl$Chla.SA, na.rm=TRUE)/mean(compareChl$Chla.SA, na.rm=TRUE)## [1] 0.7726447CVc <- ddply(compareChl, c("Treatment", "Species", "LifeStage"), summarise,
mean = mean(Chla.SA, na.rm=TRUE),
sd = sd(Chla.SA, na.rm=TRUE),
CV = (sd / mean)
)
knitr::kable(CVc)| Treatment | Species | LifeStage | mean | sd | CV |
|---|---|---|---|---|---|
| Ambient | Montipora capitata | Adult | 0.1560069 | 0.1020444 | 0.6541018 |
| Ambient | Montipora capitata | Juvenile | NaN | NA | NA |
| Ambient | Pocillopora acuta | Adult | 0.1160071 | 0.0782472 | 0.6745039 |
| Ambient | Pocillopora acuta | Juvenile | 0.1250053 | 0.0683186 | 0.5465257 |
| High | Montipora capitata | Adult | 0.0800205 | 0.0262755 | 0.3283600 |
| High | Montipora capitata | Juvenile | NaN | NA | NA |
| High | Pocillopora acuta | Adult | 0.0273666 | 0.0100348 | 0.3666814 |
| High | Pocillopora acuta | Juvenile | 0.0414499 | 0.0110001 | 0.2653826 |
B) Dissecting Scope Tissue Thickness
Assemble data frame of dissecting scope tissue thickness for each coral at Ambient and High treatments.
#chlorophyll
compareThick<-read.csv("Data/Bleaching/Bleaching_Scope_Thickness.csv", header=TRUE, sep=",", na.strings="NA")
compareThick$Colony<-as.factor(compareThick$Colony)
compareThick<-compareThick[c("Species", "LifeStage", "Colony", "Thickness", "Treatment", "Sample")]
compareThick_mcap <- compareThick[which(compareThick$Species == "Montipora capitata"),]
compareThick_pacu <- compareThick[which(compareThick$Species == "Pocillopora acuta"),]Summarize and display a figure of mean tissue thickness by dissecting microscope.
compare4<- ddply(compareThick, c("Species", "LifeStage", "Treatment"), summarise,
N = length(Thickness[!is.na(Thickness)]),
mean = mean(Thickness, na.rm=TRUE),
sd = sd(Thickness, na.rm=TRUE),
se = sd / sqrt(N)
)
compare4## Species LifeStage Treatment N mean sd se
## 1 Montipora capitata Adult Ambient 3 1695.0000 111.25674 64.23411
## 2 Montipora capitata Adult High 3 1916.8333 81.15738 46.85624
## 3 Montipora capitata Juvenile Ambient 3 1584.3333 252.42775 145.73923
## 4 Montipora capitata Juvenile High 5 1303.4500 218.63457 97.77635
## 5 Pocillopora acuta Adult Ambient 3 602.3333 77.75214 44.89022
## 6 Pocillopora acuta Adult High 3 571.9167 149.79576 86.48462
## 7 Pocillopora acuta Juvenile Ambient 3 617.1667 149.06004 86.05985
## 8 Pocillopora acuta Juvenile High 3 561.7500 59.71861 34.47856CompareBleach4 <- ggplot(data=compare4, aes(x=Treatment, y=mean, group=LifeStage, shape=LifeStage, color=Treatment)) +
scale_shape_manual(name = "Life Stage", values = c(19,21))+
facet_wrap(~Species)+
scale_color_manual(values=c("blue", "red"))+
geom_point(aes(color=Treatment), size=3, position=position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black", width=0, linetype="solid", position=position_dodge(0.3), size=0.8) +
xlab("Treatment") +
ylab(expression(paste("Tissue Thickness (", mu, "m)"))) +
ylim(0,2800) +
theme_bw()+
theme_classic() +
theme(strip.text.x=element_text(face="italic")) +
theme(panel.background = element_blank()) +
theme(text = element_text(size=20),
panel.border = element_blank(),
plot.background =element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.y = element_text(size=18, colour="black"),
axis.text.x = element_text(size=18, colour="black", face="plain"),
legend.title=element_text(),
legend.position="none",
plot.title=element_text(),
legend.text = element_text(size = 20),
axis.title.y=element_text(margin=margin(t = 0, r = 5, b = 0, l = 0), face="bold"),
axis.title.x=element_text(margin=margin(t = 5, r = 0, b = 0, l = 0), face="bold"))
CompareBleach4
ggsave(filename="Figures/FigS1a.pdf", plot=CompareBleach4, dpi=300, width=6.5, height=8, units="in")Analyze effect of treatment and lifestage on tissue thickness with a separate model for each species.
Montipora capitata:
modelB4a<-lmer(Thickness ~ Treatment * LifeStage + (1|Colony), data=compareThick_mcap)
summary(modelB4a)
anova(modelB4a, type=2)Assumptions of analysis are not met.
qqPlot(residuals(modelB4a))
## 9 8
## 3 2compareThick_mcap$Group<-paste(compareThick_mcap$Treatment, compareThick_mcap$LifeStage)
compareThick_mcap$Group<-as.factor(compareThick_mcap$Group)
bartlett.test(residuals(modelB4a)~compareThick_mcap$Group) #not violated##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelB4a) by compareThick_mcap$Group
## Bartlett's K-squared = 5.5791, df = 3, p-value = 0.134kruskal.test(compareThick_mcap$Thickness ~ compareThick_mcap$Treatment)##
## Kruskal-Wallis rank sum test
##
## data: compareThick_mcap$Thickness by compareThick_mcap$Treatment
## Kruskal-Wallis chi-squared = 0.26667, df = 1, p-value = 0.6056kruskal.test(compareThick_mcap$Thickness ~ compareThick_mcap$LifeStage)##
## Kruskal-Wallis rank sum test
##
## data: compareThick_mcap$Thickness by compareThick_mcap$LifeStage
## Kruskal-Wallis chi-squared = 6.6667, df = 1, p-value = 0.009823thick_ph2 = dunnTest(Thickness ~ Treatment,
data=compareThick_mcap,
method="bonferroni")
thick_ph2## Comparison Z P.unadj P.adj
## 1 Ambient - High 0.5163978 0.6055766 0.6055766thick_ph3 = dunnTest(Thickness ~ LifeStage,
data=compareThick_mcap,
method="bonferroni")
thick_ph3## Comparison Z P.unadj P.adj
## 1 Adult - Juvenile 2.581989 0.009823275 0.009823275There is an effect of lifestage in M. capitata.
Pocillopora acuta:
kruskal.test(compareThick_pacu$Thickness ~ compareThick_pacu$Treatment)##
## Kruskal-Wallis rank sum test
##
## data: compareThick_pacu$Thickness by compareThick_pacu$Treatment
## Kruskal-Wallis chi-squared = 0.23077, df = 1, p-value = 0.631kruskal.test(compareThick_pacu$Thickness ~ compareThick_pacu$LifeStage)##
## Kruskal-Wallis rank sum test
##
## data: compareThick_pacu$Thickness by compareThick_pacu$LifeStage
## Kruskal-Wallis chi-squared = 0, df = 1, p-value = 1thick_ph2 = dunnTest(Thickness ~ Treatment,
data=compareThick_pacu,
method="bonferroni")
thick_ph2## Comparison Z P.unadj P.adj
## 1 Ambient - High 0.4803845 0.630954 0.630954thick_ph3 = dunnTest(Thickness ~ LifeStage,
data=compareThick_pacu,
method="bonferroni")
thick_ph3## Comparison Z P.unadj P.adj
## 1 Adult - Juvenile 0 1 1There was no effect of treatment or lifestage on tissue thickness in P. acuta.
Coefficient of variation (CV):
sd(compareThick$Thickness, na.rm=TRUE)/mean(compareThick$Thickness, na.rm=TRUE)## [1] 0.4907275CVd <- ddply(compareThick, c("Treatment", "Species", "LifeStage"), summarise,
mean = mean(Thickness, na.rm=TRUE),
sd = sd(Thickness, na.rm=TRUE),
CV = (sd / mean)
)
knitr::kable(CVd)| Treatment | Species | LifeStage | mean | sd | CV |
|---|---|---|---|---|---|
| Ambient | Montipora capitata | Adult | 1695.0000 | 111.25674 | 0.0656382 |
| Ambient | Montipora capitata | Juvenile | 1584.3333 | 252.42775 | 0.1593274 |
| Ambient | Pocillopora acuta | Adult | 602.3333 | 77.75214 | 0.1290849 |
| Ambient | Pocillopora acuta | Juvenile | 617.1667 | 149.06004 | 0.2415232 |
| High | Montipora capitata | Adult | 1916.8333 | 81.15738 | 0.0423393 |
| High | Montipora capitata | Juvenile | 1303.4500 | 218.63457 | 0.1677353 |
| High | Pocillopora acuta | Adult | 571.9167 | 149.79576 | 0.2619189 |
| High | Pocillopora acuta | Juvenile | 561.7500 | 59.71861 | 0.1063082 |
C) Coral Pigmentation / Color Scores
Assemble data frame of coral pigmentation measured as a ratio of coral color : white in Fiji.
compareColor<-read.csv("Data/Bleaching/Bleaching_Coral_Color.csv", header=TRUE, sep=",", na.strings="NA")
compareColor$Colony<-as.factor(compareColor$Colony)
compareColor<-compareColor[c("Species", "LifeStage", "Colony", "Coral.Percent.White", "Treatment", "Sample")]
compareColor_mcap <- compareColor[which(compareColor$Species == "Montipora capitata"),]
compareColor_pacu <- compareColor[which(compareColor$Species == "Pocillopora acuta"),]Summarize and display a chart of mean pigmentation (coral:white).
compare5<- ddply(compareColor, c("Species", "LifeStage", "Treatment"), summarise,
N = length(Coral.Percent.White[!is.na(Coral.Percent.White)]),
mean = mean(Coral.Percent.White, na.rm=TRUE),
sd = sd(Coral.Percent.White, na.rm=TRUE),
se = sd / sqrt(N)
)
compare5## Species LifeStage Treatment N mean sd se
## 1 Montipora capitata Adult Ambient 6 0.2823965 0.05953676 0.02430578
## 2 Montipora capitata Adult High 7 0.4716676 0.13651399 0.05159744
## 3 Montipora capitata Juvenile Ambient 7 0.2442736 0.06182641 0.02336819
## 4 Montipora capitata Juvenile High 9 0.3280025 0.07647571 0.02549190
## 5 Pocillopora acuta Adult Ambient 7 0.4279750 0.08226357 0.03109271
## 6 Pocillopora acuta Adult High 7 0.5998219 0.08748308 0.03306550
## 7 Pocillopora acuta Juvenile Ambient 7 0.4596533 0.06023771 0.02276772
## 8 Pocillopora acuta Juvenile High 7 0.5009551 0.04131316 0.01561490CompareBleach5 <- ggplot(data=compare5, aes(x=Treatment, y=mean, group=LifeStage, shape=LifeStage, color=Treatment)) +
scale_shape_manual(name = "Life Stage", values = c(19,21))+
facet_wrap(~Species)+
scale_color_manual(values=c("blue", "red"))+
geom_point(aes(color=Treatment), size=3, position=position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black", width=0, linetype="solid", position=position_dodge(0.3), size=0.8) +
xlab("Treatment") +
ylab(expression(paste("Bleaching Score (Coral : White)"))) +
ylim(0,0.75) +
theme_bw()+
theme_classic() +
theme(strip.text.x=element_text(face="italic")) +
theme(panel.background = element_blank()) +
theme(text = element_text(size=20),
panel.border = element_blank(),
plot.background =element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.y = element_text(size=18, colour="black"),
axis.text.x = element_text(size=18, colour="black", face="plain"),
legend.title=element_text(),
legend.position="right",
plot.title=element_text(),
legend.text = element_text(size = 20),
axis.title.y=element_text(margin=margin(t = 0, r = 5, b = 0, l = 0), face="bold"),
axis.title.x=element_text(margin=margin(t = 5, r = 0, b = 0, l = 0), face="bold"))
CompareBleach5
ggsave(filename="Figures/FigS1c.pdf", plot=CompareBleach5, dpi=300, width=8, height=8, units="in")Analyze effect of treatment and lifestage on color with a separate model for each species.
Montipora capitata:
modelB5a<-lmer(Coral.Percent.White ~ Treatment * LifeStage + (1|Colony), data=compareColor_mcap)
summary(modelB5a)
anova(modelB5a, type=2)
qqPlot(residuals(modelB5a)) #violated 
Assumptions of analysis are not met in a linear model (normality) and are not improved with transformation.
compareColor_mcap$tdata<-sqrt(compareColor_mcap$Coral.Percent.White)
modelB5aT<-lmer(tdata ~ Treatment * LifeStage + (1|Colony), data=compareColor_mcap)
qqPlot(residuals(modelB5aT)) #violated
compareColor_mcap$Group<-paste(compareColor_mcap$Treatment, compareColor_mcap$LifeStage)
compareColor_mcap$Group<-as.factor(compareColor_mcap$Group)
bartlett.test(residuals(modelB5aT)~compareColor_mcap$Group) #not violatedConduct a non-parametric analysis.
kruskal.test(compareColor_mcap$Coral.Percent.White ~ compareColor_mcap$Treatment)##
## Kruskal-Wallis rank sum test
##
## data: compareColor_mcap$Coral.Percent.White by compareColor_mcap$Treatment
## Kruskal-Wallis chi-squared = 10.817, df = 1, p-value = 0.001006kruskal.test(compareColor_mcap$Coral.Percent.White ~ compareColor_mcap$LifeStage)##
## Kruskal-Wallis rank sum test
##
## data: compareColor_mcap$Coral.Percent.White by compareColor_mcap$LifeStage
## Kruskal-Wallis chi-squared = 3.7231, df = 1, p-value = 0.05367color_ph2 = dunnTest(Coral.Percent.White ~ Treatment,
data=compareColor_mcap,
method="bonferroni")
color_ph2## Comparison Z P.unadj P.adj
## 1 Ambient - High -3.288968 0.001005556 0.001005556color_ph3 = dunnTest(Coral.Percent.White ~ LifeStage,
data=compareColor_mcap,
method="bonferroni")
color_ph3## Comparison Z P.unadj P.adj
## 1 Adult - Juvenile 1.929528 0.05366539 0.05366539In M. capitata, treatment is significant.
Pocillopora acuta:
modelB5b<-lmer(Coral.Percent.White~LifeStage * Treatment + (1|Colony), data=compareColor_pacu)
summary(modelB5b)
anova(modelB5b, type=2)
emm = emmeans(modelB5b, ~ Treatment, adjust="tukey")
cld(emm, Letters=c(LETTERS)) #letter display
pairs(emm)
emm = emmeans(modelB5b, ~ LifeStage, adjust="tukey")
cld(emm, Letters=c(LETTERS)) #letter display
pairs(emm)
emm = emmeans(modelB5b, ~ LifeStage * Treatment, adjust="tukey")
cld(emm, Letters=c(LETTERS)) #letter displayqqPlot(residuals(modelB5b))
compareColor_pacu$Group<-paste(compareColor_pacu$Treatment, compareColor_pacu$LifeStage)
bartlett.test(residuals(modelB5b)~compareColor_pacu$Group)Conduct non-parametric test.
kruskal.test(compareColor_pacu$Coral.Percent.White ~ compareColor_pacu$Treatment)##
## Kruskal-Wallis rank sum test
##
## data: compareColor_pacu$Coral.Percent.White by compareColor_pacu$Treatment
## Kruskal-Wallis chi-squared = 10.051, df = 1, p-value = 0.001522kruskal.test(compareColor_pacu$Coral.Percent.White ~ compareColor_pacu$LifeStage)##
## Kruskal-Wallis rank sum test
##
## data: compareColor_pacu$Coral.Percent.White by compareColor_pacu$LifeStage
## Kruskal-Wallis chi-squared = 0.93103, df = 1, p-value = 0.3346color_ph2 = dunnTest(Coral.Percent.White ~ Treatment,
data=compareColor_pacu,
method="bonferroni")
color_ph2## Comparison Z P.unadj P.adj
## 1 Ambient - High -3.17039 0.001522345 0.001522345color_ph3 = dunnTest(Coral.Percent.White ~ LifeStage,
data=compareColor_pacu,
method="bonferroni")
color_ph3## Comparison Z P.unadj P.adj
## 1 Adult - Juvenile 0.9649013 0.3345943 0.3345943There is an effect of treatment in P. acuta.
Coefficient of variation (CV):
sd(compareColor$Coral.Percent.White, na.rm=TRUE)/mean(compareColor$Coral.Percent.White, na.rm=TRUE)## [1] 0.3268894CVe <- ddply(compareColor, c("Treatment", "Species", "LifeStage"), summarise,
mean = mean(Coral.Percent.White, na.rm=TRUE),
sd = sd(Coral.Percent.White, na.rm=TRUE),
CV = (sd / mean)
)
knitr::kable(CVe)| Treatment | Species | LifeStage | mean | sd | CV |
|---|---|---|---|---|---|
| Ambient | Montipora capitata | Adult | 0.2823965 | 0.0595368 | 0.2108269 |
| Ambient | Montipora capitata | Juvenile | 0.2442736 | 0.0618264 | 0.2531031 |
| Ambient | Pocillopora acuta | Adult | 0.4279750 | 0.0822636 | 0.1922158 |
| Ambient | Pocillopora acuta | Juvenile | 0.4596533 | 0.0602377 | 0.1310503 |
| High | Montipora capitata | Adult | 0.4716676 | 0.1365140 | 0.2894284 |
| High | Montipora capitata | Juvenile | 0.3280025 | 0.0764757 | 0.2331559 |
| High | Pocillopora acuta | Adult | 0.5998219 | 0.0874831 | 0.1458484 |
| High | Pocillopora acuta | Juvenile | 0.5009551 | 0.0413132 | 0.0824688 |
D) Correlations Between Factors
Assemble correlation data.
#generate a datasheet with all values summarized by sample number at the final timepoint
corr_red <- bleach2 %>%
filter(Timepoint == "Final")
corr_thick <- bleach1 %>%
filter(Timepoint == "Final")
#remove unnecessary columns
corr_red<-corr_red[c("Sample", "Red")]
corr_thick<-corr_thick[c("Sample", "Thickness")]
corr_chl<-compareChl[c("Sample", "Chla.SA")]
corr_scope<-compareThick[c("Sample", "Thickness")]
corr_scope<-rename(corr_scope, c("Sample"="Sample", "Thickness"="Scope"))
corr_color<-compareColor[c("Sample", "Coral.Percent.White")]
correlations<-left_join(corr_red, corr_thick, by=c("Sample"))
correlations<-left_join(correlations, corr_scope, by=c("Sample"))
correlations<-left_join(correlations, corr_color, by=c("Sample"))
correlations<-left_join(correlations, corr_chl, by=c("Sample"))
#add in identifying information
correlations$Species<-bleach1$Species[match(correlations$Sample,bleach1$Sample)]
correlations$LifeStage<-bleach1$LifeStage[match(correlations$Sample,bleach1$Sample)]
correlations$Treatment<-bleach1$Treatment[match(correlations$Sample,bleach1$Sample)]Check for normality of data. Use a Spearman non-parametric test in analysis of dissecting microscope tissue thickness and chlorophyll data as they are not normal. Use a Pearson correlation for normally distributed confocal thickness, confocal fluorescence, and coral color.
shapiro.test(correlations$Chla.SA) #non normal##
## Shapiro-Wilk normality test
##
## data: correlations$Chla.SA
## W = 0.85088, p-value = 0.002823shapiro.test(correlations$Thickness) #normal##
## Shapiro-Wilk normality test
##
## data: correlations$Thickness
## W = 0.9736, p-value = 0.2457shapiro.test(correlations$Red) #normal##
## Shapiro-Wilk normality test
##
## data: correlations$Red
## W = 0.98918, p-value = 0.8905shapiro.test(correlations$Scope) #non normal##
## Shapiro-Wilk normality test
##
## data: correlations$Scope
## W = 0.87274, p-value = 0.00406shapiro.test(correlations$Coral.Percent.White) #normal##
## Shapiro-Wilk normality test
##
## data: correlations$Coral.Percent.White
## W = 0.97671, p-value = 0.3377####1. Conduct correlation of Chlorophyll and Symbiodiniaceae Fluorescence.
cor.test(correlations$Chla.SA, correlations$Red, method = "spearman")##
## Spearman's rank correlation rho
##
## data: correlations$Chla.SA and correlations$Red
## S = 670, p-value = 0.0006615
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6689723There is a significant correlation between Chlorophyll and Red values. Plot this relationship.
CorBleachRed<-ggplot(correlations, aes(x=Chla.SA, y=Red)) +
geom_point(aes(color=Treatment, shape=Species), size=3) +
scale_color_manual(name="Treatment",
values=c("blue", "red"))+
scale_shape_manual(name="Species",
values=c(19,21))+
theme_classic()+
ylim(0,60) +
geom_smooth(method=lm, se=FALSE, color="black", linetype="solid") +
theme(text = element_text(size = 18))+
theme(axis.text = element_text(size = 18, color="black"))+
theme(axis.title = element_text(size = 18, face="bold"))+
theme(legend.title = element_text(size = 18, face="bold"))+
theme(legend.text = element_text(size = 18, color="black", face="italic"))+
ylab("Symbiodiniaceae Fluorescence (% RI)") +
geom_text(x=0.22, y=2, fontface="plain", label="r=0.67, p<0.01", size=6) +
xlab(expression(bold(paste("Chlorophyll ", alpha," ("*mu*"g cm"^-2,")"))))
CorBleachRed
ggsave(filename="Figures/Fig2e.pdf", plot=CorBleachRed, dpi=300, width=10, height=8, units="in")####2. Conduct correlation of scope tissue thickness and confocal tissue thickness.
cor.test(correlations$Scope, correlations$Thickness, method = "spearman")##
## Spearman's rank correlation rho
##
## data: correlations$Scope and correlations$Thickness
## S = 308, p-value = 1.596e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8947009There is a significant correlation between Scope and confocal thickness values. Plot this relationship.
CorBleachThick<-ggplot(correlations, aes(x=Scope, y=Thickness)) +
geom_point(aes(color=Treatment, shape=Species), size=3) +
scale_color_manual(name="Treatment",
values=c("blue", "red"))+
scale_shape_manual(name="Species",
values=c(19,21))+
theme_classic()+
geom_smooth(method=lm, se=FALSE, color="black", linetype="solid") +
theme(text = element_text(size = 18))+
theme(axis.text = element_text(size = 18, color="black"))+
theme(axis.title = element_text(size = 18, face="bold"))+
theme(legend.title = element_text(size = 18, face="bold"))+
theme(legend.text = element_text(size = 18, color="black", face="italic"))+
ylab(expression(bold(paste("Depth of Fluorescence ("*mu*"m)")))) +
geom_text(x=1800, y=400, fontface="plain", label="r=0.89, p<0.01", size=6) +
xlab(expression(bold(paste("Tissue Thickness ("*mu*"m)"))))
CorBleachThick
ggsave(filename="Figures/Fig2d.pdf", plot=CorBleachThick, dpi=300, width=10, height=8, units="in")####3. Conduct correlation of confocal fluorescence and coral color.
cor.test(correlations$Red, correlations$Coral.Percent.White, method = "pearson") #significantly correlated##
## Pearson's product-moment correlation
##
## data: correlations$Red and correlations$Coral.Percent.White
## t = -5.1678, df = 55, p-value = 3.397e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7243689 -0.3656062
## sample estimates:
## cor
## -0.5717098There is a significant correlation between Scope and confocal thickness values. Plot this relationship. We expect an inverse relationship as proportion white indicates increasing paling.
CorBleachColor<-ggplot(correlations, aes(x=Coral.Percent.White, y=Red)) +
geom_point(aes(color=Treatment, shape=Species), size=3) +
scale_color_manual(name="Treatment",
values=c("blue", "red"))+
scale_shape_manual(name="Species",
values=c(19,21))+
theme_classic()+
ylim(0,65)+
geom_smooth(method=lm, se=FALSE, color="black", linetype="solid") +
theme(text = element_text(size = 18))+
theme(axis.text = element_text(size = 18, color="black"))+
theme(axis.title = element_text(size = 18, face="bold"))+
theme(legend.title = element_text(size = 18, face="bold"))+
theme(legend.text = element_text(size = 18, color="black", face="italic"))+
ylab(expression(bold(paste("Symbiodiniaceae Fluorescence (% RI)")))) +
geom_text(x=0.6, y=1, fontface="plain", label="r=-0.57, p<0.01", size=6) +
xlab(expression(bold(paste("Bleaching Score (Coral Color : White)"))))
CorBleachColor
ggsave(filename="Figures/Fig2f.pdf", plot=CorBleachColor, dpi=300, width=10, height=8, units="in")Conclusions:
(1) LSCM red values are significantly correlated with chlorophyll content.
(2) LSCM tissue thickness is significantly correlated to dissecting microscope tissue thickness.
(3) LSCM red values are significantly correlated to coral pigmentation.
(4) LSCM and comparative methodologies produce similar results, but with slight differences discussed in manuscript.