######## Shepherd et al., 2022 Peat moss inoculation depends on water table depth

rm(list=ls())
library(tidyverse); library(ggplot2); library(vegan); library(phyloseq); library(plotrix) 
library(ggpubr); library(FD); library(Matching); library(ggfortify)

# setwd("~/Documents/Dutch_innoculation_exp/Data") ### set to your working directory to data location
setwd("C:/Users/k2256287/OneDrive - King's College London/Documents/PhD_work/Inoc_study/Repository_data/PURE_round1")
#### Aboveground taxonomic and functional analysis

## vegetation surveys
veg <- read.csv("Shepherd_Inoc_ABG_surveys.csv") ## vegetation surveys 
sp <- read.csv("Shepherd_Inoc_Vegetation_des.csv") ## Vegetation descriptions

### remove non-plant species and mosses
sp %>% filter(Group == "Plant") -> plants
pl_list <- plants$abbrev
veg %>% dplyr::select(Inoculated:Group) -> vab_trts ### gets treatments for analysis

veg %>% dplyr::select(pl_list) %>% 
  dplyr::select(-Epi_ang, -Era_mul, -Mat_cha, -Bet_pen, -Pri_vul) -> pl_ab ### removed rare species for analysis
NMDS.pl=metaMDS(pl_ab, k=2) 
NMDS.pl.2<- metaMDS(pl_ab, k=2) #
pro <- procrustes(NMDS.pl, NMDS.pl.2); pro ### low: all good

adonis(dist(decostand(pl_ab, method = "hellinger"))~ Water_table*Inoculated, data=vab_trts)
adonis(dist(decostand(pl_ab, method = "hellinger"))~ Inoculated*Water_table, data=vab_trts)

### prep for plotting

scores <- scores(NMDS.pl)$sites; scores <- as.data.frame(scores)
colnames(scores) <- c("NMDS1", "NMDS2")
sp_scores <- NMDS.pl$species
as.data.frame(sp_scores) %>% rownames_to_column("Species") %>% dplyr::filter(MDS1 != "NaN") -> spsc_df

