---
title: "Microplastics_Bromus_hordeaceus"
author: "Alonso-Crespo, I.M & Mateos-Cárdenas, A."
date: "2024-01-17"
output: html_document
---

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

# Set working directory

```{r}
setwd("C:/Users/Crespo/Desktop/Inés/Microplastic and plants/M&P/R")
# Put here your wd
#setwd("")

rm(list=ls())
```


## Librarys

```{r setup, include=FALSE}
library(readxl)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(tidyr)
library(viridis)
library(multcomp)
library(ggpubr)
library(DHARMa)
library(emmeans)
library(lme4)
```

## Load data: In Vitro experiment - Phase I: Germination

```{r}
#### In Vitro experiment: Germination test in Petri dishes of spiked seeds data ####

# Load and prepare data for germination test of spiked seeds
# This dataframe will be use to analyse the effect of plastic in sprout and radicle length development
# in Bromus hordeaceus seeds in Petri dishes

GT <- read_excel("In_Vitro_exp_germination.xlsx", na="NA")

GT$Plate<-as.factor(GT$Plate)
is.factor(GT$Plate)
levels(GT$Plate)

GT$Treatment<-as.factor(GT$Treatment)
is.factor(GT$Treatment)
levels(GT$Treatment)

GT$Replicate<-as.factor(GT$Replicate)
is.factor(GT$Replicate)
levels(GT$Replicate)

# Data frame for number of seeds successfully germinated per treatment
# This data frame will be used to test the effect of plastic presence on germination success
# of seeds growing in Petri dishes 

Data_Germination <- GT[, c(1,2,3,4)]
Data_Germination <- Data_Germination %>% 
  pivot_wider(names_from = "Replicate", 
               values_from = "Germination")
Data_Germination$GerminationSuccess <- rowSums(Data_Germination[, c(3:17)], na.rm = TRUE)


```

## Load data: In Vitro experiment - Phase II: Establishment

```{r, include=FALSE}
#### In Vitro experiment: Establishment test in soil of plastic-spiked germinated seeds data #### 

# Load and prepare data frame for establishment success of spiked seeds.
# This data frame will be used to plot the evolution of the establishment over time.

ES <- read_excel("In_Vitro_exp_establishment.xlsx", na="NA")

# Create one data frame with the information for establishment, second and third leave appearance.

EL1 <- ES[,c(1,2,3)]
EL2 <- ES[,c(1,2,4)]
EL3 <- ES[,c(1,2,5)]
colnames(EL1) <- c("Treatment", "Date", "N_Pots")
colnames(EL2) <- c("Treatment", "Date", "N_Pots")
colnames(EL3) <- c("Treatment", "Date", "N_Pots")
EL1$type <- "First_leaf"
EL2$type <- "Second_leaf"
EL3$type <- "Third_leaf"

# Join both dataframes and eliminate usefulness data
Data_Establishment <- rbind(EL1,EL2,EL3) # Total number of pots established and with a 3rd leaf.
rm(EL1,EL2,EL3)

Data_Establishment$Treatment<-as.factor(Data_Establishment$Treatment)
is.factor(Data_Establishment$Treatment)
levels(Data_Establishment$Treatment)

# Data frame to analyse the emergence speed of germinated seeds spiked with plastic particles

Establishment_success <- read_excel("In_Vitro_exp_establishment_dates.xlsx", na="NA")

Establishment_success$Treatment<-as.factor(Establishment_success$Treatment)
is.factor(Establishment_success$Treatment)
levels(Establishment_success$Treatment)

Establishment_success$Pot_ID<-as.factor(Establishment_success$Pot_ID)
is.factor(Establishment_success$Pot_ID)
levels(Establishment_success$Pot_ID)

Establishment_success$Replicate<-as.factor(Establishment_success$Replicate)
is.factor(Establishment_success$Replicate)
levels(Establishment_success$Replicate)

Establishment_success$Day_First_leaf<-as.numeric(Establishment_success$Day_First_leaf)
is.numeric(Establishment_success$Day_First_leaf)

Establishment_success$Day_Second_leaf<-as.numeric(Establishment_success$Day_Second_leaf)
is.numeric(Establishment_success$Day_Second_leaf)

Establishment_success$Day_Third_leaf<-as.numeric(Establishment_success$Day_Third_leaf)
is.numeric(Establishment_success$Day_Third_leaf)

Establishment_success$EstablishmentSuccess<-as.numeric(Establishment_success$EstablishmentSuccess)
is.numeric(Establishment_success$EstablishmentSuccess)

```

