---
title: "Interactive effects of drought and deforestation on multitrophic communities and aquatic ecosystem functions in the Neotropics – a test using tank bromeliads"
author: "Marie Séguigne"
date: "2023-03-21"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#Packages to load
```{r}
library(SciViews)
library(ggplot2)
library(ggfortify)
library(tidyverse)
library(emmeans)
library(effects)
library(sjPlot)
library(multcomp)
library(patchwork)
library(dplyr)
library(rcompanion)
```

#I) Organise data set 
### I.1) Load data 
```{r}
mesh <- read.table("Data1_microorg.csv", header = TRUE, sep= ";", encoding="UTF-8")
mesh$Mesh.size <- as.factor(mesh$Mesh.size)
mesh$ID_plant <- as.character(mesh$ID_plant)
mesh <- dplyr::select(mesh, Mesh.size, ID_plant)

rawdata <- readxl::read_excel("C:/Users/mseguigne/Desktop/Interactif/Interactif_R/Data2_general.xlsx")
rawdata <- rawdata[,c(-5, -7:-9)]

#rawdata2[,2] <- lapply(rawdata2[,2], as.character)
rawdata[,c(1:4, 10)] <- lapply(rawdata[,c(1:4, 10)], as.factor)
rawdata[,c(5:9, 12)] <- lapply(rawdata[,c(5:9, 12)], as.numeric)
colnames(rawdata)[1] <- "Treatment"
colnames(rawdata)[5] <- "NL"
colnames(rawdata)[10] <- "Block"
colnames(rawdata)[11] <- "Vcol"
colnames(rawdata)[12] <- "Light"

test <- subset(rawdata, Habitat == "Forest")
test$Effects <- as.character(test$Effects)
test$Effects[test$Effects == "control"] <- "Forest control"
test$Effects[test$Effects == "deforestation"] <- "Deforestation"
test$Effects[test$Effects == "drought"] <- "Forest drought"
test$Effects[test$Effects == "deforestation + drought"] <- "Deforestation + Drought"
test$Effects[test$Effects == "deforestation_drought"] <- "Deforestation_Drought"


test2 <- subset(rawdata, Habitat == "Open")
test2$Effects <- as.character(test2$Effects)
test2$Effects[test2$Effects == "control"] <- "Open control"
test2$Effects[test2$Effects == "drought"] <- "Open drought"
test2$Effects <- as.factor(test2$Effects)

rawdata <- rbind(test,test2)
rm(test,test2)

rawdata <- merge(rawdata, mesh)
rawdata <- rawdata[,c(1:4,13,5:12)]
```

### I.2) Decomposition
```{r}
decomp <- read.csv("Data3_decomposition.csv", header = TRUE, sep =";", encoding = "UTF-8")
decomp <- decomp[c(-71:-999),]
colnames(decomp)[1] <- "Treatment"
decomp[,c(1,3,4)] <- lapply(decomp[,c(1,3,4)], as.factor)

#decomp rate 
T1 <- subset(decomp, Treatment == "T1")
T1$rate_fine <- (-ln(T1$final.dry.mass.fine.g/T1$initial.dry.mass.fine.g))/170
T1$rate_coarse <- (-ln(T1$final.dry.mass.coarse.g/T1$initial.dry.mass.coarse.g))/170

T2 <- subset(decomp, Treatment == "T2")
T2$rate_fine <- (-ln(T2$final.dry.mass.fine.g/T2$initial.dry.mass.fine.g))/170
T2$rate_coarse <- (-ln(T2$final.dry.mass.coarse.g/T2$initial.dry.mass.coarse.g))/170
decompo <- rbind(T1, T2)

T3 <- subset(decomp, Treatment == "T3")
T3$rate_fine <- (-ln(T3$final.dry.mass.fine.g/T3$initial.dry.mass.fine.g))/170
T3$rate_coarse <- (-ln(T3$final.dry.mass.coarse.g/T3$initial.dry.mass.coarse.g))/170
decompo <- rbind(decompo, T3)

T4 <- subset(decomp, Treatment == "T4")
T4$rate_fine <- (-ln(T4$final.dry.mass.fine.g/T4$initial.dry.mass.fine.g))/170
T4$rate_coarse <- (-ln(T4$final.dry.mass.coarse.g/T4$initial.dry.mass.coarse.g))/170
decompo <- rbind(decompo, T4)

T5 <- subset(decomp, Treatment == "T5")
T5$rate_fine <- (-ln(T5$final.dry.mass.fine.g/T5$initial.dry.mass.fine.g))/170
T5$rate_coarse <- (-ln(T5$final.dry.mass.coarse.g/T5$initial.dry.mass.coarse.g))/170
decompo <- rbind(decompo, T5)

T6 <- subset(decomp, Treatment == "T6")
T6$rate_fine <- (-ln(T6$final.dry.mass.fine.g/T6$initial.dry.mass.fine.g))/170
T6$rate_coarse <- (-ln(T6$final.dry.mass.coarse.g/T6$initial.dry.mass.coarse.g))/170
decompo <- rbind(decompo, T6)

T7 <- subset(decomp, Treatment == "T7")
T7$rate_fine <- (-ln(T7$final.dry.mass.fine.g/T7$initial.dry.mass.fine.g))/170
T7$rate_coarse <- (-ln(T7$final.dry.mass.coarse.g/T7$initial.dry.mass.coarse.g))/170
decompo <- rbind(decompo, T7)

rm(decomp, T1, T2, T3, T4, T5, T6, T7)

decompo <- decompo[,c(2, 11, 12)]
decompo$ID_plant <- as.character(decompo$ID_plant)
fine <- dplyr::select(decompo, ID_plant, rate_fine)
fine[3] <- c("Fine")
colnames(fine)[2] <- "decomp_rate"
colnames(fine)[3] <- "Mesh.size"
fine$Mesh.size <- as.factor(fine$Mesh.size)

coarse <- dplyr::select(decompo, ID_plant, rate_coarse)
coarse[3] <- c("Coarse")
colnames(coarse)[2] <- "decomp_rate"
colnames(coarse)[3] <- "Mesh.size"
coarse$Mesh.size <- as.factor(coarse$Mesh.size)

test <- rbind(fine, coarse)
rawdata <- merge(rawdata, test, by = c("ID_plant", "Mesh.size"))

rm(fine, coarse, test, mesh, decompo)
```

### I.3) Microorganisms
```{r}
microorg <- read.table("Data1_microorg.csv", header = TRUE, sep= ";", encoding="UTF-8")
microorg$Mesh.size <- as.factor(microorg$Mesh.size)
microorg$ID_plant <- as.character(microorg$ID_plant)
microorg[,c(8,10)] <- lapply(microorg[,c(8,10)], as.numeric)

fine <- subset(microorg, Mesh.size  == "Fine")
colnames(fine)[8] <- "Ergosterol"
colnames(fine)[10] <- "Attached_bact"
fine <- fine[,c(3,8, 10)]
rawdata_fine <- subset(rawdata, Mesh.size == "Fine")
test1 <- merge(rawdata_fine, fine)


coarse <- subset(microorg, Mesh.size  == "Coarse")
colnames(coarse)[8] <- "Ergosterol"
colnames(coarse)[10] <- "Attached_bact"
coarse <- coarse[,c(3,8, 10)]
rawdata_coarse <- subset(rawdata, Mesh.size == "Coarse")
test3 <- merge(rawdata_coarse, coarse)
rawdata <- rbind(test1, test3)

rm(microorg, fine, coarse, rawdata_coarse, rawdata_fine, test1, test3)
```

