library(openxlsx)
library(ggplot2)
library(ggpubr)
library(lemon)

dados<-read.xlsx("Suplemental Table S1.xlsx")

dados$GROUP<-as.factor(dados$TAXON_GROUP)
dados$DIET<-as.factor(dados$DIET)
dados$CONTINENT<-as.factor(dados$CONTINENT)
dados$SYSTEM<-as.factor(dados$SYSTEM)

boxplot.stats(dados$Mean_d13C[dados$SYSTEM=="captive"])$out
boxplot.stats(dados$Mean_d13C[dados$SYSTEM=="wild"])$out
boxplot.stats(dados$Mean_d15N[dados$SYSTEM=="captive"])$out
boxplot.stats(dados$Mean_d15N[dados$SYSTEM=="wild"])$out
boxplot.stats(dados$Mean_d2H[dados$SYSTEM=="captive"])$out
boxplot.stats(dados$Mean_d2H[dados$SYSTEM=="wild"])$out
boxplot.stats(dados$Mean_d18O[dados$SYSTEM=="captive"])$out
boxplot.stats(dados$Mean_d18O[dados$SYSTEM=="wild"])$out
boxplot.stats(dados$Mean_d34S[dados$SYSTEM=="captive"])$out
boxplot.stats(dados$Mean_d34S[dados$SYSTEM=="wild"])$out

tapply(dados$Mean_d13C,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Mean_d13C,dados$SYSTEM,sd,na.rm=T)
tapply(dados$Mean_d15N,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Mean_d15N,dados$SYSTEM,sd,na.rm=T)
tapply(dados$Mean_d2H,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Mean_d2H,dados$SYSTEM,sd,na.rm=T)
tapply(dados$Mean_d18O,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Mean_d18O,dados$SYSTEM,sd,na.rm=T)
tapply(dados$Mean_d34S,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Mean_d34S,dados$SYSTEM,sd,na.rm=T)

tapply(dados$SD_d13C,dados$SYSTEM,mean,na.rm=T)
tapply(dados$SD_d13C,dados$SYSTEM,sd,na.rm=T)
tapply(dados$SD_d15N,dados$SYSTEM,mean,na.rm=T)
tapply(dados$SD_d15N,dados$SYSTEM,sd,na.rm=T)
tapply(dados$SD_d2H,dados$SYSTEM,mean,na.rm=T)
tapply(dados$SD_d2H,dados$SYSTEM,sd,na.rm=T)
tapply(dados$SD_d18O,dados$SYSTEM,mean,na.rm=T)
tapply(dados$SD_d18O,dados$SYSTEM,sd,na.rm=T)
tapply(dados$SD_d34S,dados$SYSTEM,mean,na.rm=T)
tapply(dados$SD_d34S,dados$SYSTEM,sd,na.rm=T)

tapply(dados$Range_d13C,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Range_d13C,dados$SYSTEM,sd,na.rm=T)
tapply(dados$Range_d15N,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Range_d15N,dados$SYSTEM,sd,na.rm=T)
tapply(dados$Range_d2H,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Range_d2H,dados$SYSTEM,sd,na.rm=T)
tapply(dados$Range_d18O,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Range_d18O,dados$SYSTEM,sd,na.rm=T)
tapply(dados$Range_d34S,dados$SYSTEM,mean,na.rm=T)
tapply(dados$Range_d34S,dados$SYSTEM,sd,na.rm=T)

t.test(Mean_d13C~SYSTEM,dados)
t.test(Mean_d15N~SYSTEM,dados)
t.test(Mean_d2H~SYSTEM,dados)
t.test(Mean_d18O~SYSTEM,dados)
t.test(Mean_d34S~SYSTEM,dados)

t.test(SD_d13C~SYSTEM,dados)
t.test(SD_d15N~SYSTEM,dados)
t.test(SD_d2H~SYSTEM,dados)
t.test(SD_d18O~SYSTEM,dados)
t.test(SD_d34S~SYSTEM,dados)

t.test(Range_d13C~SYSTEM,dados)
t.test(Range_d15N~SYSTEM,dados)
t.test(Range_d2H~SYSTEM,dados)
t.test(Range_d18O~SYSTEM,dados)
t.test(Range_d34S~SYSTEM,dados)

Afr <- subset(dados, CONTINET =="Africa")
Asi <- subset(dados, CONTINET =="Asia")
Eur <- subset(dados, CONTINET =="Europe")
Nor <- subset(dados, CONTINET =="North America")
Sou <- subset(dados, CONTINET =="South America")
Oce <- subset(dados, CONTINET =="Oceania")