## Load data: In Soil experiment: Pot assay

```{r}
#### In Soil experiment: Pot assay #### 

# Load and prepare data for pot experiment: development of roots and shoots in 
# the presence of different plastic types and concentrations in soil.
# This dataframe will be use to test effect of the presence of plastics in soil in 
# plant development

Data_Biomass <- read_excel("In_Soil_exp.xlsx", na="NA")
Data_Biomass <- Data_Biomass %>%
  separate(col=Barcode, 
           into=c("PotID", 
                  "Particle", 
                  "Plastic",
                  "Replicate",  
                  "Concentration"),
           sep="_") %>%
  add_column(PDW_g=NA, .after="RDW_g")  %>%
  add_column(RMF_g=NA, .after="PDW_g")  

Data_Biomass$PDW_g<-(Data_Biomass$RDW_g+Data_Biomass$SDW_g)
Data_Biomass$RMF_g<-(Data_Biomass$RDW_g/Data_Biomass$PDW_g)

Data_Biomass$PotID<-as.factor(Data_Biomass$PotID)
is.factor(Data_Biomass$PotID)
levels(Data_Biomass$PotID)

Data_Biomass$Plastic<-as.factor(Data_Biomass$Plastic)
is.factor(Data_Biomass$Plastic)
levels(Data_Biomass$Plastic)

Data_Biomass$Particle<-as.factor(Data_Biomass$Particle)
is.factor(Data_Biomass$Particle)
levels(Data_Biomass$Particle)

Data_Biomass$Replicate<-as.factor(Data_Biomass$Replicate)
is.factor(Data_Biomass$Replicate)
levels(Data_Biomass$Replicate)

Data_Biomass$Concentration<-as.factor(Data_Biomass$Concentration)
is.factor(Data_Biomass$Concentration)
levels(Data_Biomass$Concentration)

#Data exploration (detect possible outliers)

ggplot(Data_Biomass, aes(x=SDW_g, y=PotID,  color=Concentration))+
  geom_point()+
  theme_bw()+
  facet_wrap(~Particle, labeller = label_wrap_gen(width=10))+
  xlab("Shoot dry weight (g)")+
  theme(legend.position = "none")

# Filter by Particle == microPE
microPE_data <- Data_Biomass %>% filter(Particle == "microPE")

# Calculate descriptive stats
mean_root <- mean(microPE_data$SDW_g, na.rm = TRUE)
sd_root <- sd(microPE_data$SDW_g, na.rm = TRUE)
iqr_root <- IQR(microPE_data$SDW_g, na.rm = TRUE)

# Calculate IQR (Interquartil range)
q1 <- quantile(microPE_data$SDW_g, 0.25, na.rm = TRUE)
q3 <- quantile(microPE_data$SDW_g, 0.75, na.rm = TRUE)
lower_bound <- q1 - 1.5 * iqr_root
upper_bound <- q3 + 1.5 * iqr_root

# Filter pot 7
pot_7_value <- Data_Biomass %>% filter(Particle == "microPE" & PotID == 7) %>% pull(SDW_g)

# Verify if pot 7 is on outlier
is_outlier <- pot_7_value < lower_bound | pot_7_value > upper_bound
is_outlier # TRUE %

ggplot(Data_Biomass, aes(x=RDW_g, y=PotID, color=Concentration))+
  geom_point()+
  theme_bw()+
  facet_wrap(~Particle, labeller = label_wrap_gen(width=10))+
  xlab("Root dry weight (g)")+
  theme(legend.position = "none")

# Filter by Particle == microPE
microPE_data <- Data_Biomass %>% filter(Particle == "microPE")

# Calculate descriptive stats
mean_root <- mean(microPE_data$RDW_g, na.rm = TRUE)
sd_root <- sd(microPE_data$RDW_g, na.rm = TRUE)
iqr_root <- IQR(microPE_data$RDW_g, na.rm = TRUE)

# Calculate IQR (Interquartil range)
q1 <- quantile(microPE_data$RDW_g, 0.25, na.rm = TRUE)
q3 <- quantile(microPE_data$RDW_g, 0.75, na.rm = TRUE)
lower_bound <- q1 - 1.5 * iqr_root
upper_bound <- q3 + 1.5 * iqr_root

# Filter pot 7
pot_7_value <- Data_Biomass %>% filter(Particle == "microPE" & PotID == 7) %>% pull(RDW_g)

# Verify if pot 7 is on outlier
is_outlier <- pot_7_value < lower_bound | pot_7_value > upper_bound
is_outlier # TRUE %

#Remove pot 7
Data_Biomass<-Data_Biomass[-c(7),]
```