### I.4) Invertebrates biomass 
```{r}
inv_biomass <- read.csv("Data4_invertebrates.csv", header = FALSE, sep =";", encoding = "UTF-8")
inv_biomass <- inv_biomass[c(1:48),]
colnames(inv_biomass) <- inv_biomass[1,]
inv_biomass <- inv_biomass[-1,]
inv_biomass[,c(2:72)] <- lapply(inv_biomass[,c(2:72)], as.numeric)
inv_biomass <- inv_biomass[c(-1,-4, -6, -21:-23, -26, -27, -33, -40),] #take off outliers + species present in less than 2 brom
inv_biomass$grp <- c("collector", "collector", "collector", "filter feeder", "filter feeder","filter feeder","filter feeder",
                     "filter feeder","filter feeder","filter feeder","filter feeder", "mesopredator", "mesopredator",
                     "mesopredator", "collector", "collector","collector","collector", "mesopredator", "shredder", "shredder",
                     "shredder","shredder","shredder", "scraper", "scraper","scraper","scraper","scraper","scraper",
                     "predator","predator","predator","predator","predator","predator","predator")
inv_biomass$grp <- as.factor(inv_biomass$grp)
inv_biomass <- inv_biomass[,c(1,2,73, 3:72)]
colnames(inv_biomass)[1] <- "Sp_stage"
colnames(inv_biomass)[2] <- "biomass"

ID_plant <- dplyr::select(rawdata, ID_plant, Mesh.size)
ID_plant <- subset(rawdata, Mesh.size == "Fine")
ID_plant <- ID_plant$ID_plant

prey_biom <- inv_biomass[inv_biomass$grp %in% c("shredder", "collector", "filter feeder", "scraper"), ]
prey_biom$biom <- prey_biom[,2]*prey_biom[,4:73]
prey_biom <- prey_biom[,c(-1:-3)]
prey_biom <- apply(prey_biom, 2, sum)
prey_biom <- as.data.frame(prey_biom)
prey_biom <- prey_biom[c(71:140),]
prey_biom <- as.data.frame(prey_biom)
prey_biom <- cbind(ID_plant, prey_biom)
rawdata <- merge(rawdata, prey_biom, by = "ID_plant")
rm(prey_biom)

pred_biom <- inv_biomass[inv_biomass$grp %in% c("predator", "mesopredator"), ]
pred_biom$biom <- pred_biom[,2]*pred_biom[,4:73]
pred_biom <- pred_biom[,c(-1:-3)]
pred_biom <- apply(pred_biom, 2, sum)
pred_biom <- as.data.frame(pred_biom)
pred_biom <- pred_biom[c(71:140),]
pred_biom <- as.data.frame(pred_biom)
pred_biom <- cbind(ID_plant, pred_biom)
rawdata <- merge(rawdata, pred_biom, by = "ID_plant")
rm(pred_biom)

shred_biom <- inv_biomass[inv_biomass$grp %in% c("shredder"),]
shred_biom$biom <- shred_biom[,2]*shred_biom[,4:73]
shred_biom <- shred_biom[,c(-1:-3)]
shred_biom <- apply(shred_biom, 2, sum)
shred_biom <- as.data.frame(shred_biom)
shred_biom <- shred_biom[c(71:140),]
shred_biom <- as.data.frame(shred_biom)
shred_biom <- cbind(ID_plant, shred_biom)
rawdata <- merge(rawdata, shred_biom, by = "ID_plant")
rm(shred_biom)

scrap_biom <- inv_biomass[inv_biomass$grp %in% c("scraper"),]
scrap_biom$biom <- scrap_biom[,2]*scrap_biom[,4:73]
scrap_biom <- scrap_biom[,c(-1:-3)]
scrap_biom <- apply(scrap_biom, 2, sum)
scrap_biom <- as.data.frame(scrap_biom)
scrap_biom <- scrap_biom[c(71:140),]
scrap_biom <- as.data.frame(scrap_biom)
scrap_biom <- cbind(ID_plant, scrap_biom)
rawdata <- merge(rawdata, scrap_biom, by = "ID_plant")
rm(scrap_biom)

fil_biom <- inv_biomass[inv_biomass$grp %in% c("filter feeder"),]
fil_biom$biom <- fil_biom[,2]*fil_biom[,4:73]
fil_biom <- fil_biom[,c(-1:-3)]
fil_biom <- apply(fil_biom, 2, sum)
fil_biom <- as.data.frame(fil_biom)
fil_biom <- fil_biom[c(71:140),]
fil_biom <- as.data.frame(fil_biom)
fil_biom <- cbind(ID_plant, fil_biom)
rawdata <- merge(rawdata, fil_biom, by = "ID_plant")
rm(fil_biom)

col_biom <- inv_biomass[inv_biomass$grp %in% c("collector"),]
col_biom$biom <- col_biom[,2]*col_biom[,4:73]
col_biom <- col_biom[,c(-1:-3)]
col_biom <- apply(col_biom, 2, sum)
col_biom <- as.data.frame(col_biom)
col_biom <- col_biom[c(71:140),]
col_biom <- as.data.frame(col_biom)
col_biom <- cbind(ID_plant, col_biom)
rawdata <- merge(rawdata, col_biom, by = "ID_plant")
rm(col_biom)

benthic_biom <- inv_biomass[inv_biomass$grp %in% c("shredder", "collector", "scraper"),]
benthic_biom$biom <- benthic_biom[,2]*benthic_biom[,4:73]
benthic_biom <- benthic_biom[,c(-1:-3)]
benthic_biom <- apply(benthic_biom, 2, sum)
benthic_biom <- as.data.frame(benthic_biom)
benthic_biom <- benthic_biom[c(71:140),]
benthic_biom <- as.data.frame(benthic_biom)
benthic_biom <- cbind(ID_plant, benthic_biom)
rawdata <- merge(rawdata, benthic_biom, by = "ID_plant")
rm(benthic_biom)

FPOM_biom <- inv_biomass[inv_biomass$grp %in% c("filter feeder", "collector"),]
FPOM_biom$biom <- FPOM_biom[,2]*FPOM_biom[,4:73]
FPOM_biom <- FPOM_biom[,c(-1:-3)]
FPOM_biom <- apply(FPOM_biom, 2, sum)
FPOM_biom <- as.data.frame(FPOM_biom)
FPOM_biom <- FPOM_biom[c(71:140),]
FPOM_biom <- as.data.frame(FPOM_biom)
FPOM_biom <- cbind(ID_plant, FPOM_biom)
rawdata <- merge(rawdata, FPOM_biom, by = "ID_plant")
rm(FPOM_biom)

CPOM_biom <- inv_biomass[inv_biomass$grp %in% c("shredder", "scraper"),]
CPOM_biom$biom <- CPOM_biom[,2]*CPOM_biom[,4:73]
CPOM_biom <- CPOM_biom[,c(-1:-3)]
CPOM_biom <- apply(CPOM_biom, 2, sum)
CPOM_biom <- as.data.frame(CPOM_biom)
CPOM_biom <- CPOM_biom[c(71:140),]
CPOM_biom <- as.data.frame(CPOM_biom)
CPOM_biom <- cbind(ID_plant, CPOM_biom)
rawdata <- merge(rawdata, CPOM_biom, by = "ID_plant")
rm(CPOM_biom)
```

