---
title: "Van Wert et al 2023 Fig 4"
author: "Jacey C. Van Wert"
output: html_document
date: "2023-03-28"
---

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

library(patchwork) 
library(ggplot2)
library(vegan)
library(dplyr)
library(data.table)
library(plyr)
library(DescTools)
library(RColorBrewer)
library(devtools)
library(factoextra) 
library(here)
library(tidyverse)
library(ggfortify) 
library(Hmisc)
library(ecodist)
library(vegan)
library(pairwiseAdonis) #devtools::install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
library(knitr)
library(kableExtra)
```

```{r}
new_palette_10_fresh <- c("#c8cb7e","#9a9c61","#6e6f45","#66A61E","#bad80a","#91d3c5", "#35a5d5", "#ff90b5","#ffad2d","#E6AB02","#A6761D", "#EE7E09")
new_palette_8_fresh2 <- c("#9E0142", "#D23B4E", "#FB9A99","#A6CEE3", "#38A1A0", "#3990B9", "#9a99e6")
new_palette_4_fresh <- c("#c8cb7e", "#9a9c61","#ffad2d","#6e6f45")
new_palette_8_fresh3 <- c("#D23B4E", "#FB9A99","#A6CEE3", "#38A1A0", "#3990B9", "#9a99e6")
```

## Upload data
```{r}
j_dat<- read.csv(here::here("data","VanWert_etal_2023_fishfeces.csv"))

j_dat2 <- j_dat %>% 
  dplyr::rename(Ca = ca_ppm, Mg = mg_ppm, Cu = cu_ppm, Mn = mn_ppm, Zn = zn_ppm, Fe = fe_ppm, Protein = protein_percent, Carb = carb_percent, Lipid = lipid_percent) %>% 
  dplyr::mutate(genus_species_sh = case_when(
    genus_species %in% "Acanthurus nigricans" ~ "A. nigricans",
    genus_species %in% "Acanthurus nigrofuscus" ~ "A. nigrofuscus",
    genus_species %in% "Acanthurus triostegus" ~ "A. triostegus",
    genus_species %in% "Acanthurus olivaceus" ~ "A. olivaceus",
    genus_species %in% "Chlorurus spilurus" ~ "C. spilurus",
    genus_species %in% "Amanses scopas" ~ "A. scopas",
    genus_species %in% "Melichthys vidua" ~ "M. vidua",
    genus_species %in% "Naso lituratus" ~ "N. lituratus",
    genus_species %in% "Plectroglyphidodon lacrymatus" ~ "P. lacrymatus",
    genus_species %in% "Scarus psittacus" ~ "S. psittacus",
    genus_species %in% "Stegastes nigricans" ~ "S. nigricans",
    genus_species %in% "Zebrasoma scopas" ~ "Z. scopas",
    genus_species %in% "Chaetodon citrinellus" ~ "C. citrinellus",
    genus_species %in% "Chaetodon ornatissimus" ~ "C. ornatissimus",
    genus_species %in% "Chaetodon quadrimaculatus" ~ "C. quadrimaculatus",
    genus_species %in% "Chaetodon reticulatus" ~ "C. reticulatus",
    genus_species %in% "Chaetodon vagabundus" ~ "C. vagabundus",
    genus_species %in% "Chaetodon vagabundus" ~ "C. vagabundus",
    genus_species %in% "Heniochus chrysostomus" ~ "H. chrysostomus",
    genus_species %in% "Chaetodon auriga" ~ "C. auriga",
  ))

j_dat2$new_genus_species <- j_dat2$genus_species
j_dat2 <- separate(j_dat2, new_genus_species, into = c("genus", "species"), sep = "\\s")
j_dat2 = as.data.frame(j_dat2)