## In Soil experiment: Nitrogen uptake

```{r}
#### N15 uptake data #### 

# Load and prepare data for main experiment: N15 uptake.

Data_N15 <- read_excel("In_Soil_exp_N15.xlsx", na="NA")

#Remove pot number 7

Data_N15<-Data_N15[-c(7),]

#Remove pot number 16. Nitrogen values were outside the detection limits of the elemental analyzer

Data_N15<-Data_N15[-c(15),]

Data_N15$PotID<-as.factor(Data_N15$PotID)
is.factor(Data_N15$PotID)
levels(Data_N15$PotID)

Data_N15$Replicate<-as.factor(Data_N15$Replicate)
is.factor(Data_N15$Replicate)
levels(Data_N15$Replicate)

Data_N15$Particle<-as.factor(Data_N15$Particle)
is.factor(Data_N15$Particle)
levels(Data_N15$Particle)

Data_N15$PlasticType<-as.factor(Data_N15$PlasticType)
is.factor(Data_N15$PlasticType)
levels(Data_N15$PlasticType)

Data_N15$Treatment<-as.factor(Data_N15$Treatment)
is.factor(Data_N15$Treatment)
levels(Data_N15$Treatment)
```

## In Vitro experiment - Phase I: Germination test in Petri dishes

```{r, include=FALSE}
###Germination success 

# Fit General Linear Model for germination success (Discrete variable)

glm <- glm(GerminationSuccess ~ Treatment, family = gaussian, data=Data_Germination)
summary(glm)
anova(glm, test = "Chisq") # NOT SIGNIFICANT

# Check model residuals

simulation <- simulateResiduals(fittedModel = glm, n = 250)
plot(simulation)

# Plot data germination success of seeds spiked with plastic in Petri dishes

FigureS4<-ggplot(data = Data_Germination, na.rm=TRUE, aes(x = Treatment, y= GerminationSuccess, group = Treatment, color= Treatment)) +
  geom_point(alpha= 0.2,na.rm=TRUE, position = position_jitter(width = 0.05, height = 0.05))+
  labs(title="",
       x = "Treatment",
       y = "Germinated seeds") +
theme_bw()+
stat_summary(fun.data = "mean_cl_boot", size=0.5)+
scale_color_viridis(discrete=T,
                    name="Treatment",
                    direction=-1,
                    option="C")+
scale_y_continuous(limit = c(0,15), breaks=seq(0,15, by=1))+
theme(legend.position = "right",
axis.text=element_text(size=10,color="black"),
axis.title.y = element_text(margin=margin(t=0,b=0,r=10,l=0)),
axis.text.x=element_text(angle=0))

FigureS4

#ggsave("FigureS4.jpg",
#       dpi=1200,
#       width=16,
#       height=9,
#       units="cm")
```

## In Vitro experiment - Phase I: Sprout length in Germinated seeds