t.test(Mean_d13C~SYSTEM,Afr) 
t.test(Mean_d15N~SYSTEM,Afr) 
t.test(Mean_d13C~SYSTEM,Asi) 
t.test(Mean_d15N~SYSTEM,Asi) 
t.test(Mean_d13C~SYSTEM,Eur) 
t.test(Mean_d15N~SYSTEM,Eur) 
t.test(Mean_d13C~SYSTEM,Nor) 
t.test(Mean_d15N~SYSTEM,Nor) 
t.test(Mean_d13C~SYSTEM,Sou)
t.test(Mean_d15N~SYSTEM,Sou) 
t.test(Mean_d13C~SYSTEM,Oce) 
t.test(Mean_d15N~SYSTEM,Oce)

t.test(SD_d13C~SYSTEM,Afr) 
t.test(SD_d15N~SYSTEM,Afr) 
t.test(SD_d13C~SYSTEM,Asi) 
t.test(SD_d15N~SYSTEM,Asi) 
t.test(SD_d13C~SYSTEM,Eur) 
t.test(SD_d15N~SYSTEM,Eur) 
t.test(SD_d13C~SYSTEM,Nor) 
t.test(SD_d15N~SYSTEM,Nor) 
t.test(SD_d13C~SYSTEM,Sou) 
t.test(SD_d15N~SYSTEM,Sou) 
t.test(SD_d13C~SYSTEM,Oce) 
t.test(SD_d15N~SYSTEM,Oce)

t.test(Range_d13C~SYSTEM,Afr) 
t.test(Range_d15N~SYSTEM,Afr) 
t.test(Range_d13C~SYSTEM,Asi) 
t.test(Range_d15N~SYSTEM,Asi) 
t.test(Range_d13C~SYSTEM,Eur) 
t.test(Range_d15N~SYSTEM,Eur) 
t.test(Range_d13C~SYSTEM,Nor) 
t.test(Range_d15N~SYSTEM,Nor) 
t.test(Range_d13C~SYSTEM,Sou) 
t.test(Range_d15N~SYSTEM,Sou) 
t.test(Range_d13C~SYSTEM,Oce) 
t.test(Range_d15N~SYSTEM,Oce)

Fish <- subset(dados, GROUP =="Fish")
Mam <- subset(dados, GROUP =="Mammal")
Bird <- subset(dados, GROUP =="Bird")
Rep <- subset(dados, GROUP =="Reptile")
Amph<-subset(dados, GROUP =="Amphibian")

t.test(Mean_d13C~SYSTEM,Mam)
t.test(Mean_d15N~SYSTEM,Mam)
t.test(Mean_d13C~SYSTEM,Rep)
t.test(Mean_d15N~SYSTEM,Rep)
t.test(Mean_d13C~SYSTEM,Bird)
t.test(Mean_d15N~SYSTEM,Bird)
t.test(Mean_d13C~SYSTEM,Amph)
t.test(Mean_d15N~SYSTEM,Amph)
t.test(Mean_d13C~SYSTEM,Fish)
t.test(Mean_d15N~SYSTEM,Fish)

t.test(SD_d13C~SYSTEM,Mam)
t.test(SD_d15N~SYSTEM,Mam)
t.test(SD_d13C~SYSTEM,Rep)
t.test(SD_d15N~SYSTEM,Rep)
t.test(SD_d13C~SYSTEM,Bird) 
t.test(SD_d15N~SYSTEM,Bird) 
t.test(SD_d13C~SYSTEM,Amph) 
t.test(SD_d15N~SYSTEM,Amph) 
t.test(SD_d13C~SYSTEM,Fish) 
t.test(SD_d15N~SYSTEM,Fish)

t.test(Range_d13C~SYSTEM,Mam) 
t.test(Range_d15N~SYSTEM,Mam) 
t.test(Range_d13C~SYSTEM,Rep) 
t.test(Range_d15N~SYSTEM,Rep) 
t.test(Range_d13C~SYSTEM,Bird)
t.test(Range_d15N~SYSTEM,Bird)
t.test(Range_d13C~SYSTEM,Amph)
t.test(Range_d15N~SYSTEM,Amph)
t.test(Range_d13C~SYSTEM,Fish)
t.test(Range_d15N~SYSTEM,Fish)

fishE<-subset(dados,GROUP=="Fish"&CONTINET=="Europe")
t.test(Mean_d13C~SYSTEM,fishE)
t.test(Mean_d15N~SYSTEM,fishE)
summary(aov(Mean_d15N~SYSTEM+Error(DIET),fishE))
fishO<-subset(dados,GROUP=="Fish"&CONTINET=="South America"|
                GROUP=="Fish"&CONTINET=="North America"|
                GROUP=="Fish"&CONTINET=="Asia"|
                GROUP=="Fish"&CONTINET=="Oceania"|
                GROUP=="Fish"&CONTINET=="Africa")

t.test(Mean_d13C~SYSTEM,fishO) 
t.test(Mean_d15N~SYSTEM,fishO)

RAMB<-subset(dados,GROUP=="Reptile"|GROUP=="Amphibian"|GROUP=="Mammal"|GROUP=="Bird")

