##Author: Jill R. Bourque, U.S. Geological Survey Wetland and Aquatic Research Center
##Gainesville, FL, 32653
##Contact: jbourque@usgs.gov

##R code for univariate statistics for PeeJ manuscript
##"The influence of different deep-sea coral habitats on sediment macrofaunal community structure and function"
##Revision date: 13-June-2018

##All required data to run below code can be downloaded from Supplemental Tables 1 & 2
##Definitions of column headers are provided on Metadata tab

##Full raw data can be downloaded at:
##https://doi.org/10.5066/F7TH8KXD

##First convert the data tabs in the supplemental tables to .csv files and save them locally. 

##Set your working directory to wherever you download the .csv files to

##Load the "car" library
library(car)

##1. Sediment infauna data analyses

##Import the infauna data
##Change file name to user defined parameters
infauna<-read.csv("Bourque_etal_GOM_Corals_Supp1.csv")

##Create data subsets for assumption testing
noloph<-infauna[c(which(!infauna$Coral=="Lophelia")),]    #All non-Lophelia cores
loph<-infauna[c(which(infauna$Coral=="Lophelia")),]       #Lophelia cores only
mad<-infauna[c(which(infauna$Coral=="Madrepora")),]       #Madrepora cores only
octo<-infauna[c(which(infauna$Coral=="Octocoral")),]      #Octocoral cores only
near<-infauna[c(which(infauna$Location=="Near")),]        #All near-coral cores only
back<-infauna[c(which(infauna$Location=="Background")),]  #All background cores only
lophnear<-loph[c(which(loph$Location=="Near")),]          #Lophelia near-coral cores only
madnear<-mad[c(which(mad$Location=="Near")),]             #Madrepora near-coral cores only
octonear<-octo[c(which(octo$Location=="Near")),]          #Octocoral near-coral cores only
lophback<-loph[c(which(loph$Location=="Background")),]    #Lophelia background cores only
madback<-mad[c(which(mad$Location=="Background")),]       #Madrepora background cores only
octoback<-octo[c(which(octo$Location=="Background")),]    #Octocoral background cores only
at<-infauna[c(which(infauna$Site=="AT47")),]              #Cores from site AT47 only
atnear<-at[c(which(at$Location=="Near")),]                #Near-coral cores at AT47 only
atback<-at[c(which(at$Location=="Background")),]          #Background cores at AT47 only
gb<-infauna[c(which(infauna$Site=="GB299")),]             #Cores from site GB299 only
gbnear<-gb[c(which(gb$Location=="Near")),]                #Near-coral cores at GB299 only
gbback<-gb[c(which(gb$Location=="Background")),]          #Background cores at GB299 only

##1a. Density statistics

##Between coral habitats

##Test assumptions
##Normality of each coral habitat
shapiro.test(lophnear$Density)       #Normal distribution
shapiro.test(madnear$Density)        #Significantly non-normal distribution
qqnorm(madnear$Density)              #Visual inspection of Madrepora density distribution
shapiro.test(octonear$Density)       #Normal distribution
##Test for homogeneity of variance
leveneTest(near$Density~near$Coral)  #Variances are homogeneous between coral habitats

##As ANOVA analysis is generally robust to small deviations in normality, subsequent analysis was performed
##on un-transformed data

##One-way ANOVA analysis of infaunal density among coral habitat types
anova(lm(near$Density~near$Coral))

##Post-hoc TUkey's Honest Squares Difference test for pairwise comparisons between coral habitat types
test<-aov(near$Density~near$Coral)
TukeyHSD(test)

##Relationship of density with depth

##Visual inspection of data distribution
plot(near$Density~near$Depth)
##Test for a polynomial relationship based on visual data inspection
summary(lm(near$Density~poly(near$Depth, 2, raw=TRUE)))

##1b. Shannon Diversity (H'loge) analyses

##Remove single large core from diversity analyses
infauna<-infauna[c(which(!infauna$CoreID=="SJ09-GOM-3725-PC05")),]

##Rerun lines 27 and 35 to exclude large core from data frames

##Test assumptions
##Normality of each coral habitat
shapiro.test(lophnear$Diversity)       #Normal distribution
shapiro.test(madnear$Diversity)        #Normal distribution
shapiro.test(octonear$Diversity)       #Significantly non-normal distribution
qqnorm(octonear$Diversity)             #Visual inspection of octocoral diversity distribution
##Test for homogeneity of variance
leveneTest(near$Diversity~near$Coral)  #Variances are homogeneous between coral habitats

##As ANOVA analysis is generally robust to small deviations in normality, subsequent analysis was performed
##on un-transformed data

##One-way ANOVA analysis of infaunal density among coral habitat types
anova(lm(near$Diversity~near$Coral))
##Post-hoc TUkey's Honest Squares Difference test for pairwise comparisons between coral habitat types
test<-aov(near$Diversity~near$Coral)
TukeyHSD(test)