cbind(scores, vab_trts) %>%
  mutate_if(is.character, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") %>%
  rename()-> nmds_scores2

# plot fig 2a
fig2a <- ggplot() +
  geom_point(aes(x=NMDS1, y=NMDS2, colour=Water_table, shape=Inoculated), data=nmds_scores2, size=2) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_shape_manual(values=c(1, 19)) +
  guides(color=guide_legend("Water table")) +
  xlab("NMDS1") +
  ylab("NMDS2") +
  theme_bw() +
  ggtitle("(a) Vascular plant \n     community composition") +
  theme(axis.text = element_text(size=12), 
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

# Functional composition

read.csv("Shepherd_Inoc_ABG_trts.csv") %>%
  column_to_rownames("X") %>%
  rename("H" = "Plant_height") %>%
  rownames_to_column("Species") %>%
  dplyr::filter(Species %in% c("Mol_cae", "Jun_bul", "Siu_lat", "Bid_fer", "Epi_hir", "Eri_tet")) %>%
  column_to_rownames("Species") -> vptrts

pl_ab %>% dplyr::select(-Pop_del, -Cys_fra, -Lys_vul) -> pl_ab2
vegmat <- as.matrix(pl_ab2)

functcomp(vptrts, vegmat, CWM.type = "all") -> cwm_trts
cwm_pca <- prcomp(cwm_trts, center=TRUE, scale=TRUE)
summary(cwm_pca)

vab_trts %>%
  mutate_if(is.character, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") -> pca_trt

fig2b <- autoplot(cwm_pca, data=pca_trt, colour = "Water_table", shape = "Inoculated", loadings=TRUE, loadings.colour = "gray55",
                  loadings.label = TRUE, loadings.label.size = 4, label.size=5, 
                  loadings.label.vjust = 1.2, loadings.label.hjust = 0.4,
                  loadings.label.colour = "gray1", size=2) +  
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_shape_manual(values=c(1, 19)) +
  guides(color=guide_legend("Water table")) + 
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text=element_text(face="bold", size=10))+
  ggtitle("(b) Vascular plant \n     functional composition") +
  ylab("PC2 (26%)") +
  xlab("PC1 (70%)") +
  theme(axis.text = element_text(size=12), 
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) +
  xlim(-0.35, 0.35)

########## Bryophyte composition

sp %>% dplyr::filter(Group == "Sphagnum" | Group == "Moss") -> sp_bryo
bryo_list <- sp_bryo$abbrev

veg %>% dplyr::select(bryo_list) -> br_ab

br_ab2 <- br_ab %>% 
  rowwise() %>% 
  filter(sum(c(Pol_com, Bry_spp, Sph_fal, Sph_cus, Hyp_cup, Dic_sco)) != 0)

NMDS.br=metaMDS(br_ab2, k=2) 
NMDS.br2 <- metaMDS(br_ab2, k=2) #
pro <- procrustes(NMDS.br, NMDS.br2); pro ### low: all good

bry.scores <- scores(NMDS.br)$sites; bry.scores <- as.data.frame(bry.scores)

bry_sp_scores <- NMDS.br$species
as.data.frame(bry_sp_scores) %>% rownames_to_column("Species") %>% dplyr::filter(MDS1 != "NaN") -> bry_spsc_df

cbind(br_ab, vab_trts) ### get rid of YL4 and YL1 whhere there is no moss cover

vab_trts %>% unite(trt, c("Inoculated", "Water_table", "Group"), remove=FALSE) %>%
  dplyr::filter(trt != "Y_L_4" & trt != "Y_L_1") -> bry_trts

cbind(bry.scores, bry_trts) -> bry_scoresb

bry_scoresb %>%
  mutate_if(is.character, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") -> bry_scores2

### permanova analysis
adonis(dist(decostand(br_ab2, method = "hellinger"))~ Water_table*Inoculated, data=bry_scoresb)
adonis(dist(decostand(br_ab2, method = "hellinger"))~ Inoculated*Water_table, data=bry_scoresb)


fig2c <- ggplot() +
  geom_point(aes(x=NMDS1, y=NMDS2, colour=Water_table, shape=Inoculated), data=bry_scores2, size = 2) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_shape_manual(values=c(1, 19)) +
  guides(color=guide_legend("Water table")) +
  xlab("NMDS1") +
  ylab("NMDS2") +
  theme_bw()+
  ggtitle("(c) Bryophyte \n     community composition") +
  theme(axis.text = element_text(size=12), 
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

##### Sphagnum cover ####

sp %>% dplyr::filter(Group == "Sphagnum") -> sp_sphag
sphag_sp_nms <- sp_sphag$abbrev
veg %>% dplyr::select(sphag_sp_nms)-> sphag_ab
as.data.frame(rowSums(sphag_ab)) -> sphag_cover
colnames(sphag_cover) <- c("Sphagnum")

cbind(sphag_cover, vab_trts) %>% 
  unite(Treatment, c("Inoculated", "Water_table"), remove = FALSE) -> sphag_cover2

sphag_cover2 %>% dplyr::select(Sphagnum, Inoculated, Water_table) %>%
  mutate_if(is.character, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") %>%
  unite(Treatment, c("Inoculated", "Water_table"), remove=FALSE) -> sphag_cover3

sphag_cover3 %>% group_by(Inoculated, Water_table, Treatment) %>%
  dplyr::summarise(across(c("Sphagnum"), ~ mean(.x, na.rm = TRUE))) -> ms_means

sphag_cover3 %>% group_by(Inoculated, Water_table, Treatment) %>%
  dplyr::summarise(across(c("Sphagnum"), ~ std.error(.x, na.rm = TRUE))) %>%
  rename(ci_sphag = Sphagnum) %>% ungroup() %>%
  dplyr::select(Treatment, ci_sphag) %>% full_join(ms_means, by="Treatment") -> ms_values

ms_values$Water_table <- factor(ms_values$Water_table, levels=c("Low", "High"))
sphag_cover3$Water_table<- factor(sphag_cover3$Water_table, levels=c("Low", "High"))

### Using bootstrapped KS boot: just for inoculated samples as no establishment when not inoculated
sphag_cover3 %>% filter(Water_table == "High" & Inoculated == "Yes") -> sphag_inhigh
sphag_cover3 %>% filter(Water_table == "Low" & Inoculated == "Yes") -> sphag_inlow
ks.boot(Tr = sphag_inhigh$Sphagnum, Co = sphag_inlow$Sphagnum, nboot=1000)
## value is variable due to bootstrapping

fig2d <- ggplot(aes(x=Inoculated, y=Sphagnum, fill=Water_table, colour=Water_table), data=ms_values) +
  geom_pointrange(aes(ymax=Sphagnum+ci_sphag,
                      ymin=Sphagnum-ci_sphag), position = position_dodge(0.4)) +
  geom_jitter(aes(x=Inoculated, y=Sphagnum, fill=Water_table, shape=Inoculated), data=sphag_cover3,
              alpha=0.6, position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.4)) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  guides(color=guide_legend("Water table"), fill = "none") +
  scale_shape_manual(values=c(1, 19)) +
  ylab("Cover (%)") +
  theme_bw() +
  ggtitle(expression(paste("(d) ",italic("Sphagnum"),"\n cover"))) +
  theme(axis.text = element_text(size=12), 
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) 

#### Make final figure
fig2<- ggarrange(fig2a, fig2b, fig2c, fig2d, ncol=2, nrow=2, common.legend = TRUE, legend = "right",
                 align = "hv")
fig2
# export 800 x 560



########## Appendix figures: Individial species covers + community weighted means for each trait


sp %>% dplyr::filter(Group == "Moss") -> sp_moss
ms_sp_nms <- sp_moss$abbrev
veg %>%
  mutate_if(is.factor, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") %>%
  unite(Treatment, c("Inoculated", "Water_table"), remove=FALSE) %>%
  dplyr::select(Inoculated, Water_table, Treatment, pl_list, ms_sp_nms, sphag_sp_nms) -> veg2

veg2 %>% group_by(Inoculated, Water_table, Treatment) %>%
  dplyr::summarise(across(c("Mol_cae":"Sph_cus"), ~ mean(.x, na.rm = TRUE)))  %>%
  gather("species", "mean", 4:23) %>% unite(trt2, c("Treatment", "species")) %>% ungroup() %>%
  dplyr::select(trt2, mean) -> veg_means

sp %>% dplyr::select(Species_t2, abbrev) %>% rename(species = abbrev) -> sp_abbrev

veg2 %>% group_by(Inoculated, Water_table, Treatment) %>%
  dplyr::summarise(across(c("Mol_cae":"Sph_cus"), ~ std.error(.x, na.rm = TRUE)))  %>%
  gather("species", "se", 4:23) %>% unite(trt2, c("Treatment", "species"), remove=FALSE) %>% 
  full_join(veg_means, by="trt2") %>% dplyr::select(Inoculated, Water_table, species, mean, se) %>%
  dplyr::filter(species != "Epi_ang") %>% left_join(sp_abbrev, by="species") -> veg_vals

veg   %>%
  mutate_if(is.character, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") %>%
  dplyr::select(Inoculated, Water_table, Group, pl_list, ms_sp_nms, sphag_sp_nms) %>%
  gather("species", "value", 4:23) %>% dplyr::filter(species != "Epi_ang")  %>% 
  left_join(sp_abbrev, by="species") -> veg_inds

### Individual species level covers
ggplot(aes(x=Inoculated, y=mean, fill=Water_table, colour=Water_table), data= veg_vals) +
  geom_pointrange(aes(ymax=mean+se,
                      ymin=mean-se), position = position_dodge(0.4)) +
  geom_jitter(aes(x=Inoculated, y=value, fill=Water_table, colour=Water_table, shape=Inoculated), data=veg_inds,
              alpha=0.6, position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.4)) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  guides(color=guide_legend("Water table"), fill = "none") +
  scale_shape_manual(values=c(1, 19)) +
  ylab("Cover (%)") +
  theme_bw() +
  facet_wrap(~Species_t2, scales = "free_y", strip.position = "top") +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text=element_text(face="italic", size=12))

###### community weighted means

cbind(veg2, cwm_trts)  %>% dplyr::select(Inoculated:Water_table, EMV:Plant.height) %>%
  unite(Treatment, c("Inoculated", "Water_table"), remove=FALSE) -> cwms 

cwms %>% group_by(Inoculated, Water_table, Treatment) %>%
  dplyr::summarise(across(c("EMV":"Plant.height"), ~ mean(.x, na.rm = TRUE)))  %>%
  gather("trait", "mean", 4:8) %>% unite(trait2, c("Treatment", "trait")) %>% ungroup() %>%
  dplyr::select(trait2, mean) -> cwm_means

trt1 <- c("EMV", "LDMC", "Plant.height", "Seed.mass", "SLA")
trt2 <- c("Ellenberg moisture value", "Leaf dry matter content", "Plant height",
          "Seed mass", "Specific leaf area")
as.data.frame(cbind(trt1, trt2)) %>% dplyr::rename(trait = trt1) -> trt2

cwms %>% group_by(Inoculated, Water_table, Treatment) %>%
  dplyr::summarise(across(c("EMV":"Plant.height"), ~ std.error(.x, na.rm = TRUE)))  %>%
  gather("trait", "se", 4:8) %>% unite(trait2, c("Treatment", "trait"), remove=FALSE) %>% 
  full_join(cwm_means, by="trait2") %>% dplyr::select(Inoculated, Water_table, trait, mean, se) %>%
  left_join(trt2, by="trait") -> cwm_vals

cwms %>% gather("trait", "value", 4:8) %>% left_join(trt2, by="trait")-> cwm_inds

ggplot(aes(x=Inoculated, y=mean, fill=Water_table, colour=Water_table), data= cwm_vals) +
  geom_pointrange(aes(ymax=mean+se,
                      ymin=mean-se), position = position_dodge(0.4)) +
  geom_jitter(aes(x=Inoculated, y=value, fill=Water_table, colour=Water_table, shape=Inoculated), data=cwm_inds,
              alpha=0.6, position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.4)) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  guides(color=guide_legend("Water table"), fill = "none") +
  scale_shape_manual(values=c(1, 19)) +
  ylab("Community weighted mean") +
  theme_bw() +
  facet_wrap(~trt2, scales = "free_y", strip.position = "top") +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text=element_text(size=12))