### I.5) Microorganisms biomass
```{r}
#Estimating bacteria biomass
bact_biovol = (4/3)*pi*(0.4^3) #We consider a bacteria as a sphere with a diameter = 0.8µm, so a radius of 0.4µm
rawdata$bact_biomass <- ((rawdata$Attached_bact)*bact_biovol*(4*10^(-7)))/10 #To express in µgDW
#rawdata2_count$bact_biomass <- ((rawdata2_count$Attached_bact)*bact_biovol*0.4)/10 #To express in pgDW
rm(bact_biovol)
```

### I.6) Final data frames
```{r}
data_mesh <- rawdata
data_mesh$Effects <- as.factor(data_mesh$Effects)
data_mesh$Effects_f <- factor(data_mesh$Effects, levels=c('Forest control', 'Deforestation','Forest drought','Deforestation + Drought', 'Deforestation_Drought', 'Open control', 'Open drought'))

#Dataframe for invertebrates and microorganisms 
data <- subset(data_mesh, Mesh.size == "Fine")

control <- data[data$Treatment %in% c("T1", "T6"), ]
control <- droplevels(control)

data_for <- data[!data$Treatment %in% c("T6", "T7"), ]
data_for <- droplevels(data_for)

data_open <- data[data$Treatment %in% c("T6", "T7"), ]
data_open <- droplevels(data_open)

#Dataframe for decomposition and attached microorganisms 
data_mesh <- data_mesh[c(-5,-11,-18,-19, -37,-45, -49, -52, -53, -85, -90, -54),] #clean dataset 

control_mesh <- data_mesh[data_mesh$Treatment %in% c("T1", "T6"), ]
control_mesh <- droplevels(control_mesh)

data_for_mesh <- data_mesh[!data_mesh$Treatment %in% c("T6", "T7"), ]
data_for_mesh <- droplevels(data_for_mesh)

data_open_mesh <- data_mesh[data_mesh$Treatment %in% c("T6", "T7"), ]
data_open_mesh <- droplevels(data_open_mesh)
```



#II) Habitat conditions
Water volume
```{r}
# Mean and SE 
Rmisc::summarySE(data = data, measurevar = "Vcol", groupvars = "Treatment")

#Aligned rank transformed Anova
mod<- ARTool::art(Vcol  ~ Effects_f, data = data)
summary(mod) #ok 
anova(mod)
marginal = ARTool::art.con(mod, "Effects_f", adjust="bonferroni")
Sum = as.data.frame(marginal)
letters  <- cldList(p.value ~ contrast, data=Sum)
myletters_df <- data.frame(Effects_f=levels(data$Effects_f),letters=letters)

#Plot
p1 <- ggplot(data, aes(x = Effects_f, y = Vcol, color = Treatment)) +
      geom_boxplot(outlier.alpha = 0, alpha=0.5) +
      geom_jitter(width=0.25)+
        theme(axis.text.x = element_blank(), 
              axis.title.y = element_text(vjust = 2, size = 12, face = "bold"),
              legend.position="none", 
              strip.text.x = element_text(face = "bold")) +
          geom_text(data = myletters_df, aes(label = letters.Letter, y = 0 + 3020), colour="black", size= 4, fontface = "bold") + 
  ylim(0,3050) + 
  labs(x = "", y = "Water volume collected (mL)") 

```

Light input
```{r}
##mean and SE
Rmisc::summarySE(data = data, measurevar = "Light", groupvars = "Treatment")

#Aligned rank transformed Anova
mod<- ARTool::art(Light  ~ Effects_f, data = data)
summary(mod) #ok 
anova(mod)
marginal = ARTool::art.con(mod, "Effects_f", adjust="bonferroni")
Sum = as.data.frame(marginal)
letters  <- cldList(p.value ~ contrast, data=Sum)
myletters_df <- data.frame(Effects_f=levels(data$Effects_f),letters=letters)

#Plot
p2 <- ggplot(data, aes(x = Effects_f, y = Light, color = Treatment)) +
      geom_boxplot(outlier.alpha = 0, alpha=0.5) +
      geom_jitter(width=0.25)+
        theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1, size = 10, face = "bold"), 
              axis.title.y = element_text(vjust = 2, size = 12, face = "bold"),
              legend.position="none", 
              strip.text.x = element_text(face = "bold")) +
          geom_text(data = myletters_df, aes(label = letters.Letter, y = 0 + 84), colour="black", size= 4, fontface = "bold") + 
  ylim(0,85) + 
  labs(x = "", y = "Light input (%)") 
```

Figure 2
```{r}
theme_set(
  theme_bw() +
    theme(legend.position = "right")
  )
p <- p1 / p2 +  plot_layout(guides = "collect")

ggsave("Figure2.png", p, width = 20, height = 20, units = "cm")

rm(mod, marginal, Sum, letters, myletters_df, p1, p2, p)
```


#II) Predator and prey biomass
Comparison of both controls 
```{r}
#Non parametric test --> small sample size 
wilcox.test(control$pred_biom ~ control$Treatment,paired = FALSE, detailed = TRUE) 
wilcox.test(control$prey_biom ~ control$Treatment,paired = FALSE, detailed = TRUE) 
```

Impact of treatments in forest 
```{r}
#Predators
mod <- ARTool::art(pred_biom  ~ Effects_f, data = data_for)
summary(mod) #ok 
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")


#Prey
mod<- ARTool::art(prey_biom  ~ Effects_f, data = data_for)
summary(mod) #ok 
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")

rm(mod)
```

Impact of drought in open area
```{r}
#Non parametric test --> assumptions not verified + small sample size + unbalanced design 
wilcox.test(data_open$pred_biom~ data_open$Treatment,paired = FALSE, detailed = TRUE) 
wilcox.test(data_open$prey_biom~ data_open$Treatment,paired = FALSE, detailed = TRUE) 
```



#III) FPOM (col + fil) and CPOM (shred + scrap) consumers biomass
Comparison of both controls 
```{r}
#Non parametric test --> small sample size 
wilcox.test(control$FPOM_biom ~ control$Treatment,paired = FALSE, detailed = TRUE) 
wilcox.test(control$CPOM_biom ~ control$Treatment,paired = FALSE, detailed = TRUE) 
```

Impact of treatments in forest 
```{r}
#FPOM consumers
mod <- ARTool::art(FPOM_biom  ~ Effects_f, data = data_for)
summary(mod) #ok 
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")

rm(mod)

#CPOM consumers
mod <- ARTool::art(CPOM_biom  ~ Effects_f, data = data_for)
summary(mod) #ok 
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")

rm(mod)
```

Impact of drought in open area 
```{r}
#Non parametric test --> small sample size 
wilcox.test(data_open$FPOM_biom~ data_open$Treatment,paired = FALSE, detailed = TRUE) 
wilcox.test(data_open$CPOM_biom~ data_open$Treatment,paired = FALSE, detailed = TRUE) 
```