boxplot(Mean_d13C~SYSTEM,RAMB)
boxplot(Mean_d15N~SYSTEM,RAMB)

t.test(Mean_d13C~SYSTEM,RAMB) 
t.test(Mean_d15N~SYSTEM,RAMB)

summary(aov(Mean_d15N~SYSTEM+Error(DIET),RAMB))

t.test(Mean_d13C~SYSTEM,Car)
t.test(Mean_d15N~SYSTEM,Car)
t.test(Mean_d13C~SYSTEM,Her)
t.test(Mean_d15N~SYSTEM,Her)
t.test(Mean_d13C~SYSTEM,Omn)
t.test(Mean_d15N~SYSTEM,Omn)

library(ggplot2)
library(ggpubr)

CC=data.frame(dados$SYSTEM,dados$CONTINET,dados$Mean_d13)
CCc=CC[complete.cases(CC),]
BoxCC<-ggplot(CCc, aes(x=dados.CONTINET, y=dados.Mean_d13, colour= dados.SYSTEM, fill=dados.SYSTEM)) + 
  geom_boxplot(alpha = 0.3) +
  #geom_jitter(alpha = 0.6) +
  scale_colour_manual(values = c( "#d7191c", "#2c7bb6")) +
  scale_fill_manual(values = c("#d7191c", "#2c7bb6")) + 
  ylab(expression('Mean '~{delta}^13*C~'(\u2030)'))+
  xlab(label = "")+
  #scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