```{r}
# Fit Linear Model for sprout length in germinated seeds
lm2<-lm(Sprout_Length~Treatment,data=GT[GT$Sprout_Length>0,],na.action = na.omit)
summary(lm2)
anova(lm2)

residlm3<-resid(lm2)
plot(residlm2) 
abline(0,0) 
qqnorm(residlm2) 
qqline(residlm2) 
plot(density(residlm2)) 
hist(residlm2, breaks=20)

# PostHoc

summary(glht(lm2,linfct=mcp(Treatment="Tukey"))) # microPBAT =! Control nanoPE =! Control microPE =! Control microPBAT =! microPE

table_lm2<-data.frame(Treatment=rep(c("Control","microPBAT" ,"microPE", "nanoPE"),5),
                       Label=rep(c("a", "b", "c", "bc"), 1),
                       y=c(rep(3, 5))) 

# Plot data average sprout development in seeds spiked with plastic in Petri dishes

TotalSprout <-ggplot(data = GT[GT$Sprout_Length>0,], na.rm=TRUE, aes(x = Treatment, y= Sprout_Length, group = Treatment, color= Treatment)) +
  geom_point(alpha=0.1,na.rm=TRUE, position = position_jitter(width = 0.05, height = 0.05))+
  labs(title="",
       x = "Treatment",
       y = "Sprout length in germinated seeds (cm)") +
theme_bw()+
stat_summary(fun.data = "mean_cl_boot", size=0.5)+
scale_color_viridis(discrete=T,
                    name="Treatment",
                    direction=-1,
                    option="C")+
scale_y_continuous(limit = c(0,3), breaks=seq(0,3, by=0.5))+
geom_text(data=table_lm2, mapping=aes(x=Treatment, y=y, label=Label), inherit.aes = FALSE)+
theme(legend.position = "right",
axis.text=element_text(size=10,color="black"),
axis.title.y = element_text(margin=margin(t=0,b=0,r=10,l=0)),
axis.text.x=element_text(angle=45,
                         hjust=1,
                         vjust=1))
TotalSprout


#ggsave("C:/Users/Crespo/Desktop/Inés/Microplastic and plants/Microplastic and Plants/GT_Total_Root_Development.jpg",
#       dpi=1200,
#       width=16,
#       height=9,
#       units="cm")
```

## In Vitro experiment - Phase I: Radicle lenght development in Germinated seeds

```{r}
# Fit Linear Model for radicle development (Summary of the length of the radicles produced by seed)
lm3<-lm(Radicle_Length~Treatment,data=GT[GT$Radicle_Length>0,],na.action = na.omit)
summary(lm3)
anova(lm3)

residlm3<-resid(lm3)
plot(residlm3) 
abline(0,0) 
qqnorm(residlm3) 
qqline(residlm3) 
plot(density(residlm3)) 
hist(residlm2, breaks=10)

# PostHoc

summary(glht(lm3,linfct=mcp(Treatment="Tukey"))) # microPBAT =! Control microPE =! Control nanoPE =! Control microPBAT =! microPE
table_lm3<-data.frame(Treatment=rep(c("Control","microPBAT" ,"microPE", "nanoPE"),5),
                       Label=rep(c("a", "b", "c", "cb"), 1),
                       y=c(rep(4, 5)))

# Plot radicle length develpment in seeds spiked with plastic in Petri dishes

TotalRadicle <-ggplot(data = GT[GT$Radicle_Length>0,], na.rm=TRUE, aes(x = Treatment, y= Radicle_Length, group = Treatment, color= Treatment)) +
  geom_point(alpha=0.1,na.rm=TRUE, position = position_jitter(width = 0.05, height = 0.05))+
  labs(title="",
       x = "Treatment",
       y = "Radicle length in germinated seeds (cm)") +
theme_bw()+
stat_summary(fun.data = "mean_cl_boot", size=0.5)+
scale_color_viridis(discrete=T,
                    name="Treatment",
                    direction=-1,
                    option="C")+
scale_y_continuous(limit = c(0,4), breaks=seq(0,4, by=0.5))+
geom_text(data=table_lm3, mapping=aes(x=Treatment, y=y, label=Label), inherit.aes = FALSE)+
theme(legend.position = "right",
axis.text=element_text(size=10,color="black"),
axis.title.y = element_text(margin=margin(t=0,b=0,r=10,l=0)),
axis.text.x=element_text(angle=45,
                         hjust=1,
                         vjust=1))
TotalRadicle



#ggsave("C:/Users/Crespo/Desktop/Inés/Microplastic and plants/Microplastic and Plants/GT_Total_Root_Development.jpg",
#       dpi=1200,
#       width=16,
#       height=9,
#       units="cm")
```