### IV.4) FPOM and CPOM comparison 
```{r}
#Control forest 
CF <- subset(data, Treatment == "T1")
shapiro.test(CF$FPOM_biom)
shapiro.test(CF$CPOM_biom)

test <- tidyr::pivot_longer(CF, cols = c(35,36), names_to ="FPOM_CPOM",
values_to = "biom")
test$FPOM_CPOM <- as.factor(test$FPOM_CPOM)

wilcox.test(test$biom ~ test$FPOM_CPOM,paired = FALSE, detailed = TRUE) 

#drought
CF <- subset(data, Treatment == "T2")
shapiro.test(CF$FPOM_biom)
shapiro.test(CF$CPOM_biom)

test <- tidyr::pivot_longer(CF, cols = c(35,36), names_to ="FPOM_CPOM",
values_to = "biom")
test$FPOM_CPOM <- as.factor(test$FPOM_CPOM)

wilcox.test(test$biom ~ test$FPOM_CPOM,paired = FALSE, detailed = TRUE) 


#deforestation
CF <- subset(data, Treatment == "T3")
shapiro.test(CF$FPOM_biom)
shapiro.test(CF$CPOM_biom)

test <- tidyr::pivot_longer(CF, cols = c(35,36), names_to ="FPOM_CPOM",
values_to = "biom")
test$FPOM_CPOM <- as.factor(test$FPOM_CPOM)

wilcox.test(test$biom ~ test$FPOM_CPOM,paired = FALSE, detailed = TRUE) 


#deforestation_drought
CF <- subset(data, Treatment == "T5")
shapiro.test(CF$FPOM_biom)
shapiro.test(CF$CPOM_biom)

test <- tidyr::pivot_longer(CF, cols = c(35,36), names_to ="FPOM_CPOM",
values_to = "biom")
test$FPOM_CPOM <- as.factor(test$FPOM_CPOM)

wilcox.test(test$biom ~ test$FPOM_CPOM,paired = FALSE, detailed = TRUE) 

#control open
CF <- subset(data, Treatment == "T6")
shapiro.test(CF$FPOM_biom)
shapiro.test(CF$CPOM_biom)

test <- tidyr::pivot_longer(CF, cols = c(35,36), names_to ="FPOM_CPOM",
values_to = "biom")
test$FPOM_CPOM <- as.factor(test$FPOM_CPOM)

wilcox.test(test$biom ~ test$FPOM_CPOM,paired = FALSE, detailed = TRUE) 

#drought open 
CF <- subset(data, Treatment == "T7")
shapiro.test(CF$FPOM_biom)
shapiro.test(CF$CPOM_biom)

test <- tidyr::pivot_longer(CF, cols = c(35,36), names_to ="FPOM_CPOM",
values_to = "biom")
test$FPOM_CPOM <- as.factor(test$FPOM_CPOM)

wilcox.test(test$biom ~ test$FPOM_CPOM,paired = FALSE, detailed = TRUE) 
```

#IV) Bacteria and fungi biomass 
Comparison of both controls 
```{r}
##Attached bacteria biomass
#Aligned rank transformed Anova with Mesh size and effects as factor
mod <- ARTool::art(bact_biomass  ~ Effects_f + Mesh.size + Effects_f:Mesh.size, data = control_mesh)
summary(mod) #Not at 0, let's subset dataset by the type of mesh

test <- subset(control_mesh, Mesh.size =="Fine")
wilcox.test(test$bact_biomass ~ test$Treatment,paired = FALSE, detailed = TRUE) 
test <- subset(control_mesh, Mesh.size =="Coarse")
wilcox.test(test$bact_biomass ~ test$Treatment,paired = FALSE, detailed = TRUE) 

rm(mod)

##Attached fungi biomass
#Aligned rank transformed Anova with Mesh size and effects as factor
mod <- ARTool::art(Ergosterol  ~ Effects_f + Mesh.size + Effects_f:Mesh.size, data = control_mesh)
summary(mod) #Not at 0, let's subset dataset by the type of mesh

test <- subset(control_mesh, Mesh.size =="Fine")
wilcox.test(test$Ergosterol ~ test$Treatment,paired = FALSE, detailed = TRUE) 
test <- subset(control_mesh, Mesh.size =="Coarse")
wilcox.test(test$Ergosterol ~ test$Treatment,paired = FALSE, detailed = TRUE) 

rm(mod)
```

Impact of treatments in forest 
```{r}
##Attacged bacteria
#Aligned rank transformed Anova with Mesh size and effects as factor
mod <- ARTool::art(bact_biomass  ~ Effects_f + Mesh.size + Effects_f:Mesh.size, data = data_for_mesh)
summary(mod) #Not at 0, let's subset dataset by the type of mesh

test <- subset(data_for_mesh, Mesh.size =="Fine")
mod <- ARTool::art(bact_biomass  ~ Effects_f, data = test)
summary(mod) #ok 
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")


test <- subset(data_for_mesh, Mesh.size =="Coarse")
mod <- ARTool::art(bact_biomass  ~ Effects_f, data = test)
summary(mod) #ok 
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")

rm(mod, test)

##Attacged fungi
#Aligned rank transformed Anova with Mesh size and effects as factor
mod <- ARTool::art(Ergosterol  ~ Effects_f + Mesh.size + Effects_f:Mesh.size, data = data_for_mesh)
summary(mod) #Not at 0, let's subset dataset by the type of mesh

test <- subset(data_for_mesh, Mesh.size =="Fine")
mod <- ARTool::art(Ergosterol  ~ Effects_f, data = test)
summary(mod) #ok 
anova(mod)

test <- subset(data_for_mesh, Mesh.size =="Coarse")
mod <- ARTool::art(Ergosterol  ~ Effects_f, data = test)
summary(mod) #ok 
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")

rm(mod, test)
```

Impact of drought in open area 
```{r}
##Attacged bacteria
test <- subset(data_open_mesh, Mesh.size =="Fine")
wilcox.test(test$bact_biomass ~ test$Treatment,paired = FALSE, detailed = TRUE) 
test <- subset(data_open_mesh, Mesh.size =="Coarse")
wilcox.test(test$bact_biomass ~ test$Treatment,paired = FALSE, detailed = TRUE) 

rm(test)

#Attached fungi
test <- subset(data_open_mesh, Mesh.size =="Fine")
wilcox.test(test$Ergosterol ~ test$Treatment,paired = FALSE, detailed = TRUE) 
test <- subset(data_open_mesh, Mesh.size =="Coarse")
wilcox.test(test$Ergosterol ~ test$Treatment,paired = FALSE, detailed = TRUE) 

rm(test)
```



#V) Decomposition rate
Comparison of both controls 
```{r}
#Aligned rank transformed Anova with Mesh size and effects as factor
mod <- ARTool::art(decomp_rate  ~ Habitat + Mesh.size + Habitat:Mesh.size, data = control_mesh)
summary(mod)#Not at 0, let's subset dataset by the type of mesh

test <- subset(control_mesh, Mesh.size =="Fine")
wilcox.test(test$decomp_rate ~ test$Treatment,paired = FALSE, detailed = TRUE) 
test <- subset(control_mesh, Mesh.size =="Coarse")
wilcox.test(test$decomp_rate ~ test$Treatment,paired = FALSE, detailed = TRUE) 

rm(mod, test)
```