##1c. Taxa Evenness (Pielou's J') analyses

##Test assumptions
##Normality of each coral habitat
shapiro.test(lophnear$Evenness)       #Normal distribution
shapiro.test(madnear$Evenness)        #Significantly non-normal distribution
shapiro.test(octonear$Evenness)       #Significantly non-normal distribution

##Test for homogeneity of variance
leveneTest(near$Evenness~near$Coral)  #Variances are homogeneous between coral habitats

##Since more than one group has a non-normal distribution, a non-parametric test was used. 
##Kruskal Wallis analysis of infaunal density among coral habitat types due to non-normal distributions
kruskal.test(near$Evenness~near$Coral)

##Post-hoc pairwise Wilcox test with Holm correction for comparisons between coral habitat types
pairwise.wilcox.test(near$Evenness, near$Coral, p.adjust.method="holm")

##1d. Comparison of diversity metrics with depth

##Shannon diversity (H'loge) 
shapiro.test(near$Diversity)         #Significantly non-normal distribution
shapiro.test(near$Depth)             #Significantly non-normal distribution

##Spearman correlation required due to non-normal distributions
cor.test(near$Diversity, near$Depth, method="spearman")

##Test of linear relationship of diversity and depth
summary(lm(near$Diversity~near$Depth))

##A polynomial relationship was tested based on known diversity patterns in the Gulf of Mexico
summary(lm(near$Diversity~poly(near$Depth, 2, raw=TRUE)))

##Taxa Evenness (Pielou's J')
shapiro.test(near$Evenness)         #Significantly non-normal distribution
shapiro.test(near$Depth)            #Significantly non-normal distribution

##Spearman correlation required due to non-normal distributions
cor.test(near$Evenness, near$Depth, method="spearman")

##Test of linear relationship of diversity and depth
summary(lm(near$Evenness~near$Depth))

##A polynomial relationship was tested based on known diversity patterns in the Gulf of Mexico
summary(lm(near$Evenness~poly(near$Depth, 2, raw=TRUE)))

##ES(n) rarefaction, where the value represents the estimated number of taxa that would be encountered
##when n number of individuals are sampled. 

##ES(5)
shapiro.test(near$ES5)           #Significantly non-normal distribution
shapiro.test(near$Depth)         #Significantly non-normal distribution

##Spearman correlation required due to non-normal distributions
cor.test(near$ES5, near$Depth, method="spearman")

##Test of linear relationship of diversity and depth
summary(lm(near$ES5~near$Depth))

##A polynomial relationship was tested based on known diversity patterns in the Gulf of Mexico
summary(lm(near$ES5~poly(near$Depth, 2, raw=TRUE)))

##ES(10)
shapiro.test(near$ES10)          #Significantly non-normal distribution
shapiro.test(near$Depth)         #Significantly non-normal distribution

##Spearman correlation required due to non-normal distributions
cor.test(near$ES10, near$Depth, method="spearman")

##Test of linear relationship of diversity and depth
summary(lm(near$ES10~near$Depth))

##A polynomial relationship was tested based on known diversity patterns in the Gulf of Mexico
summary(lm(near$ES10~poly(near$Depth, 2, raw=TRUE)))

##ES(20)
shapiro.test(near$ES20)          #Significantly non-normal distribution
shapiro.test(near$Depth)         #Significantly non-normal distribution

##Spearman correlation required due to non-normal distributions
cor.test(near$ES20, near$Depth, method="spearman")

##Test of linear relationship of diversity and depth
summary(lm(near$ES20~near$Depth))

##A polynomial relationship was tested based on known diversity patterns in the Gulf of Mexico
summary(lm(near$ES20~poly(near$Depth, 2, raw=TRUE)))

##2. Background habitats at coral sites

##Density

##Test for homogeneity of variance among background coral communities
leveneTest(back$Density~back$Coral)         #Significantly heterogeneous
##Test for normality
shapiro.test(lophback$Density)              #Significantly non-normal distribution
shapiro.test(madback$Density)   
shapiro.test(octoback$Density)

##Test for homogeneity of log-transformed density data
leveneTest(log(back$Density)~back$Coral)
##Test for normality of log-transformed data
shapiro.test(log(lophback$Density))
shapiro.test(log(madback$Density))   
shapiro.test(log(octoback$Density))

##Log-transformation of density data meets ANOVA assumptions

##One-way ANOVA analysis of background infaunal density among coral habitat types
anova(lm(log(back$Density)~back$Coral))
##Post-hoc TUkey's Honest Squares Difference test for pairwise comparisons between coral habitat types
test<-aov(log(back$Density)~back$Coral)
TukeyHSD(test)

##At AT47 only, the single Madrepora site with near and background cores