```{r}
pp1 <- get_legend(TotalRadicle)

Figure3 <- ggarrange(TotalSprout, TotalRadicle, ncol = 2, nrow = 1,  labels = c("a.","b."),
          legend	= "right",legend.grob=  pp1)

Figure3

#ggsave("Figure3.jpg",
#       dpi=1200,
#       width=25,
#       height=15,
#       units="cm")
```

#### In Vitro experiment - Phase II: Establishment in cleaned soil of spiked germinated seeds #### 

```{r}
# Data exploration
Figure4 <-ggplot() + 
  geom_point(data = Data_Establishment,
           mapping = aes(x=Date , y=N_Pots, color=Treatment, shape = type), na.rm=TRUE,  position = position_jitterdodge(dodge.width = 0.2, jitter.height = 0.1), size= 1.5, alpha = 0.3)+
  geom_smooth(data = Data_Establishment,
            mapping = aes(x =Date, y=N_Pots, color = Treatment, fill= Treatment, group = Treatment), alpha = 0.1) +
  xlab("Day") + ylab("Number of established plants") +
  theme_bw()+
    scale_fill_viridis(discrete=T,
                    name="Treatment",
                    direction=-1,
                    option="C")+
     scale_color_viridis(discrete=T,
                    name="Treatment",
                    direction=-1,
                    option="C")+
  scale_y_continuous(limit = c(-1,5.5), breaks=seq(0,5, by=1))

Figure4

#ggsave("Figure4.jpg",
#       dpi=1200,
#       width=16,
#       height=9,
#       units="cm")

# Fit Linear Model for the speed in the appearance of the first leaf

lm1leaf<-lm(First_leaf~ Treatment,data=Establishment_success,na.action = na.omit)
summary(lm1leaf)
anova(lm1leaf) # No significant


# Fit Linear Model for the speed in the appearance of the second leaf

lm2leaf<-lm(Second_leaf~ Treatment,data=Establishment_success,na.action = na.omit)
summary(lm2leaf)
anova(lm2leaf) # No significant

# Fit Linear Model for the speed in the appearance of the third leaf

lm3leaf<-lm(Third_leaf~ Treatment,data=Establishment_success,na.action = na.omit)
summary(lm3leaf)
anova(lm3leaf) # No significant

# Fit General linear model for establishment success.

Estab_succ <- data.frame(
  Treatment = c("Control", "microPE", "microPBAT", "nanoPE"),
  Pots_established = c(3, 4, 2, 4),
  Total_pots = c(5, 5, 5, 5))


model <- glm(cbind(Pots_established, Total_pots - Pots_established) ~ Treatment, data = Estab_succ, family = binomial)
summary(model) # No significant
```

# In Soil experiment - Pot assay: Plant development in soil spiked with plastic particules

## In Soil experiment: Shoot dry weight

```{r}
# Fit Linear Model for shoot dry weight

lm4<-lm(SDW_g~Particle*Concentration,data=Data_Biomass,na.action = na.omit)
summary(lm4)
anova(lm4)

residlm4<-resid(lm4)
plot(residlm4) 
abline(0,0) 
qqnorm(residlm4) 
qqline(residlm4) 
plot(density(residlm4)) 
hist(residlm4, breaks=5)

# PostHoc

postHoclm4 <- lsmeans(lm4, pairwise~Particle:Concentration)
summary(postHoclm4, type="response")

# Create a data frame for annotation
Treatment <- c("microPBAT", "microPE", "nanoPE")
Concentration <- c(rep(c("High", "Low"), each = 3))
y_valueslm4 <- c(1,2.5,2.5,2.2,2.5,2.5,2.5)  # Replace with your actual y-values
significant_letterslm4 <- c("b","a","ac","c","ac","ac","ac") 

annotation_datalm4 <- data.frame(Treatment, Concentration)
control <- c("Control", "Control")
annotation_datalm4 <- rbind(annotation_datalm4, control)
annotation_datalm4 <- cbind(annotation_datalm4,y_valueslm4,significant_letterslm4)

# Plot data shoot dry weight

SDW_ME <- ggplot(data = Data_Biomass, aes(x = Particle, y = SDW_g, color = Concentration)) +
  geom_jitter(position = position_dodge(width = 0.8), size = 0.1) +
  stat_summary(fun.data = "mean_cl_boot", 
               position = position_dodge(width = 0.8),
               size = 0.3) +
  geom_text(data = annotation_datalm4, 
            aes(x = Treatment, y = y_valueslm4, label = significant_letterslm4),
            position = position_dodge(width = 0.8),
            vjust = -0.5, hjust = 0.5, size = 4,
            show.legend = FALSE) +  
  labs(title = "",
       x = "",
       y = "Shoot dry weight (g)") +
  theme_bw() +
  scale_color_viridis(discrete = TRUE,
                      name = "Concentration",
                      direction = -1,
                      option = "D") +
  scale_y_continuous(limit = c(0, 4.5), breaks = seq(0, 4.5, by = 0.5)) +
  theme(legend.position = "right",
        axis.text = element_text(size = 12, color = "black"),
        axis.title.y = element_text(margin = margin(t = 0, b = 0, r = 10, l = 0)),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

SDW_ME
```
## In Soil experiment: Root dry weight