Impact of treatments in forest 
```{r}
test <- subset(data_for_mesh, Mesh.size =="Fine")
mod <-  ARTool::art(decomp_rate ~ Effects_f, data = test)
summary(mod) #ok
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")


test <- subset(data_for_mesh, Mesh.size =="Coarse")
mod <-  ARTool::art(decomp_rate ~ Effects_f, data = test)
summary(mod)#OK
anova(mod)
contrast(emmeans(ARTool::artlm(mod, "Effects_f"), ~ Effects_f), method = "dunnett",infer = TRUE, adjust = "bonferroni")


rm(mod, test)
```

Impact of drought in open area 
```{r}
test <- subset(data_open_mesh, Mesh.size =="Fine")
wilcox.test(test$decomp_rate ~ test$Treatment,paired = FALSE, detailed = TRUE) 
test <- subset(data_open_mesh, Mesh.size =="Coarse")
wilcox.test(test$decomp_rate ~ test$Treatment,paired = FALSE, detailed = TRUE) 

rm(test)
```


#VI) Figure 3
```{r}
test <- dplyr::select(data, ID_plant, Effects, prey_biom, pred_biom)
test$Effects_f <- factor(test$Effects, levels=c('Open drought', 'Open control', 'Deforestation_Drought', 'Deforestation + Drought', 'Forest drought','Deforestation','Forest control'))

prey <- dplyr::select(test, ID_plant, Effects_f, prey_biom)
prey$group <- "Prey"
colnames(prey)[3] <- "Biomass"

pred <- dplyr::select(test, ID_plant, Effects_f, pred_biom)
pred$group <- "Predator"
colnames(pred)[3] <- "Biomass"

test <- rbind(prey, pred)
test$group <- as.factor(test$group)

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
 return(data_sum)
}
df2 <- data_summary(test, varname="Biomass", 
                    groupnames=c("group", "Effects_f"))

df2$group_f=factor(df2$group, levels = c("Prey", "Predator"))


pred_prey <- ggplot(df2, aes(x=Effects_f, y=Biomass, fill=group_f)) +  labs(fill = "") + ylab("Biomass (DW/mg)") +
  geom_bar(stat="identity", color = "black", position=position_dodge(),  alpha=0.8) + 
  scale_fill_manual(values = c("#009E73", "#D55E00")) + 
  geom_errorbar(aes(ymin=Biomass, ymax= Biomass+sd), width=.2,
                 position=position_dodge(.9)) + 
  theme(axis.title.y = element_blank(), 
        axis.title.x = element_text(face = "bold", size = 14), 
        axis.text.y = element_text(face = "bold", size = 14),
        axis.text.x = element_text(face = "bold", size = 14),
        legend.text = element_text(face = "bold", size = 14),
        legend.position = "top", 
        plot.title = element_text(face = "bold"))  + 
  coord_flip() +
  ggtitle('(a)')


rm(test, pred, prey, df2, data_summary)

test <- dplyr::select(data, ID_plant, Effects, FPOM_biom, CPOM_biom)
test$Effects_f <- factor(test$Effects, levels=c('Open drought', 'Open control', 'Deforestation_Drought', 'Deforestation + Drought', 'Forest drought','Deforestation','Forest control'))

FPOM <- dplyr::select(test, ID_plant, Effects_f, FPOM_biom)
FPOM$group <- "FPOM c"
colnames(FPOM)[3] <- "Biomass"

CPOM <- dplyr::select(test, ID_plant, Effects_f, CPOM_biom)
CPOM$group <- "CPOM c"
colnames(CPOM)[3] <- "Biomass"

test <- rbind(FPOM, CPOM)
test$group <- as.factor(test$group)

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
 return(data_sum)
}

df2 <- data_summary(test, varname="Biomass", 
                    groupnames=c("group", "Effects_f"))

df2$group_f=factor(df2$group, levels = c("FPOM c", "CPOM c"))


CPOM_FPOM <- ggplot(df2, aes(x=Effects_f, y=Biomass, fill=group_f)) +  labs(fill = "") + ylab("Biomass (DW/mg)") +
  geom_bar(stat="identity", color = "black", position=position_dodge(),  alpha=0.8) + 
  scale_fill_manual(values = c("#009292", "#FFC34C")) + 
  geom_errorbar(aes(ymin=Biomass, ymax= Biomass+sd), width=.2,
                 position=position_dodge(.9)) + 
  theme(axis.title.y = element_blank(), 
        axis.title.x = element_text(face = "bold", size = 14), 
        axis.text.y = element_text(face = "bold", size = 14),
        axis.text.x = element_text(face = "bold", size = 14),
        legend.text = element_text(face = "bold", size = 14),
        legend.position = "top", 
        plot.title = element_text(face = "bold"))  + 
  coord_flip() +
  ggtitle('(b)')

rm(test, CPOM, FPOM, df2, data_summary)

test <- dplyr::select(data_mesh, ID_plant, Effects, bact_biomass, Ergosterol, Mesh.size)
test$Effects_f <- factor(test$Effects, levels=c('Open drought', 'Open control', 'Deforestation_Drought', 'Deforestation + Drought', 'Forest drought','Deforestation','Forest control'))

bact <- dplyr::select(test, ID_plant, Effects_f, bact_biomass, Mesh.size)
bact$group <- "Bacteria"
colnames(bact)[3] <- "Biomass"

fungi <- dplyr::select(test, ID_plant, Effects_f, Ergosterol, Mesh.size)
fungi$group <- "Fungi"
colnames(fungi)[3] <- "Biomass"

test <- rbind(bact, fungi)
test$group <- as.factor(test$group)

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
 return(data_sum)
}
df2 <- data_summary(test, varname="Biomass", 
                    groupnames=c("group", "Effects_f", "Mesh.size"))

df2$group_f=factor(df2$group, levels = c("Bacteria", "Fungi"))


bact_fung <- ggplot(df2, aes(x=Effects_f, y=log(Biomass +1), fill=group_f)) +  labs(fill = "") + ylab(expression(bold(atop("Log(Biomass +1)", paste("(B:μgDW/gDW; F:μg/gDW)"))))) +
  geom_bar(stat="identity", color = "black", position=position_dodge(),  alpha=0.8) + 
  scale_fill_manual(values = c("#56B4E9", "grey")) + 
  facet_grid(rows = vars(Mesh.size)) +
  geom_errorbar(aes(ymin=log(Biomass+1), ymax= log(Biomass +sd+1)), width=.2, position=position_dodge(.9)) + 
  theme(axis.title.y = element_blank(), 
        axis.title.x = element_text(face = "bold", size = 14), 
        axis.text.y = element_text(face = "bold", size = 14),
        axis.text.x = element_text(face = "bold", size = 14),
        legend.text = element_text(face = "bold", size = 14), 
        legend.position = "top", 
        plot.title = element_text(face = "bold"), 
        strip.text = element_text(size = 14))  + 
  coord_flip() +
  ggtitle('(c)')


rm(test, fungi, bact, df2, data_summary)

test <- dplyr::select(data_mesh, ID_plant, Effects, decomp_rate, Mesh.size)
test$Effects_f <- factor(test$Effects, levels=c('Open drought', 'Open control', 'Deforestation_Drought', 'Deforestation + Drought', 'Forest drought','Deforestation','Forest control'))

decomp <- ggplot(test, aes(x=Effects_f, y=decomp_rate, fill=Mesh.size)) + 
  geom_boxplot() +  
  labs(fill = "") +
  scale_fill_manual(values = c("orange4", "wheat"))  + 
 theme(axis.title.y = element_blank(), 
        axis.title.x = element_text(face = "bold", size = 14), 
        axis.text.y = element_text(face = "bold", size = 14),
        axis.text.x = element_text(face = "bold", size = 14, angle = 20, vjust = 1, hjust = 1),
        legend.text = element_text(face = "bold", size = 14), 
        legend.position = "top", 
        plot.title = element_text(face = "bold"))+ 
  coord_flip() + 
  ylab(expression(bold(paste("Decomposition rate k", (d^{-1})))))+
  ggtitle('(d)')


library(patchwork)
theme_set(
  theme_bw() +
    theme(legend.position = "top")
  )

p <- pred_prey + CPOM_FPOM + bact_fung + decomp 
ggsave("figure3.png", p, width = 25, height = 25, units = "cm")
rm(pred_prey, CPOM_FPOM, bact_fung, decomp, p)
```