## Merging the datasets for the micronutrients
j_data_merged <- subset(j_dat2, study == "ucsb_2018" | study == "nina_ucsb" | study=="nina_umich")
```



##### HERBIVORE #####
```{r}
j_data_herbi <- subset(j_data_merged, trophic_group == "Herbivore")
j_data_herb2 <- j_data_herbi[c("genus_species_sh", "site", "Ca", "Cu", "Fe", "Mg", "Mn", "Zn")]
j_data_herb2 <- j_data_herb2 %>%
  mutate(ca_trans = log(Ca + 1)) %>%
  mutate(mg_trans = log(Mg + 1)) %>%
  mutate(cu_trans = log(Cu + 1)) %>%
  mutate(mn_trans = log(Mn + 1)) %>%
  mutate(zn_trans = log(Zn + 1)) %>%
  mutate(fe_trans = log(Fe + 1)) %>%
  select("genus_species_sh", "site", "ca_trans", "mg_trans", "cu_trans", "mn_trans", "zn_trans", "fe_trans") %>%
  dplyr::rename(Ca = ca_trans, Mg = mg_trans, Cu = cu_trans, Mn = mn_trans, Zn = zn_trans, Fe = fe_trans)

## Remove na values to compute the analyses
df_na_herb <- (na.omit(j_data_herb2))
com_herb <- df_na_herb[c("Ca", "Cu", "Fe", "Mg", "Mn", "Zn")]
com_herb_sub <- com_herb[rowSums(com_herb[]) > 0, ] ## remove rows containing 0 values (row 494)
env_herb <- df_na_herb %>% select(genus_species_sh)
envsub_herb <- subset(env_herb, rownames(env_herb) %in% rownames(com_herb_sub)) # subset the same rows as com data

m_com_herb = as.matrix(com_herb_sub)
set.seed(123)
nmds_herb = metaMDS(m_com_herb, distance = "bray", trymax=999) 
nmds_herb$stress #0.02477057

herb.spp.fit <- envfit(nmds_herb, (m_com_herb), permutations = 999)
head(herb.spp.fit)


##with ggplot
site.scrs.herb <- as.data.frame(scores(nmds_herb, display = "sites")) 
site.scrs.herb <- cbind(site.scrs.herb, genus_species_sh = envsub_herb$genus_species_sh)

## spp here means micronut
spp.scrs.herb <- as.data.frame(scores(herb.spp.fit, display = "vectors")) 
spp.scrs.herb <- cbind(spp.scrs.herb, Species = rownames(spp.scrs.herb))
spp.scrs.herb <- cbind(spp.scrs.herb, pval = herb.spp.fit$vectors$pvals) 
sig.spp.scrs.herb <- subset(spp.scrs.herb, pval<=0.05) 


## add polygons
grp.a <- site.scrs.herb[site.scrs.herb$genus_species_sh == "A. nigricans", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                    "A. nigricans", c("NMDS1", "NMDS2")]), ] # hull values for grp A
grp.b <- site.scrs.herb[site.scrs.herb$genus_species_sh == "A. nigrofuscus", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                      "A. nigrofuscus", c("NMDS1", "NMDS2")]), ] # hull values for grp B
grp.c <- site.scrs.herb[site.scrs.herb$genus_species_sh == "A. triostegus", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                     "A. triostegus", c("NMDS1", "NMDS2")]), ] # hull values for grp C
grp.d <- site.scrs.herb[site.scrs.herb$genus_species_sh == "C. spilurus", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                   "C. spilurus", c("NMDS1", "NMDS2")]), ] # hull values for grp D
grp.e <- site.scrs.herb[site.scrs.herb$genus_species_sh == "M. vidua", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                "M. vidua", c("NMDS1", "NMDS2")]), ] # hull values for grp E
grp.f <- site.scrs.herb[site.scrs.herb$genus_species_sh == "N. lituratus", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                    "N. lituratus", c("NMDS1", "NMDS2")]), ] # hull values for grp F
grp.g <- site.scrs.herb[site.scrs.herb$genus_species_sh == "P. lacrymatus", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                     "P. lacrymatus", c("NMDS1", "NMDS2")]), ] # hull values for grp G
grp.h <- site.scrs.herb[site.scrs.herb$genus_species_sh == "S. nigricans", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                    "S. nigricans", c("NMDS1", "NMDS2")]), ] # hull values for grp H
grp.i <- site.scrs.herb[site.scrs.herb$genus_species_sh == "S. psittacus", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                    "S. psittacus", c("NMDS1", "NMDS2")]), ] # hull values for grp I
grp.j <- site.scrs.herb[site.scrs.herb$genus_species_sh == "Z. scopas", ][chull(site.scrs.herb[site.scrs.herb$genus_species_sh ==
                                                                                                 "Z. scopas", c("NMDS1", "NMDS2")]), ] # hull values for grp J