```{r}
# Fit Linear Model for root dry weight

lm5<-lm(RDW_g~Particle*Concentration,data=Data_Biomass,na.action = na.omit)
summary(lm5)
anova(lm5)

residlm5<-resid(lm5)
plot(residlm5) 
abline(0,0) 
qqnorm(residlm5) 
qqline(residlm5) 
plot(density(residlm5)) 
hist(residlm5, breaks=5)

# PostHoc

postHoclm5 <- lsmeans(lm5, pairwise~Particle:Concentration)
summary(postHoclm5, type="response")

# Create a data frame for annotation
Treatment <- c("microPBAT", "microPE", "nanoPE")
Concentration <- c(rep(c("High", "Low"), each = 3))
y_valueslm5 <- c(0.5,2.5,2.5,2,2.5,2.5,2)  # Replace with your actual y-values
significant_letterslm5 <- c("b","ac","a","c","a","a","ac") 

annotation_datalm5 <- data.frame(Treatment, Concentration)
control <- c("Control", "Control")
annotation_datalm5 <- rbind(annotation_datalm5, control)
annotation_datalm5 <- cbind(annotation_datalm5,y_valueslm5,significant_letterslm5)

# Plot data root dry weight

RDW_ME<-ggplot(data = Data_Biomass, na.rm=TRUE, aes(x = Particle, y= RDW_g, color= Concentration)) +
  geom_jitter(na.rm=TRUE, position=position_dodge(width=0.8), 
             size=0.1)+
  stat_summary(fun.data="mean_cl_boot", 
               position=position_dodge(width=0.8),
               size=0.3)+
  labs(title="",
       x = "",
       y = "Root dry weight (g)") +
theme_bw()+
scale_color_viridis(discrete=T,
                    name="Concentration",
                    direction=-1,
                    option="D")+
scale_y_continuous(limit = c(0,4.5), breaks=seq(0,4.5, by=0.5))+  
geom_text(data = annotation_datalm5, 
            aes(x = Treatment, y = y_valueslm5, label = significant_letterslm5),
            position = position_dodge(width = 0.8),
            vjust = -0.5, hjust = 0.5, size = 4,
            show.legend = FALSE) +  
theme(legend.position = "right",
axis.text=element_text(size=12,color="black"),
axis.title.y = element_text(margin=margin(t=0,b=0,r=10,l=0)),
axis.text.x=element_text(angle=45,
                         hjust=1,
                         vjust=1))


RDW_ME
```

## In Soil experiment: Plant dry weight