#VII) Figures S1 and S2
Create a new general df 
```{r}
coreth_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Corethrella.sp_st2", "Corethrella.sp_st3", "Corethrella.sp_st4"), ]
coreth_biom$biom <- coreth_biom[,2]*coreth_biom[,4:73]
coreth_biom <- coreth_biom[,c(-1:-3)]
coreth_biom <- apply(coreth_biom, 2, sum)
coreth_biom <- as.data.frame(coreth_biom)
coreth_biom <- coreth_biom[c(71:140),]
coreth_biom <- as.data.frame(coreth_biom)
coreth_biom <- cbind(ID_plant, coreth_biom)
test <- merge(data, coreth_biom, by = "ID_plant")
rm(coreth_biom)

Tanyp_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Tanypodinae.sp"), ]
Tanyp_biom$biom <- Tanyp_biom[,2]*Tanyp_biom[,4:73]
Tanyp_biom <- Tanyp_biom[,c(-1:-3)]
Tanyp_biom <- apply(Tanyp_biom, 2, sum)
Tanyp_biom <- as.data.frame(Tanyp_biom)
Tanyp_biom <- Tanyp_biom[c(71:140),]
Tanyp_biom <- as.data.frame(Tanyp_biom)
Tanyp_biom <- cbind(ID_plant, Tanyp_biom)
test <- merge(test, Tanyp_biom, by = "ID_plant")
rm(Tanyp_biom)

Eryth_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Erythrodiplax_small", "Erythrodiplax_medium", "Erythrodiplax_large"), ]
Eryth_biom$biom <- Eryth_biom[,2]*Eryth_biom[,4:73]
Eryth_biom <- Eryth_biom[,c(-1:-3)]
Eryth_biom <- apply(Eryth_biom, 2, sum)
Eryth_biom <- as.data.frame(Eryth_biom)
Eryth_biom <- Eryth_biom[c(71:140),]
Eryth_biom <- as.data.frame(Eryth_biom)
Eryth_biom <- cbind(ID_plant, Eryth_biom)
test <- merge(test, Eryth_biom, by = "ID_plant")
rm(Eryth_biom)

Eryth_small <- inv_biomass[inv_biomass$Sp_stage %in% c("Erythrodiplax_small"), ]
Eryth_small$biom <- Eryth_small[,2]*Eryth_small[,4:73]
Eryth_small <- Eryth_small[,c(-1:-3)]
Eryth_small <- apply(Eryth_small, 2, sum)
Eryth_small <- as.data.frame(Eryth_small)
Eryth_small <- Eryth_small[c(71:140),]
Eryth_small <- as.data.frame(Eryth_small)
Eryth_small <- cbind(ID_plant, Eryth_small)
test <- merge(test, Eryth_small, by = "ID_plant")
rm(Eryth_small)

Eryth_med <- inv_biomass[inv_biomass$Sp_stage %in% c("Erythrodiplax_medium"), ]
Eryth_med$biom <- Eryth_med[,2]*Eryth_med[,4:73]
Eryth_med <- Eryth_med[,c(-1:-3)]
Eryth_med <- apply(Eryth_med, 2, sum)
Eryth_med <- as.data.frame(Eryth_med)
Eryth_med <- Eryth_med[c(71:140),]
Eryth_med <- as.data.frame(Eryth_med)
Eryth_med <- cbind(ID_plant, Eryth_med)
test <- merge(test, Eryth_med, by = "ID_plant")
rm(Eryth_med)

Eryth_large <- inv_biomass[inv_biomass$Sp_stage %in% c("Erythrodiplax_large"), ]
Eryth_large$biom <- Eryth_large[,2]*Eryth_large[,4:73]
Eryth_large <- Eryth_large[,c(-1:-3)]
Eryth_large <- apply(Eryth_large, 2, sum)
Eryth_large <- as.data.frame(Eryth_large)
Eryth_large <- Eryth_large[c(71:140),]
Eryth_large <- as.data.frame(Eryth_large)
Eryth_large <- cbind(ID_plant, Eryth_large)
test <- merge(test, Eryth_large, by = "ID_plant")
rm(Eryth_large)

Lept1_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Leptagrion.sp1_small", "Leptagrion.sp1_medium", "Leptagrion.sp1_large"), ]
Lept1_biom$biom <- Lept1_biom[,2]*Lept1_biom[,4:73]
Lept1_biom <- Lept1_biom[,c(-1:-3)]
Lept1_biom <- apply(Lept1_biom, 2, sum)
Lept1_biom <- as.data.frame(Lept1_biom)
Lept1_biom <- Lept1_biom[c(71:140),]
Lept1_biom <- as.data.frame(Lept1_biom)
Lept1_biom <- cbind(ID_plant, Lept1_biom)
test <- merge(test, Lept1_biom, by = "ID_plant")
rm(Lept1_biom)

Lept1_small <- inv_biomass[inv_biomass$Sp_stage %in% c("Leptagrion.sp1_small"), ]
Lept1_small$biom <- Lept1_small[,2]*Lept1_small[,4:73]
Lept1_small <- Lept1_small[,c(-1:-3)]
Lept1_small <- apply(Lept1_small, 2, sum)
Lept1_small <- as.data.frame(Lept1_small)
Lept1_small <- Lept1_small[c(71:140),]
Lept1_small <- as.data.frame(Lept1_small)
Lept1_small <- cbind(ID_plant, Lept1_small)
test <- merge(test, Lept1_small, by = "ID_plant")
rm(Lept1_small)

Lept1_med <- inv_biomass[inv_biomass$Sp_stage %in% c("Leptagrion.sp1_medium"), ]
Lept1_med$biom <- Lept1_med[,2]*Lept1_med[,4:73]
Lept1_med <- Lept1_med[,c(-1:-3)]
Lept1_med <- apply(Lept1_med, 2, sum)
Lept1_med <- as.data.frame(Lept1_med)
Lept1_med <- Lept1_med[c(71:140),]
Lept1_med <- as.data.frame(Lept1_med)
Lept1_med <- cbind(ID_plant, Lept1_med)
test <- merge(test, Lept1_med, by = "ID_plant")
rm(Lept1_med)

Lept1_large <- inv_biomass[inv_biomass$Sp_stage %in% c("Leptagrion.sp1_large"), ]
Lept1_large$biom <- Lept1_large[,2]*Lept1_large[,4:73]
Lept1_large <- Lept1_large[,c(-1:-3)]
Lept1_large <- apply(Lept1_large, 2, sum)
Lept1_large <- as.data.frame(Lept1_large)
Lept1_large <- Lept1_large[c(71:140),]
Lept1_large <- as.data.frame(Lept1_large)
Lept1_large <- cbind(ID_plant, Lept1_large)
test <- merge(test, Lept1_large, by = "ID_plant")
rm(Lept1_large)

Lept2_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Leptagrion.sp2_large"), ]
Lept2_biom$biom <- Lept2_biom[,2]*Lept2_biom[,4:73]
Lept2_biom <- Lept2_biom[,c(-1:-3)]
Lept2_biom <- apply(Lept2_biom, 2, sum)
Lept2_biom <- as.data.frame(Lept2_biom)
Lept2_biom <- Lept2_biom[c(71:140),]
Lept2_biom <- as.data.frame(Lept2_biom)
Lept2_biom <- cbind(ID_plant, Lept2_biom)
test <- merge(test, Lept2_biom, by = "ID_plant")
rm(Lept2_biom)

Elpi_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Elpidium.bromeliarum"), ]
Elpi_biom$biom <- Elpi_biom[,2]*Elpi_biom[,4:73]
Elpi_biom <- Elpi_biom[,c(-1:-3)]
Elpi_biom <- apply(Elpi_biom, 2, sum)
Elpi_biom <- as.data.frame(Elpi_biom)
Elpi_biom <- Elpi_biom[c(71:140),]
Elpi_biom <- as.data.frame(Elpi_biom)
Elpi_biom <- cbind(ID_plant, Elpi_biom)
test <- merge(test, Elpi_biom, by = "ID_plant")
rm(Elpi_biom)

Aulo_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Aulophorus.superterrenus"), ]
Aulo_biom$biom <- Aulo_biom[,2]*Aulo_biom[,4:73]
Aulo_biom <- Aulo_biom[,c(-1:-3)]
Aulo_biom <- apply(Aulo_biom, 2, sum)
Aulo_biom <- as.data.frame(Aulo_biom)
Aulo_biom <- Aulo_biom[c(71:140),]
Aulo_biom <- as.data.frame(Aulo_biom)
Aulo_biom <- cbind(ID_plant, Aulo_biom)
test <- merge(test, Aulo_biom, by = "ID_plant")
rm(Aulo_biom)

Enchy_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Enchytraeidae"), ]
Enchy_biom$biom <- Enchy_biom[,2]*Enchy_biom[,4:73]
Enchy_biom <- Enchy_biom[,c(-1:-3)]
Enchy_biom <- apply(Enchy_biom, 2, sum)
Enchy_biom <- as.data.frame(Enchy_biom)
Enchy_biom <- Enchy_biom[c(71:140),]
Enchy_biom <- as.data.frame(Enchy_biom)
Enchy_biom <- cbind(ID_plant, Enchy_biom)
test <- merge(test, Enchy_biom, by = "ID_plant")
rm(Enchy_biom)

micro_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Microculex.stonei_st1", "Microculex.stonei_st2", "Microculex.stonei_st3", "Microculex.stonei_st4", "Microculex.stonei_st5"), ]
micro_biom$biom <- micro_biom[,2]*micro_biom[,4:73]
micro_biom <- micro_biom[,c(-1:-3)]
micro_biom <- apply(micro_biom, 2, sum)
micro_biom <- as.data.frame(micro_biom)
micro_biom <- micro_biom[c(71:140),]
micro_biom <- as.data.frame(micro_biom)
micro_biom <- cbind(ID_plant, micro_biom)
test <- merge(test, micro_biom, by = "ID_plant")
rm(micro_biom)

wyeo_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Wyeomyia.sp_st2", "Wyeomyia.sp_st3", "Wyeomyia.sp_st4"), ]
wyeo_biom$biom <- wyeo_biom[,2]*wyeo_biom[,4:73]
wyeo_biom <- wyeo_biom[,c(-1:-3)]
wyeo_biom <- apply(wyeo_biom, 2, sum)
wyeo_biom <- as.data.frame(wyeo_biom)
wyeo_biom <- wyeo_biom[c(71:140),]
wyeo_biom <- as.data.frame(wyeo_biom)
wyeo_biom <- cbind(ID_plant, wyeo_biom)
test <- merge(test, wyeo_biom, by = "ID_plant")
rm(wyeo_biom)

tany_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Tanytarsus.sp_st2", "Tanytarsus.sp_st3", "Tanytarsus.sp_st4"), ]
tany_biom$biom <- tany_biom[,2]*tany_biom[,4:73]
tany_biom <- tany_biom[,c(-1:-3)]
tany_biom <- apply(tany_biom, 2, sum)
tany_biom <- as.data.frame(tany_biom)
tany_biom <- tany_biom[c(71:140),]
tany_biom <- as.data.frame(tany_biom)
tany_biom <- cbind(ID_plant, tany_biom)
test <- merge(test, tany_biom, by = "ID_plant")
rm(tany_biom)

chiro_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Chironomini"), ]
chiro_biom$biom <- chiro_biom[,2]*chiro_biom[,4:73]
chiro_biom <- chiro_biom[,c(-1:-3)]
chiro_biom <- apply(chiro_biom, 2, sum)
chiro_biom <- as.data.frame(chiro_biom)
chiro_biom <- chiro_biom[c(71:140),]
chiro_biom <- as.data.frame(chiro_biom)
chiro_biom <- cbind(ID_plant, chiro_biom)
test <- merge(test, chiro_biom, by = "ID_plant")
rm(chiro_biom)

trent_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Trentepohlia.sp1_st1","Trentepohlia.sp1_st2", "Trentepohlia.sp1_st3", "Trentepohlia.sp1_st4", "Trentepohlia.sp1_st4","Trentepohlia.sp1_st5" ), ]
trent_biom$biom <- trent_biom[,2]*trent_biom[,4:73]
trent_biom <- trent_biom[,c(-1:-3)]
trent_biom <- apply(trent_biom, 2, sum)
trent_biom <- as.data.frame(trent_biom)
trent_biom <- trent_biom[c(71:140),]
trent_biom <- as.data.frame(trent_biom)
trent_biom <- cbind(ID_plant, trent_biom)
test <- merge(test, trent_biom, by = "ID_plant")
rm(trent_biom)

scirtes_biom <- inv_biomass[inv_biomass$Sp_stage %in% c("Scirtes.sp_st1","Scirtes.sp_st2", "Scirtes.sp_st3", "Scirtes.sp_st4", "Scirtes.sp_st5","Scirtes.sp_>st5" ), ]
scirtes_biom$biom <- scirtes_biom[,2]*scirtes_biom[,4:73]
scirtes_biom <- scirtes_biom[,c(-1:-3)]
scirtes_biom <- apply(scirtes_biom, 2, sum)
scirtes_biom <- as.data.frame(scirtes_biom)
scirtes_biom <- scirtes_biom[c(71:140),]
scirtes_biom <- as.data.frame(scirtes_biom)
scirtes_biom <- cbind(ID_plant, scirtes_biom)
test <- merge(test, scirtes_biom, by = "ID_plant")
rm(scirtes_biom)
```