##Test for homogeneity of variance
leveneTest(at$Density~at$Location)
##Test for normality
shapiro.test(atnear$Density)
shapiro.test(atback$Density)
##One-way ANOVA analysis of infaunal density between locations
anova(lm(at$Density~at$Location))

##At GB299 only, the single octocoral site with near and background cores

##Test for homogeneity of variance
leveneTest(gb$Density~gb$Location)
##Test for normality
shapiro.test(gbnear$Density)
shapiro.test(gbback$Density)
##One-way ANOVA analysis of infaunal density between locations
anova(lm(gb$Density~gb$Location))

##Shannon diversity (H'loge)

##Test for homogeneity of variance among corals
leveneTest(back$Diversity~back$Coral)
##Test for normality
shapiro.test(lophback$Diversity)
shapiro.test(madback$Diversity)   
shapiro.test(octoback$Diversity)

##One-way ANOVA analysis of background infaunal diversity among coral habitat types
anova(lm(back$Diversity~back$Coral))
##Post-hoc TUkey's Honest Squares Difference test for pairwise comparisons between coral habitat types
test<-aov(back$Diversity~back$Coral)
TukeyHSD(test)

##At AT47 only, the single Madrepora site with near and background cores

##Test for homogeneity of variance
leveneTest(at$Diversity~at$Location)
##Test for normality
shapiro.test(atnear$Diversity)
shapiro.test(atback$Diversity)
##One-way ANOVA analysis of infaunal diversity between locations
anova(lm(at$Diversity~at$Location))

##At GB299 only, the single octocoral site with near and background cores

##Test for homogeneity of variance
leveneTest(gb$Diversity~gb$Location)
##Test for normality
shapiro.test(gbnear$Diversity)
shapiro.test(gbback$Diversity)
##One-way ANOVA analysis of infaunal diversity between locations
anova(lm(gb$Diversity~gb$Location))

##Comparing diversity between near and background cores across all Madrepora sites

##Test for homogeneity of variance
leveneTest(mad$Diversity~mad$Location)         #Significantly heterogeneous
##Test for normality
shapiro.test(madnear$Diversity)
shapiro.test(madback$Diversity)

##Kruskal Wallis analysis of diversity between locations due to heterogeneous variance
kruskal.test(mad$Diversity~mad$Location)

##Comparing diversity between near and background cores across all octocoral sites

##Test for homogeneity of variance
leveneTest(octo$Diversity~octo$Location)
##Test for normality
shapiro.test(octonear$Diversity)               #Significantly non-normal distribution
shapiro.test(octoback$Diversity)
qqnorm(octonear$Diversity)                     #Visual inspection of octocoral diversity distribution

##As ANOVA analysis is generally robust to small deviations in normality, subsequent analysis was performed
##on un-transformed data
anova(lm(octo$Diversity~octo$Location))

##3. Sediment geochemistry data analyses

##Import the sediment geochemistry data
##Change file name to user defined parameters
sedchem<-read.csv("Bourque_etal_GOM_Corals_Supp2.csv")

##Create data subsets for assumption testing
loph<-sedchem[c(which(sedchem$Habitat=="Lophelia")),]
mad<-sedchem[c(which(sedchem$Habitat=="Madrepora")),]  
octo<-sedchem[c(which(sedchem$Habitat=="Octocoral")),]

##Mud content

##Test for homogeneity of variance
leveneTest(sedchem$Mud~sedchem$Habitat)
##Test for normality
shapiro.test(loph$Mud)
shapiro.test(mad$Mud)
shapiro.test(octo$Mud)                        #Significantly non-normal distribution

##As ANOVA analysis is generally robust to small deviations in normality, subsequent analysis was performed
##on un-transformed data
##One-way ANOVA analysis of mud content among coral habitats
anova(lm(sedchem$Mud~sedchem$Habitat))
##Post-hoc TUkey's Honest Squares Difference test for pairwise comparisons between coral habitat types
test<-aov(sedchem$Mud~sedchem$Habitat)
TukeyHSD(test)

##Organic carbon content

##Test for homogeneity of variance
leveneTest(sedchem$OC~sedchem$Habitat)         #Significantly heterogeneous

##Kruskal Wallis analysis due to heterogeneous variance
kruskal.test(sedchem$OC~sedchem$Habitat)
##Post-hoc pairwise Wilcox test with Holm correction for comparisons between coral habitat types
pairwise.wilcox.test(sedchem$OC, sedchem$Habitat, p.adjust.method="holm")

##Relationship with depth

##Test for normality
shapiro.test(sedchem$OC)                        #Significantly non-normal distribution
shapiro.test(sedchem$Depth)                     #Significantly non-normal distribution

##Spearman correlation required due to non-normal distributions
cor.test(sedchem$OC, sedchem$Depth, method="spearman")