```{r}
# Fit Linear Model for plant dry weight

lm6<-lm(PDW_g~Particle*Concentration,data=Data_Biomass,na.action = na.omit)
summary(lm6)
anova(lm6)

residlm6<-resid(lm6)
plot(residlm6) 
abline(0,0) 
qqnorm(residlm6) 
qqline(residlm6) 
plot(density(residlm6)) 
hist(residlm6, breaks=5)

# PostHoc

postHoclm6 <- lsmeans(lm6, pairwise~Particle:Concentration)
summary(postHoclm6, type="response")

# Create a data frame for annotation
Treatment <- c("microPBAT", "microPE", "nanoPE")
Concentration <- c(rep(c("High", "Low"), each = 3))
y_valueslm6 <- c(1,4.2,4.2,3.2,4.2,4.2,4)  # Replace with your actual y-values
significant_letterslm6 <- c("b","ad","d","c","ad","ad","ac") 

annotation_datalm6 <- data.frame(Treatment, Concentration)
control <- c("Control", "Control")
annotation_datalm6 <- rbind(annotation_datalm6, control)
annotation_datalm6 <- cbind(annotation_datalm6,y_valueslm6,significant_letterslm6)

# Plot data Plant dry weight

PDW_ME<-ggplot(data = Data_Biomass, na.rm=TRUE, aes(x = Particle, y= PDW_g, color= Concentration)) +
  geom_jitter(na.rm=TRUE, position=position_dodge(width=0.8), 
             size=0.1)+
  stat_summary(fun.data="mean_cl_boot", 
               position=position_dodge(width=0.8),
               size=0.3)+
  labs(title="",
       x = "",
       y = "Plant dry weight (g)") +
theme_bw()+
scale_color_viridis(discrete=T,
                    name="Concentration",
                    direction=-1,
                    option="D")+
scale_y_continuous(limit = c(0,4.5), breaks=seq(0,4.5, by=0.5))+
geom_text(data = annotation_datalm6, 
            aes(x = Treatment, y = y_valueslm6, label = significant_letterslm6),
            position = position_dodge(width = 0.8),
            vjust = -0.5, hjust = 0.5, size = 4,
            show.legend = FALSE) + 
theme(legend.position = "right",
axis.text=element_text(size=12,color="black"),
axis.title.y = element_text(margin=margin(t=0,b=0,r=10,l=0)),
axis.text.x=element_text(angle=45,
                         hjust=1,
                         vjust=1))
PDW_ME
```

```{r}
pp2 <- get_legend(SDW_ME)

Figure5 <- ggarrange(SDW_ME, RDW_ME, PDW_ME, ncol = 3, nrow = 1,  labels = c("a.","b.","c."),
          legend	= "right",legend.grob=  pp2)

Figure5

ggsave("Figure5.jpg",
       dpi=1200,
       width=32,
       height=12,
       units="cm")
```

## Nitrogen 15