Predator species per treatment
```{r}
test2 <- dplyr::select(test, ID_plant, Effects_f, coreth_biom, Tanyp_biom, Eryth_biom, Lept1_biom, Lept2_biom)
test2$Effects_f <- factor(test2$Effects_f, levels=c('Open drought', 'Open control', 'Deforestation_Drought', 'Deforestation + Drought', 'Forest drought','Deforestation','Forest control'))

coreth <- dplyr::select(test2, ID_plant, Effects_f, coreth_biom)
coreth$group <- "Corethrella sp."
colnames(coreth)[3] <- "Biomass"

Tanyp <- dplyr::select(test2, ID_plant, Effects_f, Tanyp_biom)
Tanyp$group <- "Tanypodinae sp."
colnames(Tanyp)[3] <- "Biomass"

Eryth <- dplyr::select(test2, ID_plant, Effects_f, Eryth_biom)
Eryth$group <- "Erythrodiplax sp."
colnames(Eryth)[3] <- "Biomass"

Lept1 <- dplyr::select(test2, ID_plant, Effects_f, Lept1_biom)
Lept1$group <- "Leptagrion sp1"
colnames(Lept1)[3] <- "Biomass"

Lept2 <-dplyr::select(test2, ID_plant, Effects_f, Lept2_biom)
Lept2$group <- "Leptagrion sp2"
colnames(Lept2)[3] <- "Biomass"

test3 <- rbind(coreth, Tanyp, Eryth, Lept1, Lept2)
test3$group <- as.factor(test3$group)

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
 return(data_sum)
}
df2 <- data_summary(test3, varname="Biomass", 
                    groupnames=c("group", "Effects_f"))

df2$group_f=factor(df2$group, levels = c("Corethrella sp.", "Tanypodinae sp.", "Erythrodiplax sp.", "Leptagrion sp1", "Leptagrion sp2"))


p <- ggplot(df2, aes(x=Effects_f, y=Biomass, fill=group_f)) +  labs(fill = "") + ylab("Biomass (DW/mg)") +
  geom_bar(stat="identity", color = "black", position=position_dodge(),  alpha=0.8) + 
  geom_errorbar(aes(ymin=Biomass, ymax= Biomass+sd), width=.2,
                 position=position_dodge(.9)) + 
  theme(axis.title.y = element_blank(), 
        axis.title.x = element_text(face = "bold"), 
        axis.text.y = element_text(face = "bold", size = 10),
        axis.text.x = element_text(face = "bold", size = 11),
        legend.title = element_text(face = "bold"),
        legend.text = element_text(face = "italic"), 
        legend.position = "right")  + 
  labs(fill = "Predator species", face = "bold") + 
  coord_flip() 

ggsave("FigureS1.png", p, width = 20, height = 20, units = "cm")


rm(df2, data_summary, coreth, Tanyp, Eryth, Lept1, Lept2, p)
```

