###########
#
#	PeerJ submission Code. Alex Young.  "Epiphyte type and sampling height".
#
# # #1. Read in dataset, calculate Simpsons diversity index. 

# # #2. Results for Microfauna density ANOVA

# # #3. Test for correlation of mass and microfauna density

# # #4. Figure 2. Scatter plots of microfauna Density. 

# # #5. Figure 3. NMDS on tardigrade community 

# # #6. Result for permanova (adonis) of tardigrade community composition 

# # #7. Indicator species analysis
#
###############################################################################

# #1. Read in dataset

###############################################################################

meio=read.csv("Raw_data_sheet_Young_06_22_18.csv   ", header=T)

#Calculate diversity indeces
library(vegan)
meio[1:4, 20:31] # These are the 12 species(taxa) of tardigrades.
meio$T.divSIM<-diversity(meio[ ,20:31], index="simpson")


meio[1:4,12:14] # these columns are the microfauna species richness
meio$M.divSIM<-diversity(meio[ ,12:14], index="simpson")


#############################################################

# #2. Results for microfauna Density (a blocked anova, with TukeyHSD)

#############################################################

# example blocking setup
#fit <- aov(y ~ A*B + block, data=mydataframe)  

head(meio)
# Nematode Density model
nem<- aov(nem.Density ~ position*epiphyte.type + Tree.ID, data=meio)
drop1(nem,~.,test="F")
TukeyHSD(nem)  

# Tardigrade Density model
tar<- aov(tardi.Density ~ position*epiphyte.type + Tree.ID, data=meio)
drop1(tar,~.,test="F")
TukeyHSD(tar)

# Rotifer Density model
rot<- aov(rot.Density ~ position*epiphyte.type + Tree.ID, data=meio)
drop1(rot,~.,test="F")
TukeyHSD(rot)

# Microfauna Richness model
M.richaov<-aov (Micro.Rich~ position*epiphyte.type + Tree.ID, data=meio)
drop1(M.richaov,~.,test="F")
TukeyHSD(M.richaov)

# Microfauna Diversity model
M.divaov<-aov (M.divSIM  ~ position*epiphyte.type + Tree.ID, data=meio)
drop1(M.divaov,~.,test="F")
TukeyHSD(M.divaov)

# Tar Diversity
T.div<-aov (T.divSIM  ~ position*epiphyte.type + Tree.ID, data=meio)
drop1(T.div,~.,test="F")
TukeyHSD(T.div)


##############################################################################

# #3. Correlation coefficient between mass and Density

##############################################################################

# nem Density ~ mass
cor.test(meio$nem.Density, meio$mass, use="pairwise.complete.obs", method="pearson", alternative="greater")$p.value

# rot Density~ mass
cor.test(meio$rot.Density, meio$mass, use="pairwise.complete.obs", method="pearson", alternative="greater")$p.value

# tardi Density ~ mass
cor.test(meio$tardi.Density, meio$mass, use="pairwise.complete.obs", method="pearson", alternative="greater")$p.value

# M.div ~ mass
cor.test(meio$M.divSW, meio$mass, use="pairwise.complete.obs", method="pearson", alternative="greater")$p.value

##############################################################################

# #4. Scatterplot of Microfauna Density   3x3 grid. Figure 2.

##############################################################################

library(tidyr)
library(ggplot2)

############  re-arrange dataset
micro<-gather(meio,"Microfauna","Density",12:14)
micro$epiphyte.type<-sapply(micro[, "epiphyte.type"], switch, "Bryophyte"="Bryophyte","Foliose"="Foliose lichen", "Fruticose"="Fruticose lichen")
micro$Microfauna<-sapply(micro[, "Microfauna"], switch, "nem.Density"="Nematode Density","rot.Density"="Rotifer Density","tardi.Density"="Tardigrade Density")

#Figure 2.
ggplot(micro, aes(x=Density, y=position, shape=position, col=epiphyte.type))+geom_jitter(height=.3, width=.2,size=2)+
  xlab("Density (animals/gram)")+ ylab("Canopy Sampling Height")+xlim(0,30)+
  facet_grid(epiphyte.type~ Microfauna)+theme_bw()+ scale_color_manual(values=c(' green','red','cyan'))+ 
  theme(panel.grid.minor=element_blank(), panel.grid.major=element_blank())+ theme(legend.position="none")


##############################################################################

# #5. NMDS on tardigrade community. Figure 3.

##############################################################################

library(vegan)

# For community matrices, we only want samples that had tardigrades in them. (where T.rich does not equal zero)
cm.only<-meio[!meio$T.Rich==0,20:31]
cm.sample<-meio[!meio$T.Rich==0,1:4]

#rownames(cm.only)<-cm.sample$sample.ID
pls<-cm.only
str(pls)

m1 <- metaMDS(pls, distance = "bray", k = 2, trymax = 200,
              autotransform =FALSE,  
              wascores = TRUE, expand = TRUE,
              trace = 1, plot = TRUE)

#### Figure 3
plot(m1, display = c("sites", "species"), xlim =c(-2,2),ylim=c(-1.8,1.8), choices = c(1,2),  type = "n", shrink = FALSE)
points(m1, display = "sites" , pch = c(15, 17, 16)[as.numeric(cm.sample$position)],col=c(3,2,5)[as.numeric(cm.sample$epiphyte.type)], cex=1.5)
legend("topleft", legend = paste(c("top sampling height", "mid sampling height", "low  sampling height")), pch = c(15, 17, 16))
legend("topright", legend = paste(c("Bryophytes", "Foliose Lichen", "Fruticose Lichen")),pch=8, col = c(3, 2, 5))
#ordispider( m1, cm.sample$position, display = "sites", kind = "sd", label = T)
ordiellipse(m1, cm.sample$epiphyte.type, display="sites", col=c(3,2,5),kind="sd", cex=1.5,lwd=2,label=T)


##############################################################################

# # 6. Permanova on community composition

##############################################################################

adonis(pls~cm.sample$epiphyte.type*cm.sample$position, strata=cm.sample$Tree.ID)

##############################################################################

# # 7. Indicator species analysis

##############################################################################

library(vegan)
library(indicspecies)

#test what species are indicators of epiphyte types, and canopy sampling heights.
indval.epi = multipatt(pls, cm.sample$epiphyte.type, control = how(nperm=999))
indval.pos = multipatt(pls, cm.sample$position, control = how(nperm=999))

summary(indval.epi, indvalcomp=TRUE) 
summary(indval.pos, indvalcomp=TRUE)