CN=data.frame(dados$SYSTEM,dados$CONTINET,dados$Mean_d15N)
CNc=CN[complete.cases(CN),]
BoxCN<-ggplot(CNc, aes(x=dados.CONTINET, y=dados.Mean_d15N, colour= dados.SYSTEM, fill=dados.SYSTEM)) + 
  geom_boxplot(alpha = 0.3) +
  #geom_jitter(alpha = 0.6) +
  scale_colour_manual(values = c( "#d7191c", "#2c7bb6")) +
  scale_fill_manual(values = c("#d7191c", "#2c7bb6")) + 
  ylab(expression('Mean '~{delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

GC=data.frame(dados$SYSTEM,dados$GROUP,dados$Mean_d13C)
BoxCG<-ggplot(GC, aes(x=dados$GROUP, y=dados$Mean_d13C, colour= dados.SYSTEM, fill=dados.SYSTEM)) + 
  geom_boxplot(alpha = 0.3) +
  #geom_jitter(alpha = 0.6) +
  scale_colour_manual(name = "Rearing system", values = c( "#d7191c", "#2c7bb6")) +
  scale_fill_manual(name = "Rearing system", values = c("#d7191c", "#2c7bb6")) + 
  ylab(expression('Mean '~{delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = c(.89,.87),
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

GN=data.frame(dados$SYSTEM,dados$GROUP,dados$Mean_d15N)
BoxNG<-ggplot(GN, aes(x=dados$GROUP, y=dados$Mean_d15N, colour= dados.SYSTEM, fill=dados$SYSTEM)) + 
  geom_boxplot(alpha = 0.3) +
  #geom_jitter(alpha = 0.6) +
  scale_colour_manual(values = c( "#d7191c", "#2c7bb6")) +
  scale_fill_manual(values = c("#d7191c", "#2c7bb6")) + 
  ylab(expression('Mean '~{delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

windows(8,12)
ggarrange(BoxCC, BoxCG,BoxCN,BoxNG,ncol=2,nrow=2)
savePlot("Fig._3.pdf", type = "pdf")

BoxCD<-ggplot(dados, aes(x=DIET, y=Mean_d13C,fill=DIET, colour=DIET)) + 
  geom_boxplot(alpha = 0.3) +
  #geom_jitter(alpha = 0.6) +
  scale_colour_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
  scale_fill_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) + 
  ylab(expression('Mean '~{delta}^13*C~'(\u2030)'))+
  xlab(label = "")+ 
  theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

BoxND<-ggplot(dados, aes(x=DIET, y=Mean_d15N,fill=DIET, colour=DIET)) + 
  geom_boxplot(alpha = 0.3) +
  #geom_jitter(alpha = 0.6) +
  scale_colour_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
  scale_fill_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) + 
  ylab(expression('Mean '~{delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

windows(5,10)
ggarrange(BoxCD,BoxND,ncol=2,nrow=1)

savePlot("Fig._4.pdf", type = "pdf")

windows(6,10)
ggplot(dados,aes(x=Mean_d13C,y=Mean_d15N,colour=SYSTEM))+ 
  geom_point(size=2, alpha = 0.4)+
  stat_ellipse()+
  xlab(expression('Mean '~{delta}^13*C~'(\u2030)'))+
  ylab(expression('Mean '~{delta}^15*N~'(\u2030)'))+
  scale_colour_manual(name = "Rearing system", values = c( "#d7191c", "#2c7bb6")) +
  #scale_fill_manual(name = "Rearing system", values = c("#d7191c", "#2c7bb6")) +
  facet_rep_wrap(~GROUP, nrow=2, ncol=3)+  theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = c(.85,.30),
        legend.title = element_text(color = "black", size=16),
        legend.text = element_text(color = "black", size=14),
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))
#quartz.save("Fig Supp_1.pdf", type = "pdf")
savePlot("Fig Supp_1.pdf", type = "pdf")

#### Analyses -----

data <- read.table("2022-11-28_Banco de dados_VLxCat.txt", stringsAsFactors = T, header = T, )

# Cree_1999 ----
iso.cree <- data[which(data$Reference == "Cree_1999"), ]
iso.cree <- droplevels(iso.cree)

table(iso.cree$System)

ggplot(iso.cree, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.cree$d13C[iso.cree$System == "Wild"])
shapiro.test(iso.cree$d13C[iso.cree$System == "Captive"])
var.test(d13C~System, data = iso.cree)
bartlett.test(d13C~System, data = iso.cree)

# Amostra desigual, Na0-Normal e Heterogeneo
#t.test(d13C~System, data = iso.cree, var.equal=T, paired=F)
wilcox.test(d13C~System, data = iso.cree, paired=F)

# Stoskopf_2001 ----
iso.Stoskopf <- data[which(data$Reference == "Stoskopf_2001"), ]
iso.Stoskopf <- droplevels(iso.Stoskopf)

table(iso.Stoskopf$System)

ggplot(iso.Stoskopf, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.Stoskopf$d13C[iso.Stoskopf$System == "Wild"])
shapiro.test(iso.Stoskopf$d13C[iso.Stoskopf$System == "Captive"])
var.test(d13C~System, data = iso.Stoskopf)
bartlett.test(d13C~System, data = iso.Stoskopf)

# Amostra desigual, Não-Normal e Heterogeneo
#t.test(d13C~System, data = iso.Stoskopf, var.equal=T, paired=F)
wilcox.test(d13C~System, data = iso.Stoskopf, paired=F)


# Hammershop_2004 ----
# Hammershop - Claw ----
iso.HammershopClaw <- data[which(data$Reference == "Hammershop_2004" & data$Tissue == "Claw"), ]
iso.HammershopClaw <- droplevels(iso.HammershopClaw)

table(iso.HammershopClaw$System)

ggplot(iso.HammershopClaw, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.HammershopClaw$d13C[iso.HammershopClaw$System == "Wild"])
shapiro.test(iso.HammershopClaw$d13C[iso.HammershopClaw$System == "Captive"])
var.test(d13C~System, data = iso.HammershopClaw)
bartlett.test(d13C~System, data = iso.HammershopClaw)

# Amostra desigual, Normal e Heterogeneo
t.test(d13C~System, data = iso.HammershopClaw, var.equal=F, paired=F)
#wilcox.test(d13C~System, data = iso.Stoskopf, paired=F)

#  Hammershop - Teeth ----
iso.HammershopTeeth <- data[which(data$Reference == "Hammershop_2004" & data$Tissue == "Teeth"), ]
iso.HammershopTeeth <- droplevels(iso.HammershopTeeth)

table(iso.HammershopTeeth$System)

ggplot(iso.HammershopTeeth, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.HammershopTeeth$d13C[iso.HammershopTeeth$System == "Wild"])
shapiro.test(iso.HammershopTeeth$d13C[iso.HammershopTeeth$System == "Captive"])
var.test(d13C~System, data = iso.HammershopTeeth)
bartlett.test(d13C~System, data = iso.HammershopTeeth)

# Amostra desigual, Normal e Heterogeneo
#t.test(d13C~System, data = iso.HammershopTeeth, var.equal=F, paired=F)
wilcox.test(d13C~System, data = iso.HammershopTeeth, paired=F)


# Germain_2012 ----

iso.Germain <- data[which(data$Reference == "Germain_2012"), ]
iso.Germain <- droplevels(iso.Germain)

table(iso.Germain$System)

# Carbono
ggplot(iso.Germain, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.Germain$d13C[iso.Germain$System == "Wild"])
shapiro.test(iso.Germain$d13C[iso.Germain$System == "Captive"])
var.test(d13C~System, data = iso.Germain)
bartlett.test(d13C~System, data = iso.Germain)

# Amostra desigual, Normal e Homogenea
t.test(d13C~System, data = iso.Germain, var.equal=T, paired=F)
#wilcox.test(d13C~System, data = iso.Germain, paired=F)

# Nitro
ggplot(iso.Germain, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.Germain$d15N[iso.Germain$System == "Wild"])
shapiro.test(iso.Germain$d15N[iso.Germain$System == "Captive"])
var.test(d15N~System, data = iso.Germain)
bartlett.test(d15N~System, data = iso.Germain)

# Amostra desigual, Nao-Normal e Homogenea
#t.test(d15N~System, data = iso.Germain, var.equal=T, paired=F)
wilcox.test(d15N~System, data = iso.Germain, paired=F)


# Jenkins_2020 ----

# Jenkins - F_arctica Plasma ----
iso.JenkinsFA <- data[which(data$Reference == "Jenkins_2020" & 
                              data$Tissue == "Plasma" & data$Specie == "F_arctica"), ]
iso.JenkinsFA <- droplevels(iso.JenkinsFA)

table(iso.JenkinsFA$System)

# Carbo
ggplot(iso.JenkinsFA, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.JenkinsFA$d13C[iso.JenkinsFA$System == "Wild"])
shapiro.test(iso.JenkinsFA$d13C[iso.JenkinsFA$System == "Captive"])
var.test(d13C~System, data = iso.JenkinsFA)
bartlett.test(d13C~System, data = iso.JenkinsFA)

# Amostra desigual, Nao-Normal e Homogenea
#t.test(d13C~System, data = iso.JenkinsFA, var.equal=T, paired=F)
wilcox.test(d13C~System, data = iso.JenkinsFA, paired=F)

# Nitro
ggplot(iso.JenkinsFA, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.JenkinsFA$d15N[iso.JenkinsFA$System == "Wild"])
shapiro.test(iso.JenkinsFA$d15N[iso.JenkinsFA$System == "Captive"])
var.test(d15N~System, data = iso.JenkinsFA)
bartlett.test(d15N~System, data = iso.JenkinsFA)

# Amostra desigual, Nao-Normal e Homogenea
# t.test(d15N~System, data = iso.JenkinsFA, var.equal=T, paired=F)
wilcox.test(d15N~System, data = iso.JenkinsFA, paired=F)

# Jenkins - F_arctica RBC ----
iso.JenkinsFArbc <- data[which(data$Reference == "Jenkins_2020" & 
                                 data$Tissue == "RBC" & data$Specie == "F_arctica"), ]
iso.JenkinsFArbc <- droplevels(iso.JenkinsFArbc)

table(iso.JenkinsFArbc$System)

# Carbo
ggplot(iso.JenkinsFArbc, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.JenkinsFArbc$d13C[iso.JenkinsFArbc$System == "Wild"])
shapiro.test(iso.JenkinsFArbc$d13C[iso.JenkinsFArbc$System == "Captive"])
var.test(d13C~System, data = iso.JenkinsFArbc)
bartlett.test(d13C~System, data = iso.JenkinsFArbc)

# Amostra desigual, Nao-Normal e Homogenea
t.test(d13C~System, data = iso.JenkinsFArbc, var.equal=T, paired=F)
#wilcox.test(d13C~System, data = iso.JenkinsFArbc, paired=F)

# Nitro
ggplot(iso.JenkinsFArbc, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.JenkinsFArbc$d15N[iso.JenkinsFArbc$System == "Wild"])
shapiro.test(iso.JenkinsFArbc$d15N[iso.JenkinsFArbc$System == "Captive"])
var.test(d15N~System, data = iso.JenkinsFArbc)
bartlett.test(d15N~System, data = iso.JenkinsFArbc)

# Amostra desigual, Nao-Normal e Homogenea
t.test(d15N~System, data = iso.JenkinsFArbc, var.equal=T, paired=F)
#wilcox.test(d15N~System, data = iso.JenkinsFArbc, paired=F)


# Jenkins - U_aalge Plasma ----
iso.JenkinsUA <- data[which(data$Reference == "Jenkins_2020" & 
                              data$Tissue == "Plasma" & data$Specie == "U_aalge"), ]
iso.JenkinsUA <- droplevels(iso.JenkinsUA)

table(iso.JenkinsUA$System)

# Carbo
ggplot(iso.JenkinsUA, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.JenkinsUA$d13C[iso.JenkinsUA$System == "Wild"])
shapiro.test(iso.JenkinsUA$d13C[iso.JenkinsUA$System == "Captive"])
var.test(d13C~System, data = iso.JenkinsUA)
bartlett.test(d13C~System, data = iso.JenkinsUA)

# Amostra desigual, Nao-Normal e Homogenea
#t.test(d13C~System, data = iso.JenkinsUA, var.equal=T, paired=F)
wilcox.test(d13C~System, data = iso.JenkinsUA, paired=F)

# Nitro
ggplot(iso.JenkinsUA, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.JenkinsUA$d15N[iso.JenkinsUA$System == "Wild"])
shapiro.test(iso.JenkinsUA$d15N[iso.JenkinsUA$System == "Captive"])
var.test(d15N~System, data = iso.JenkinsUA)
bartlett.test(d15N~System, data = iso.JenkinsUA)

# Amostra desigual, Nao-Normal e Homogenea
# t.test(d15N~System, data = iso.JenkinsUA, var.equal=T, paired=F)
wilcox.test(d15N~System, data = iso.JenkinsUA, paired=F)

# Jenkins - U_aagae RBC ----
iso.JenkinsUArbc <- data[which(data$Reference == "Jenkins_2020" & 
                                 data$Tissue == "RBC" & data$Specie == "U_aalge"), ]
iso.JenkinsUArbc <- droplevels(iso.JenkinsUArbc)

table(iso.JenkinsUArbc$System)

# Carbo
ggplot(iso.JenkinsUArbc, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.JenkinsUArbc$d13C[iso.JenkinsUArbc$System == "Wild"])
shapiro.test(iso.JenkinsUArbc$d13C[iso.JenkinsUArbc$System == "Captive"])
var.test(d13C~System, data = iso.JenkinsUArbc)
bartlett.test(d13C~System, data = iso.JenkinsUArbc)

# Amostra desigual, Nao-Normal e Homogenea
t.test(d13C~System, data = iso.JenkinsUArbc, var.equal=T, paired=F)
#wilcox.test(d13C~System, data = iso.JenkinsUArbc, paired=F)

# Nitro
ggplot(iso.JenkinsUArbc, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.JenkinsUArbc$d15N[iso.JenkinsUArbc$System == "Wild"])
shapiro.test(iso.JenkinsUArbc$d15N[iso.JenkinsUArbc$System == "Captive"])
var.test(d15N~System, data = iso.JenkinsUArbc)
bartlett.test(d15N~System, data = iso.JenkinsUArbc)

# Amostra desigual, Normal e Heterogenea
t.test(d15N~System, data = iso.JenkinsUArbc, var.equal=F, paired=F)
#wilcox.test(d15N~System, data = iso.JenkinsUArbc, paired=F)


# Cardona_2017 ----

iso.Cardona <- data[which(data$Reference == "Cardona_2017"), ]
iso.Cardona <- droplevels(iso.Cardona)

table(iso.Cardona$System)

# Carbono
ggplot(iso.Cardona, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

cardona.null <- lm(d13C ~ 1, iso.Cardona)
cardona.Npsr <- lm(d13C ~ System, iso.Cardona)
cardona.psr <- lmer(d13C ~ System + (1 | Id), iso.Cardona)
summary(cardona.psr)

AIC(cardona.null, cardona.Npsr, cardona.psr)
drop1(update(cardona.psr, REML = F), test = "Chisq")

# Nitro
ggplot(iso.Cardona, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

cardona.null <- lm(d15N ~ 1, iso.Cardona)
cardona.Npsr <- lm(d15N ~ System, iso.Cardona)
cardona.psr <- lmer(d15N ~ System + (1 | Id), iso.Cardona)
summary(cardona.psr)

AIC(cardona.null, cardona.Npsr, cardona.psr)
drop1(update(cardona.psr, REML = T), test = "Chisq")


# Codron_2013 ----
iso.Codron <- data[which(data$Reference == "Codron_2013"), ]
iso.Codron <- droplevels(iso.Codron)

table(iso.Codron$System)

# Carbono
ggplot(iso.Codron, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

codron.null <- lm(d13C ~ 1, iso.Codron)
codron.Npsr <- lm(d13C ~ System, iso.Codron)
codron.psr <- lmer(d13C ~ System + (1 | Id), iso.Codron)
summary(codron.psr)

AIC(codron.null, codron.Npsr, codron.psr)
drop1(update(codron.psr, REML = F), test = "Chisq")

# Nitrogenio
ggplot(iso.Codron, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

codron.null <- lm(d15N ~ 1, iso.Codron)
codron.Npsr <- lm(d15N ~ System, iso.Codron)
codron.psr <- lmer(d15N ~ System + (1 | Id), iso.Codron)
summary(codron.psr)

AIC(codron.null, codron.Npsr, codron.psr)
drop1(update(codron.psr, REML = F), test = "Chisq")

# Navarro_2009 ----

iso.Navarro <- data[which(data$Reference == "Navarro_2009"), ]
iso.Navarro <- droplevels(iso.Navarro)
iso.Navarro$System <- gsub("Wild_Agr", "Wild", iso.Navarro$System)
iso.Navarro$System <- gsub("Wild_Nat", "Wild", iso.Navarro$System)

table( iso.Navarro$Local, iso.Navarro$System, iso.Navarro$Tissue)

iso.NavBlood <- iso.Navarro[which(iso.Navarro$Tissue == "Blood"), ]
iso.NavBlood <- droplevels(iso.NavBlood)

iso.NavClaw <- iso.Navarro[which(iso.Navarro$Tissue == "Claw"), ]
iso.NavClaw <- droplevels(iso.NavClaw)

iso.NavHair <- iso.Navarro[which(iso.Navarro$Tissue == "Hair"), ]
iso.NavHair <- droplevels(iso.NavHair)

iso.NavMuscle <- iso.Navarro[which(iso.Navarro$Tissue == "Muscle"), ]
iso.NavMuscle <- droplevels(iso.NavMuscle)


# Navarro - Blood ----
# Navarro - Carb - Blood
ggplot(iso.NavBlood, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.NavBlood$d13C[iso.NavBlood$System == "Wild"])
shapiro.test(iso.NavBlood$d13C[iso.NavBlood$System == "Captive"])
var.test(d13C~System, data = iso.NavBlood)
bartlett.test(d13C~System, data = iso.NavBlood)

# Amostra desigual, Nao-Normal e Homogenea
# t.test(d13C~System, data = iso.NavBlood, var.equal=T, paired=F)
wilcox.test(d13C~System, data = iso.NavBlood, paired=F)

# Navarro - Nitro - Blood
ggplot(iso.NavBlood, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.NavBlood$d15N[iso.NavBlood$System == "Wild"])
shapiro.test(iso.NavBlood$d15N[iso.NavBlood$System == "Captive"])
var.test(d15N~System, data = iso.NavBlood)
bartlett.test(d15N~System, data = iso.NavBlood)

# Amostra desigual, Nao-Normal e Heterogenea
# t.test(dd15N13C~System, data = iso.NavBlood, var.equal=T, paired=F)
wilcox.test(d15N~System, data = iso.NavBlood, paired=F)

# Navarro - Claw ----
# Navarro - Carb - Claw
ggplot(iso.NavClaw, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.NavClaw$d13C[iso.NavClaw$System == "Wild"])
shapiro.test(iso.NavClaw$d13C[iso.NavClaw$System == "Captive"])
var.test(d13C~System, data = iso.NavClaw)
bartlett.test(d13C~System, data = iso.NavClaw)

# Amostra desigual, Normal e Homogenea
t.test(d13C~System, data = iso.NavClaw, var.equal=T, paired=F)
# wilcox.test(d13C~System, data = iso.NavClaw, paired=F)

# Navarro - Nitro - Claw
ggplot(iso.NavClaw, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.NavClaw$d15N[iso.NavClaw$System == "Wild"])
shapiro.test(iso.NavClaw$d15N[iso.NavClaw$System == "Captive"])
var.test(d15N~System, data = iso.NavClaw)
bartlett.test(d15N~System, data = iso.NavClaw)

# Amostra desigual, Nao-Normal e Heterogenea
# t.test(dd15N13C~System, data = iso.NavClaw, var.equal=T, paired=F)
wilcox.test(d15N~System, data = iso.NavClaw, paired=F)

# Navarro - Hair ----

# Navarro - Carb - Hair
ggplot(iso.NavHair, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.NavHair$d13C[iso.NavHair$System == "Wild"])
shapiro.test(iso.NavHair$d13C[iso.NavHair$System == "Captive"])
var.test(d13C~System, data = iso.NavHair)
bartlett.test(d13C~System, data = iso.NavHair)

# Amostra desigual, Normal e Homogenea
t.test(d13C~System, data = iso.NavHair, var.equal=T, paired=F)
#wilcox.test(d13C~System, data = iso.NavHair, paired=F)

# Navarro - Nitro - Hair 
ggplot(iso.NavHair, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.NavHair$d15N[iso.NavHair$System == "Wild"])
shapiro.test(iso.NavHair$d15N[iso.NavHair$System == "Captive"])
var.test(d15N~System, data = iso.NavHair)
bartlett.test(d15N~System, data = iso.NavHair)

# Amostra desigual, Nao-Normal e Heterogenea
# t.test(dd15N13C~System, data = iso.NavHair, var.equal=T, paired=F)
wilcox.test(d15N~System, data = iso.NavHair, paired=F)


# Navarro - Muscle ----

# Navarro - Carb - Muscle
ggplot(iso.NavMuscle, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.NavMuscle$d13C[iso.NavMuscle$System == "Wild"])
shapiro.test(iso.NavMuscle$d13C[iso.NavMuscle$System == "Captive"])
var.test(d13C~System, data = iso.NavMuscle)
bartlett.test(d13C~System, data = iso.NavMuscle)

# Amostra desigual, Normal e Homogenea
#t.test(d13C~System, data = iso.NavMuscle, var.equal=T, paired=F)
wilcox.test(d13C~System, data = iso.NavMuscle, paired=F)

# Navarro - Nitro - Muscle 
ggplot(iso.NavMuscle, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.NavMuscle$d15N[iso.NavMuscle$System == "Wild"])
shapiro.test(iso.NavMuscle$d15N[iso.NavMuscle$System == "Captive"])
var.test(d15N~System, data = iso.NavMuscle)
bartlett.test(d15N~System, data = iso.NavMuscle)

# Amostra desigual, Nao-Normal e Heterogenea
# t.test(dd15N13C~System, data = iso.NavMuscle, var.equal=T, paired=F)
wilcox.test(d15N~System, data = iso.NavMuscle, paired=F)


# Trembaczowski_2011 ----

iso.Trembaczowski <- data[which(data$Reference == "Trembaczowski_2011"), ]
iso.Trembaczowski <- droplevels(iso.Trembaczowski)

table( iso.Trembaczowski$Local, iso.Trembaczowski$System, iso.Trembaczowski$Specie)

iso.TrembTruta <- iso.Trembaczowski[which(iso.Trembaczowski$Specie == "S_trutta" &
                                            iso.Trembaczowski$System == "Farmed" | 
                                            iso.Trembaczowski$System == "Wild"), ]
iso.TrembTruta <- droplevels(iso.TrembTruta)

# Carbono
ggplot(iso.TrembTruta, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.TrembTruta$d13C[iso.TrembTruta$System == "Wild"])
shapiro.test(iso.TrembTruta$d13C[iso.TrembTruta$System == "Farmed"])
var.test(d13C~System, data = iso.TrembTruta)
bartlett.test(d13C~System, data = iso.TrembTruta)

# Amostra desigual, Normal e Homogenea
# t.test(d13C~System, data = iso.TrembTruta, var.equal=F, paired=F)
wilcox.test(d13C~System, data = iso.TrembTruta, paired=F)

# Enxofre
ggplot(iso.TrembTruta, aes(y=d24S, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^24*S~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

shapiro.test(iso.TrembTruta$d24S[iso.TrembTruta$System == "Wild"])
shapiro.test(iso.TrembTruta$d24S[iso.TrembTruta$System == "Farmed"])
var.test(d24S~System, data = iso.TrembTruta)
bartlett.test(d24S~System, data = iso.TrembTruta)

# Amostra desigual, Nao-Normal e Homogenea
#t.test(d24S~System, data = iso.TrembTruta, var.equal=T, paired=F)
wilcox.test(d24S~System, data = iso.TrembTruta, paired=F)


# Martinez_2015 ----

iso.Martinez <- data[which(data$Reference == "Martinez_2015"), ]
iso.Martinez <- droplevels(iso.Martinez)

table( iso.Martinez$Local, iso.Martinez$System, iso.Martinez$Tissue)

iso.MartinezEsc <- iso.Martinez[which(iso.Martinez$Tissue == "Escama"), ]
iso.MartinezEsc <- droplevels(iso.MartinezEsc)

iso.MartinezSang <- iso.Martinez[which(iso.Martinez$Tissue == "Sangue"), ]
iso.MartinezSang <- droplevels(iso.MartinezSang)

# Martinez - Escama ----
# Carb
ggplot(iso.MartinezEsc, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

mart.esc.aov <- aov(d13C~System, data = iso.MartinezEsc)
summary(mart.esc.aov)
TukeyHSD(mart.esc.aov)

plot(mart.esc.aov)
shapiro.test(mart.esc.aov$residuals)
bartlett.test(d13C~System, data = iso.MartinezEsc)

kruskal.test(d13C~System, data = iso.MartinezEsc)

# Nitro
ggplot(iso.MartinezEsc, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

mart.esc.aov <- aov(d15N~System, data = iso.MartinezEsc)
summary(mart.esc.aov)
TukeyHSD(mart.esc.aov)

plot(mart.esc.aov)
shapiro.test(mart.esc.aov$residuals)
bartlett.test(d15N~System, data = iso.MartinezEsc)


# Martinez - Sangue ----
# Carb
ggplot(iso.MartinezSang, aes(y=d13C, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^13*C~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

mart.Sang.aov <- aov(d13C~System, data = iso.MartinezSang)
summary(mart.Sang.aov)
#TukeyHSD(mart.Sang.aov)

plot(mart.Sang.aov)
shapiro.test(mart.Sang.aov$residuals)
bartlett.test(d13C~System, data = iso.MartinezSang)

# Nitro
ggplot(iso.MartinezSang, aes(y=d15N, x=System, fill=System))+
  geom_boxplot()+
  ylab(expression({delta}^15*N~'(\u2030)'))+
  xlab(label = "")+theme_classic() +
  theme(axis.text.y = element_text(color = "black", size = 14),
        axis.text.x = element_text(color = "black", size=16),
        panel.border = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.ticks.length = unit(6, "pt"),
        legend.position = "none",
        axis.title = element_text(color = "black", size=18),
        axis.line = element_line(colour = "black", linetype = "solid"))

mart.Sang.aov <- aov(d15N~System, data = iso.MartinezSang)
summary(mart.Sang.aov)
#TukeyHSD(mart.Sang.aov)

plot(mart.Sang.aov)
shapiro.test(mart.Sang.aov$residuals)
bartlett.test(d15N~System, data = iso.MartinezSang)

kruskal.test(d15N~System, data = iso.MartinezSang)