```{r}
# Average of delta N15

# Fit Linear Model for average delta N15

lm8<-lm(Avg_Delta_N15~Treatment*Concentration,data=Data_N15,na.action = na.omit)
summary(lm8)
anova(lm8) # INTERACTION NO SIGNIFICANT

lm8<-lm(Avg_Delta_N15~Treatment,data=Data_N15,na.action = na.omit)
anova(lm8) 

residlm8<-resid(lm8)
plot(residlm8) 
abline(0,0) 
qqnorm(residlm8) 
qqline(residlm8) 
plot(density(residlm8)) 
hist(residlm8, breaks=5)

# PostHoc

postHoclm8 <- lsmeans(lm8, pairwise~Treatment)
summary(postHoclm8, type="response")

# Annotation table

table_lm8<-data.frame(Treatment=rep(c("Control","microPBAT" ,"microPE", "nanoPE"),5),
                       Label=rep(c("a", "b", "a", "a"), 1),
                       y=c(rep(2000, 2000))) 

# Plot data Average delta N15

FigureS5<-ggplot(data = Data_N15, na.rm=TRUE, aes(x = Treatment, y= Avg_Delta_N15, color= Treatment)) +
  geom_jitter(alpha=0.2,na.rm=TRUE, position=position_dodge(width=0.8))+
  stat_summary(fun.data="mean_cl_boot", 
               position=position_dodge(width=0.8),
               size=0.5)+
  labs(title="",
       x = "",
       y = "Average delta (δ) ^15N (‰)") +
theme_bw()+
geom_text(data=table_lm8, mapping=aes(x=Treatment, y=y, label=Label), inherit.aes = FALSE)+
scale_color_viridis(discrete=T,
                    name="Treatment",
                    direction=-1,
                    option="C")+
theme(legend.position = "right",
axis.text=element_text(size=12,color="black"),
axis.title.y = element_text(margin=margin(t=0,b=0,r=10,l=0)),
axis.text.x=element_text(angle=45,
                         hjust=1,
                         vjust=1))

FigureS5

#ggsave("FigureS5.jpg",
#       dpi=1200,
#       width=12,
#       height=10,
#       units="cm")

# Average deta N15 x concentration

Figure6a<-ggplot(data = Data_N15, na.rm=TRUE, aes(x = Treatment, y= Avg_Delta_N15, color= Concentration)) +
  geom_point(alpha=0.2, na.rm=TRUE, position=position_dodge(width=0.8))+
  stat_summary(fun.data="mean_cl_boot", 
               position=position_dodge(width=0.8),
               size=0.5)+
  labs(title="",
       x = "",
       y = "Average delta (δ) ^15N (‰)") +
theme_bw()+
scale_color_viridis(discrete=T,
                    name="Concentration",
                    direction=-1,
                    option="D")+
theme(legend.position = "right",
axis.text=element_text(size=12,color="black"),
axis.title.y = element_text(margin=margin(t=0,b=0,r=10,l=0)),
axis.text.x=element_text(angle=45,
                         hjust=1,
                         vjust=1))

Figure6a


# Fit Linear Model for nitrogen accumulation

lm9<-lm(Naccum~Treatment*Concentration,data=Data_N15,na.action = na.omit)
summary(lm9)
anova(lm9)

residlm9<-resid(lm9)
plot(residlm9) 
abline(0,0) 
qqnorm(residlm9) 
qqline(residlm9) 
plot(density(residlm9)) 
hist(residlm9, breaks=5)

# PostHoc

postHoclm9 <- lsmeans(lm9, pairwise~Treatment*Concentration)
summary(postHoclm9, type="response")

# Create a data frame for annotation
Treatment <- c("microPBAT", "microPE", "nanoPE")
Concentration <- c(rep(c("high", "low"), each = 3))
y_valueslm9 <- c(0.027,0.027,0.027,0.027,0.027,0.027,0.027)  # Replace with your actual y-values
significant_letterslm9 <- c("d","b","ac","a","bc","ac","abc") 

annotation_datalm9 <- data.frame(Treatment, Concentration)
control <- c("Control", "Control")
annotation_datalm9 <- rbind(annotation_datalm9, control)
annotation_datalm9 <- cbind(annotation_datalm9,y_valueslm9,significant_letterslm9)

# Nitrogen accumulation effect

Figure6b<-ggplot(data = Data_N15, na.rm=TRUE, aes(x = Treatment, y= Naccum, color= Concentration)) +
  geom_point(alpha=0.2, na.rm=TRUE, position=position_dodge(width=0.8))+
  stat_summary(fun.data="mean_cl_boot", 
               position=position_dodge(width=0.8),
               size=0.5)+
  labs(title="",
       x = "",
       y = "Nitrogen accumulation (g of N x g of plant material)") +
theme_bw()+
scale_y_continuous(limit = c(0,0.030), breaks=seq(0,0.030, by=0.005))+
geom_text(data = annotation_datalm9, 
            aes(x = Treatment, y = y_valueslm9, label = significant_letterslm9),
            position = position_dodge(width = 0.8),
            vjust = -0.5, hjust = 0.5, size = 4,
            show.legend = FALSE) + 
scale_color_viridis(discrete=T,
                    name="Concentration",
                    direction=-1,
                    option="D")+
theme(legend.position = "right",
axis.text=element_text(size=12,color="black"),
axis.title.y = element_text(margin=margin(t=0,b=0,r=10,l=0)),
axis.text.x=element_text(angle=45,
                         hjust=1,
                         vjust=1))

Figure6b

```

```{r}

pp3 <- get_legend(Figure6b)

Figure6 <- ggarrange(Figure6a, Figure6b, ncol = 2, nrow = 1,  labels = c("a.","b."),
          legend	= "right",legend.grob=  pp3)

Figure6

ggsave("Figure6.jpg",
       dpi=1200,
       width=32,
       height=12,
       units="cm")
```