Prey species per treatment
```{r}
test2 <- dplyr::select(test, ID_plant, Effects_f, trent_biom, scirtes_biom, chiro_biom, tany_biom, Enchy_biom, Aulo_biom, Elpi_biom, micro_biom, wyeo_biom)
test2$Effects_f <- factor(test2$Effects_f, levels=c('Open drought', 'Open control', 'Deforestation_Drought', 'Deforestation + Drought', 'Forest drought','Deforestation','Forest control'))


trent <- dplyr::select(test2, ID_plant, Effects_f, trent_biom)
trent$group <- "Trentepohlia sp."
colnames(trent)[3] <- "Biomass"

scirt <- dplyr::select(test2, ID_plant, Effects_f, scirtes_biom)
scirt$group <- "Scirtes sp."
colnames(scirt)[3] <- "Biomass"

chiro <- dplyr::select(test2, ID_plant, Effects_f, chiro_biom)
chiro$group <- "Chironomini"
colnames(chiro)[3] <- "Biomass"

tany <- dplyr::select(test2, ID_plant, Effects_f, tany_biom)
tany$group <- "Tanytarsus sp."
colnames(tany)[3] <- "Biomass"

enchy <- dplyr::select(test2, ID_plant, Effects_f, Enchy_biom)
enchy$group <- "Enchytraeidae"
colnames(enchy)[3] <- "Biomass"

aulo <- dplyr::select(test2, ID_plant, Effects_f, Aulo_biom)
aulo$group <- "Aulophorus superterrenus"
colnames(aulo)[3] <- "Biomass"

elpi <- dplyr::select(test2, ID_plant, Effects_f, Elpi_biom)
elpi$group <- "Elpidium bromeliarum"
colnames(elpi)[3] <- "Biomass"

micro <- dplyr::select(test2, ID_plant, Effects_f, micro_biom)
micro$group <- "Microculex stonei"
colnames(micro)[3] <- "Biomass"

wyeo <- dplyr::select(test2, ID_plant, Effects_f, wyeo_biom)
wyeo$group <- "Wyeomyia sp."
colnames(wyeo)[3] <- "Biomass"


test3 <- rbind(trent, scirt, chiro, tany, enchy, aulo, elpi, micro, wyeo)
test3$group <- as.factor(test3$group)

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
 return(data_sum)
}

df2 <- data_summary(test3, varname="Biomass", 
                    groupnames=c("group", "Effects_f"))

df2$group_f=factor(df2$group, levels = c("Trentepohlia sp.", "Scirtes sp.", "Chironomini", "Tanytarsus sp.", "Enchytraeidae", "Aulophorus superterrenus", "Elpidium bromeliarum","Microculex stonei", "Wyeomyia sp."))


p <- ggplot(df2, aes(x=Effects_f, y=Biomass, fill=group_f)) +  labs(fill = "") + ylab("Biomass (DW/mg)") +
  geom_bar(stat="identity", color = "black", position=position_dodge(),  alpha=0.8) + 
  geom_errorbar(aes(ymin=Biomass, ymax= Biomass+sd), width=.2,
                 position=position_dodge(.9)) + 
  theme(axis.title.y = element_blank(), 
        axis.title.x = element_text(face = "bold"), 
        axis.text.y = element_text(face = "bold", size = 10),
        axis.text.x = element_text(face = "bold", size = 11),
        legend.title = element_text(face = "bold"),
        legend.text = element_text(face = "italic"), 
        legend.position = "right")  + 
  labs(fill = "Detritivore species", face = "bold") + 
  coord_flip() 

ggsave("FigureS2.png", p, width = 20, height = 20, units = "cm")

rm(df2, data_summary, trent, scirt, chiro, tany, enchy, aulo, elpi, micro, wyeo, p, test, test2, test3)
```
#VIII) Figure S3
```{r}
p <- ggplot(data_mesh, aes(x=Effects_f, y=Ergosterol, col =Treatment)) +
  geom_boxplot(outlier.alpha = 0, alpha=0.5) +
      geom_jitter(width=0.25)+
  facet_grid(rows = vars(Mesh.size)) + 
    labs(fill = "")  + 
 theme(axis.title.y = element_text(vjust = 2, size = 12, face = "bold"),
       axis.title.x = element_blank(),
       axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1, size = 10, face = "bold"),
       legend.position = "none") + 
  ylab("Biomass of fungi (μg/gDW)")

ggsave("FigureS3.png", p, width = 20, height = 20, units = "cm")
```