hull.data.herb <- rbind(grp.a, grp.b, grp.c, grp.d, grp.e, grp.f, grp.g, grp.h, grp.i, grp.j)
hull.data.herb


nmds.herb <- ggplot(site.scrs.herb, aes(x=NMDS1, y=NMDS2))+
  geom_point(aes(NMDS1, NMDS2, colour = factor(site.scrs.herb$genus_species_sh), shape = factor(site.scrs.herb$genus_species_sh)), size = 5)+ 
  coord_fixed()+
  theme(panel.background = element_rect(fill = NA, colour = "black", size = 1, linetype = "solid"))+
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10))+
  geom_segment(data = sig.spp.scrs.herb, aes(x = 0, xend=NMDS1/2, y=0, yend=NMDS2/2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", lwd=0.3) + 
  ggrepel::geom_text_repel(data = sig.spp.scrs.herb, aes(x=NMDS1/2, y=NMDS2/2, label = Species), cex = 6, direction = "both", segment.size = 0.25, nudge_y = 0.03, nudge_x = 0.065, segment.color = 'transparent')+ 
  scale_fill_manual(values=new_palette_10_fresh, "Species")+
  scale_color_manual(values=new_palette_10_fresh, "Species")+
  scale_shape_manual(values=c(15:25),"Species")+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        plot.title = element_text(hjust = 0.5, face = "bold"), 
        text = element_text(size=20), 
        legend.text = element_text(size = 20, face = "italic"))+
  theme(aspect.ratio=1)+
  geom_polygon(data=hull.data.herb, aes(x=NMDS1,y=NMDS2,fill=genus_species_sh,group=genus_species_sh),alpha=0.30)+
  ggtitle( "Herbivores")

#PERMANOVA
herb_veg <- vegdist((com_herb_sub), method="bray")
j_pca_herb_merged <- df_na_herb[c('genus_species_sh','site','Ca','Mg','Cu','Mn','Zn','Fe')]
j_pca_herb_merged_new <- j_pca_herb_merged[-c(72), ]
fish_herb_perma <- adonis2(herb_veg ~ genus_species_sh, data = j_pca_herb_merged_new, permutations = 999, method="bray", strata = j_pca_herb_merged_new$site)

herb_pairwise <- pairwise.adonis(herb_veg, j_pca_herb_merged_new$genus_species_sh, sim.method = "bray", p.adjust.m = "holm")
herb_pairwise_table <- kbl(herb_pairwise)
```


###### CORALLIVORE #####
```{r}
j_data_coralli <- subset(j_data_merged, trophic_group == "Corallivore") %>%
  filter(genus_species != "Chaetodon auriga") %>%
  filter(genus_species != "Arothron meleagris")
j_data_coralli2 <- j_data_coralli[c("genus_species_sh", "site", "Ca", "Cu", "Fe", "Mg", "Mn", "Zn")]
j_data_coralli2 <- j_data_coralli2 %>%
  mutate(ca_trans = log(Ca + 1)) %>%
  mutate(mg_trans = log(Mg + 1)) %>%
  mutate(cu_trans = log(Cu + 1)) %>%
  mutate(mn_trans = log(Mn + 1)) %>%
  mutate(zn_trans = log(Zn + 1)) %>%
  mutate(fe_trans = log(Fe + 1)) %>%
  select("genus_species_sh", "site", "ca_trans", "mg_trans", "cu_trans", "mn_trans", "zn_trans", "fe_trans") %>%
  dplyr::rename(Ca = ca_trans, Mg = mg_trans, Cu = cu_trans, Mn = mn_trans, Zn = zn_trans, Fe = fe_trans)

## Remove na values to compute the analyses
df_na_coralli <- (na.omit(j_data_coralli2))
com_coralli <- df_na_coralli[c("Ca", "Cu", "Fe", "Mg", "Mn", "Zn")]
#com_coralli_sub <- com_coralli[rowSums(com_coralli[]) > 0, ] ## remove rows containing 0 values, none
env_coralli <- df_na_coralli %>% select(genus_species_sh)
envsub_coralli <- subset(env_coralli, rownames(env_coralli) %in% rownames(com_coralli)) # subset the same rows as com data

m_com_coralli <- as.matrix(com_coralli)
set.seed(123)
nmds_coralli <- metaMDS(m_com_coralli, distance = "bray", trymax = 999)
nmds_coralli$stress # 0.118

coralli.spp.fit <- envfit(nmds_coralli, (m_com_coralli), permutations = 999)
head(coralli.spp.fit)


## with ggplot
site.scrs.coralli <- as.data.frame(scores(nmds_coralli, display = "sites"))
site.scrs.coralli <- cbind(site.scrs.coralli, genus_species_sh = envsub_coralli$genus_species_sh)

## spp here means micronut
spp.scrs.coralli <- as.data.frame(scores(coralli.spp.fit, display = "vectors"))
spp.scrs.coralli <- cbind(spp.scrs.coralli, Species = rownames(spp.scrs.coralli))
spp.scrs.coralli <- cbind(spp.scrs.coralli, pval = coralli.spp.fit$vectors$pvals)
sig.spp.scrs.coralli <- subset(spp.scrs.coralli, pval <= 0.05)
head(spp.scrs.coralli)


## add polygons
grp.a <- site.scrs.coralli[site.scrs.coralli$genus_species_sh == "A. scopas", ][chull(site.scrs.coralli[site.scrs.coralli$genus_species_sh ==
                                                                                                          "A. scopas", c("NMDS1", "NMDS2")]), ] # hull values for grp A
grp.b <- site.scrs.coralli[site.scrs.coralli$genus_species_sh == "C. citrinellus", ][chull(site.scrs.coralli[site.scrs.coralli$genus_species_sh ==
                                                                                                               "C. citrinellus", c("NMDS1", "NMDS2")]), ] # hull values for grp B
grp.c <- site.scrs.coralli[site.scrs.coralli$genus_species_sh == "C. ornatissimus", ][chull(site.scrs.coralli[site.scrs.coralli$genus_species_sh ==
                                                                                                                "C. ornatissimus", c("NMDS1", "NMDS2")]), ] # hull values for grp C
grp.d <- site.scrs.coralli[site.scrs.coralli$genus_species_sh == "C. quadrimaculatus", ][chull(site.scrs.coralli[site.scrs.coralli$genus_species_sh ==
                                                                                                                   "C. quadrimaculatus", c("NMDS1", "NMDS2")]), ] # hull values for grp D
grp.e <- site.scrs.coralli[site.scrs.coralli$genus_species_sh == "C. reticulatus", ][chull(site.scrs.coralli[site.scrs.coralli$genus_species_sh ==
                                                                                                               "C. reticulatus", c("NMDS1", "NMDS2")]), ] # hull values for grp E
grp.f <- site.scrs.coralli[site.scrs.coralli$genus_species_sh == "C. vagabundus", ][chull(site.scrs.coralli[site.scrs.coralli$genus_species_sh ==
                                                                                                              "C. vagabundus", c("NMDS1", "NMDS2")]), ] # hull values for grp F
grp.g <- site.scrs.coralli[site.scrs.coralli$genus_species_sh == "H. chrysostomus", ][chull(site.scrs.coralli[site.scrs.coralli$genus_species_sh ==
                                                                                                                "H. chrysostomus", c("NMDS1", "NMDS2")]), ] # hull values for grp G

hull.data.coralli <- rbind(grp.a, grp.b, grp.c, grp.d, grp.e, grp.f, grp.g)
hull.data.coralli


nmds.coralli <- ggplot(site.scrs.coralli, aes(x = NMDS1, y = NMDS2)) +
  geom_point(aes(NMDS1, NMDS2, colour = factor(site.scrs.coralli$genus_species_sh), shape = factor(site.scrs.coralli$genus_species_sh)), size = 5) +
  coord_fixed() +
  theme(panel.background = element_rect(fill = NA, colour = "black", size = 1, linetype = "solid")) +
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10)) +
  geom_segment(data = sig.spp.scrs.coralli, aes(x = 0, xend = NMDS1/2, y = 0, yend = NMDS2/2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", lwd = 0.3) +
  ggrepel::geom_text_repel(data = sig.spp.scrs.coralli, aes(x = NMDS1/2, y = NMDS2/2, label = Species), cex = 6, direction = "both", segment.size = 0.25, nudge_x = 0.05) +
  scale_fill_manual(values = new_palette_8_fresh2, "Species") +
  scale_color_manual(values = new_palette_8_fresh2, "Species") +
  scale_shape_manual(values = c(15:24), "Species") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    text = element_text(size = 20),
    legend.text = element_text(size = 20, face = "italic")
  ) +
  theme(aspect.ratio = 1) +
  geom_polygon(data = hull.data.coralli, aes(x = NMDS1, y = NMDS2, fill = genus_species_sh, group = genus_species_sh), alpha = 0.30) +
  ggtitle("Corallivores") 

#PERMANOVA
coral_veg <- vegdist((com_coralli), method="bray")
j_pca_coral_merged <- df_na_coralli[c('genus_species_sh','site','Ca','Mg','Cu','Mn','Zn','Fe')]
fish_coral_perma <- adonis2(coral_veg ~ genus_species_sh, data = j_pca_coral_merged, permutations = 999, method="bray", strata = j_pca_coral_merged$site)

coral_pairwise <- pairwise.adonis(coral_veg, j_pca_coral_merged$genus_species_sh, sim.method = "bray", p.adjust.m = "holm")
coral_pairwise_table <- kbl(coral_pairwise)
```


###### ACANTHURUS####
```{r}
j_data_acanthurus <- subset(j_data_merged, genus_species %in% c("Acanthurus nigricans", "Acanthurus nigrofuscus", "Acanthurus olivaceus", "Acanthurus triostegus"))
j_data_acanthurus2 <- j_data_acanthurus[c("genus_species_sh", "site", "Ca", "Cu", "Fe", "Mg", "Mn", "Zn")]
j_data_acanthurus2 <- j_data_acanthurus2 %>%
  mutate(ca_trans = log(Ca + 1)) %>%
  mutate(mg_trans = log(Mg + 1)) %>%
  mutate(cu_trans = log(Cu + 1)) %>%
  mutate(mn_trans = log(Mn + 1)) %>%
  mutate(zn_trans = log(Zn + 1)) %>%
  mutate(fe_trans = log(Fe + 1)) %>%
  select("genus_species_sh", "site", "ca_trans", "mg_trans", "cu_trans", "mn_trans", "zn_trans", "fe_trans") %>%
  dplyr::rename(Ca = ca_trans, Mg = mg_trans, Cu = cu_trans, Mn = mn_trans, Zn = zn_trans, Fe = fe_trans)

## Remove na values to compute the analyses
df_na_aca <- (na.omit(j_data_acanthurus2))
com_acan <- df_na_aca[c("Ca", "Cu", "Fe", "Mg", "Mn", "Zn")]
#com_acan_sub <- com_acan[rowSums(com_acan[]) > 0, ] ## remove rows containing 0 values, none
env_acan <- df_na_aca %>% select(genus_species_sh)
envsub_acan <- subset(env_acan, rownames(env_acan) %in% rownames(com_acan)) # subset the same rows as com data

m_com_acan <- as.matrix(com_acan)
set.seed(123)
nmds_acan <- metaMDS(m_com_acan, distance = "bray", trymax = 999)
nmds_acan$stress # 0.041

acan.spp.fit <- envfit(nmds_acan, (m_com_acan), permutations = 999)
head(acan.spp.fit)


## with ggplot
site.scrs.acan <- as.data.frame(scores(nmds_acan, display = "sites"))
site.scrs.acan <- cbind(site.scrs.acan, genus_species_sh = envsub_acan$genus_species_sh)

## spp here means micronut
spp.scrs.acan <- as.data.frame(scores(acan.spp.fit, display = "vectors"))
spp.scrs.acan <- cbind(spp.scrs.acan, Species = rownames(spp.scrs.acan))
spp.scrs.acan <- cbind(spp.scrs.acan, pval = acan.spp.fit$vectors$pvals)
sig.spp.scrs.acan <- subset(spp.scrs.acan, pval <= 0.05)
head(spp.scrs.acan)


## add polygons
grp.a <- site.scrs.acan[site.scrs.acan$genus_species_sh == "A. nigricans", ][chull(site.scrs.acan[site.scrs.acan$genus_species_sh ==
                                                                                                    "A. nigricans", c("NMDS1", "NMDS2")]), ] # hull values for grp A
grp.b <- site.scrs.acan[site.scrs.acan$genus_species_sh == "A. nigrofuscus", ][chull(site.scrs.acan[site.scrs.acan$genus_species_sh ==
                                                                                                      "A. nigrofuscus", c("NMDS1", "NMDS2")]), ] # hull values for grp B
grp.c <- site.scrs.acan[site.scrs.acan$genus_species_sh == "A. olivaceus", ][chull(site.scrs.acan[site.scrs.acan$genus_species_sh ==
                                                                                                    "A. olivaceus", c("NMDS1", "NMDS2")]), ] # hull values for grp C
grp.d <- site.scrs.acan[site.scrs.acan$genus_species_sh == "A. triostegus", ][chull(site.scrs.acan[site.scrs.acan$genus_species_sh ==
                                                                                                     "A. triostegus", c("NMDS1", "NMDS2")]), ] # hull values for grp D

hull.data.acan <- rbind(grp.a, grp.b, grp.c, grp.d) 
hull.data.acan


nmds.acan <- ggplot(site.scrs.acan, aes(x=NMDS1, y=NMDS2))+
  geom_point(aes(NMDS1, NMDS2, colour = factor(site.scrs.acan$genus_species_sh), shape = factor(site.scrs.acan$genus_species_sh)), size = 5)+ 
  coord_fixed()+
  theme(panel.background = element_rect(fill = NA, colour = "black", size = 1, linetype = "solid"))+
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10),)+  
  geom_segment(data = sig.spp.scrs.acan, aes(x = 0, xend=NMDS1/2, y=0, yend=NMDS2/2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", lwd=0.3) + 
  ggrepel::geom_text_repel(data = sig.spp.scrs.acan, aes(x=NMDS1/2, y=NMDS2/2, label = Species), cex = 6, direction = "both", segment.size = 0.25, nudge_y = 0.025, nudge_x = 0.04, segment.color = 'transparent')+ 
  scale_fill_manual(values=new_palette_4_fresh, "Species")+
  scale_color_manual(values=new_palette_4_fresh, "Species")+
  scale_shape_manual(values=c(15:25),"Species")+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        plot.title = element_text(hjust = 0.5, face = "bold.italic"), 
        text = element_text(size=20), 
        legend.text = element_text(size = 20, face = "italic"))+
  theme(aspect.ratio=1)+
  geom_polygon(data=hull.data.acan,aes(x=NMDS1,y=NMDS2,fill=genus_species_sh,group=genus_species_sh),alpha=0.30)+  
  ggtitle("Acanthurus")+
 scale_y_continuous(breaks = c(-0.25, 0, 0.25))+
 scale_x_continuous(breaks = c(-0.25, 0, 0.25))


#PERMANOVA
acan_veg <- vegdist((com_acan), method="bray")
j_pca_acan_merged <- df_na_aca[c('genus_species_sh','site','Ca','Mg','Cu','Mn','Zn','Fe')]
fish_acan_perma <- adonis2(acan_veg ~ genus_species_sh, data = j_pca_acan_merged, permutations = 999, method="bray", strata = j_pca_acan_merged$site)

acan_pairwise <- pairwise.adonis(acan_veg, j_pca_acan_merged$genus_species_sh, sim.method = "bray", p.adjust.m = "holm")
acan_pairwise_table <- kbl(acan_pairwise)
```


### CHAETODON #######
```{r}
j_data_chaetodon <- subset(j_data_merged, genus_species %in% c("Chaetodon ornatissimus", "Chaetodon citrinellus", "Chaetodon quadrimaculatus", "Chaetodon reticulatus", "Chaetodon vagabundus"))
j_data_chaetodon2 <- j_data_chaetodon[c("genus_species_sh", "site", "Ca", "Cu", "Fe", "Mg", "Mn", "Zn")]
j_data_chaetodon2 <- j_data_chaetodon2 %>%
  mutate(ca_trans = log(Ca + 1)) %>%
  mutate(mg_trans = log(Mg + 1)) %>%
  mutate(cu_trans = log(Cu + 1)) %>%
  mutate(mn_trans = log(Mn + 1)) %>%
  mutate(zn_trans = log(Zn + 1)) %>%
  mutate(fe_trans = log(Fe + 1)) %>%
  select("genus_species_sh", "site", "ca_trans", "mg_trans", "cu_trans", "mn_trans", "zn_trans", "fe_trans") %>%
  dplyr::rename(Ca = ca_trans, Mg = mg_trans, Cu = cu_trans, Mn = mn_trans, Zn = zn_trans, Fe = fe_trans)

## Remove na values to compute the analyses
df_na_chaet <- (na.omit(j_data_chaetodon2))
com_chaet <- df_na_chaet[c("Ca", "Mg", "Cu", "Mn", "Zn", "Fe")]
#com_chaet_sub <- com_chaet[rowSums(com_chaet[]) > 0, ] ## remove rows containing 0 values, none
env_chaet <- df_na_chaetodon_spe %>% select(genus_species_sh)
envsub_chaet <- subset(env_chaet, rownames(env_chaet) %in% rownames(com_chaet)) # subset the same rows as com data

m_com_chaet <- as.matrix(com_chaet)
set.seed(123)
nmds_chaet <- metaMDS(m_com_chaet, distance = "bray", trymax = 999)
nmds_chaet$stress # 0.078

chaet.spp.fit <- envfit(nmds_chaet, (m_com_chaet), permutations = 999)
head(chaet.spp.fit)


## with ggplot
site.scrs.chaet <- as.data.frame(scores(nmds_chaet, display = "sites")) # save NMDS results into dataframe
site.scrs.chaet <- cbind(site.scrs.chaet, genus_species_sh = envsub_chaet$genus_species_sh) # add grouping variable tropic group to dataframe

## spp here means micronut
spp.scrs.chaet <- as.data.frame(scores(chaet.spp.fit, display = "vectors"))
spp.scrs.chaet <- cbind(spp.scrs.chaet, Species = rownames(spp.scrs.chaet))
spp.scrs.chaet <- cbind(spp.scrs.chaet, pval = chaet.spp.fit$vectors$pvals)
sig.spp.scrs.chaet <- subset(spp.scrs.chaet, pval <= 0.05)


## add polygons
grp.a <- site.scrs.chaet[site.scrs.chaet$genus_species_sh == "C. citrinellus", ][chull(site.scrs.chaet[site.scrs.chaet$genus_species_sh ==
                                                                                                         "C. citrinellus", c("NMDS1", "NMDS2")]), ] # hull values for grp A
grp.b <- site.scrs.chaet[site.scrs.chaet$genus_species_sh == "C. ornatissimus", ][chull(site.scrs.chaet[site.scrs.chaet$genus_species_sh ==
                                                                                                          "C. ornatissimus", c("NMDS1", "NMDS2")]), ] # hull values for grp B
grp.c <- site.scrs.chaet[site.scrs.chaet$genus_species_sh == "C. quadrimaculatus", ][chull(site.scrs.chaet[site.scrs.chaet$genus_species_sh ==
                                                                                                             "C. quadrimaculatus", c("NMDS1", "NMDS2")]), ] # hull values for grp C
grp.d <- site.scrs.chaet[site.scrs.chaet$genus_species_sh == "C. reticulatus", ][chull(site.scrs.chaet[site.scrs.chaet$genus_species_sh ==
                                                                                                         "C. reticulatus", c("NMDS1", "NMDS2")]), ] # hull values for grp D
grp.e <- site.scrs.chaet[site.scrs.chaet$genus_species_sh == "C. vagabundus", ][chull(site.scrs.chaet[site.scrs.chaet$genus_species_sh ==
                                                                                                        "C. vagabundus", c("NMDS1", "NMDS2")]), ] # hull values for grp D
hull.data.chaet <- rbind(grp.a, grp.b, grp.c, grp.d, grp.e)
hull.data.chaet


nmds.chaet <- ggplot(site.scrs.chaet, aes(x = NMDS1, y = NMDS2)) +
  geom_point(aes(NMDS1, NMDS2, colour = factor(site.scrs.chaet$genus_species_sh), shape = factor(site.scrs.chaet$genus_species_sh)), size = 5) +
  coord_fixed() +
      ggrepel::geom_text_repel(data = sig.spp.scrs.chaet, aes(x = NMDS1/2, y = NMDS2/2, label = Species), cex = 6, direction = "both", segment.size = 0.25, nudge_y = 0.04, nudge_x = -0.01, segment.color = 'transparent') +
  theme(panel.background = element_rect(fill = NA, colour = "black", size = 1, linetype = "solid")) +
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), ) +
  geom_segment(data = sig.spp.scrs.chaet, aes(x = 0, xend = NMDS1/2, y = 0, yend = NMDS2/2), arrow = arrow(length = unit(0.1, "cm")), colour = "grey10", lwd = 0.3) +
  scale_fill_manual(values = new_palette_8_fresh3, "Species") +
  scale_color_manual(values = new_palette_8_fresh3, "Species") +
  scale_shape_manual(values = c(15:25), "Species") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"),
    plot.title = element_text(hjust = 0.5, face = "bold.italic"),
    text = element_text(size = 20),
    legend.text = element_text(size = 20, face = "italic")
  ) +
  theme(aspect.ratio = 1) +
  geom_polygon(data = hull.data.chaet, aes(x = NMDS1, y = NMDS2, fill = genus_species_sh, group = genus_species_sh), alpha = 0.30) +
  ggtitle("Chaetodon") + 
  scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4))+# limits = c(-.75, 1))+
  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4))# limits = c(-.75, 1))

#PERMANOVA
chaet_veg <- vegdist((com_chaet), method="bray")
j_pca_chaet_merged <- df_na_chaet[c('genus_species_sh','site','Ca','Mg','Cu','Mn','Zn','Fe')]
fish_chaet_perma <- adonis2(chaet_veg ~ genus_species_sh, data = j_pca_chaet_merged, permutations = 999, method="bray", strata = j_pca_chaet_merged$site)

chaet_pairwise <- pairwise.adonis(chaet_veg, j_pca_chaet_merged$genus_species_sh, sim.method = "bray", p.adjust.m = "holm")
chaet_pairwise_table <- kbl(chaet_pairwise)
```


##Combine plots
```{r}
(nmds.herb + nmds.coralli)/(nmds.acan + nmds.chaet) + plot_annotation(tag_levels = 'a') & theme(plot.tag = element_text(face = 'bold', size = 20))
ggsave("fig4.jpg",height=12,width=18)
```

