# R Script: Beta Diversity Analysis Across Sites
# Prepared for submission to *Ecology Letters*
# Author: Francesca Cerroti
# Date: 15.07.2025

# Load Required Libraries---------

rm(list= ls())

library(ade4)
library(adespatial)
library(adegraphics)
library(spdep)
library(sp)
library(vegan)
library(dplyr)
library(tidyr)
library(ggplot2)
library(gridExtra)
library(MuMIn)
library(FD)
library(BAT)
library(ape)
library(gawdis)
library(openxlsx)
library(corrplot)
library(factoextra)
library(forcats)
library(readr)

setwd("/Beta-diversity in WTH")

# Import Community and Environmental Data------------

comm_data <- read.csv("CommunitiesFinal.csv", row.names = 1, sep = ";")
comm_data[is.na(comm_data)] <- 0

df_list <- split(comm_data, comm_data$location)
df_list <- lapply(df_list, function(df) df %>% select_if(~ any(. != 0)))
df_fr <- df_list$fr
df_br <- df_list$br
df_i  <- df_list$I

# Calculate Community Metrics----------------

species_richness <- rowSums(comm_data > 0)
cat("Mean species richness:", mean(species_richness), "\n")
cat("SD of species richness:", sd(species_richness), "\n")

comm_data <- comm_data %>%
  rowwise() %>%
  mutate(total_abundance = sum(c_across(Aedes_Geniculatus:Chironomidae_A)))

summary_abundance <- comm_data %>%
  group_by(location) %>%
  summarise(mean_abundance = mean(total_abundance),
            sd_abundance = sd(total_abundance))
print(summary_abundance)

# Import and Preprocess Environmental Data-----------

env_data <- read.csv("EnvironmentalVariables.csv", row.names = 1, sep = ";")
env_data <- env_data[, -1]

env_data$Opening_ar <- pi * (env_data$Width_min_cm_1 / 2) * (env_data$Width_max_cm_1 / 2)
env_data <- env_data[, !names(env_data) %in% c("Width_max_cm_1", "Width_min_cm_1", "WTH_origin", "Water_Temp")]

env_list <- split(env_data, env_data$Location)
mapping <- list(N = 0, S = 180, E = 90, W = 270, NE = 45, NW = 315, SE = 135, SW = 225,
                NNE = 22.5, NNW = 337.5, SSE = 157.5, SSW = 202.5, 
                O = 270, SS = 180, SEE = 135, SWS = 202.5,
                NEN = 22.5, NWN = 337.5, NWW = 292.5, SWW = 247.5, SES = 157.5, "0" = 1)

env_list <- lapply(env_list, function(df) {
  df <- df[, !names(df) %in% "Location"]
  df$Orientation <- toupper(df$Orientation)
  df$Orientation <- sapply(df$Orientation, function(x) ifelse(x %in% names(mapping), mapping[[x]], 1))
  df[is.na(df)] <- 0
  df
})

env_fr <- env_list$fr
env_br <- env_list$br
env_i  <- env_list$I

# Log-transform and Standardize Environmental Data-----------

log_standardize <- function(df) {
  df %>% mutate_all(~ log(. + 1)) %>% decostand(method = "standardize")
}

env_f_s <- log_standardize(env_fr)
env_b_s <- log_standardize(env_br)
env_i_s <- log_standardize(env_i)

write.xlsx(data.frame(mean = apply(env_f_s, 2, mean), sd = apply(env_f_s, 2, sd)), "E_mean_Fr.xlsx", rowNames = TRUE)
write.xlsx(data.frame(mean = apply(env_b_s, 2, mean), sd = apply(env_b_s, 2, sd)), "E_mean_br.xlsx", rowNames = TRUE)
write.xlsx(data.frame(mean = apply(env_i_s, 2, mean), sd = apply(env_i_s, 2, sd)), "E_mean_i.xlsx", rowNames = TRUE)

# Remove correlated variables----------

remove_correlated_vars <- function(df, site, extra_remove = NULL, fine_om_vars = NULL) {
  if (!is.null(extra_remove)) {
    df <- df[, !names(df) %in% extra_remove]
  }
  
  cor_mat <- cor(df)
  corrplot::corrplot(cor_mat, method = "number")
  
  cor_thresh <- 0.69
  high_corr <- which(cor_mat > cor_thresh & cor_mat < 1, arr.ind = TRUE)
  
  if (nrow(high_corr) > 0) {
    for (i in 1:nrow(high_corr)) {
      var1 <- rownames(cor_mat)[high_corr[i, 1]]
      var2 <- colnames(cor_mat)[high_corr[i, 2]]
      cat("(", site, ") Highly correlated:", var1, "and", var2, "\n")
    }
  }
  
  if (!is.null(fine_om_vars)) {
    df <- df %>% mutate(sum_fine_organic_matter = rowSums(across(all_of(fine_om_vars))))
    df <- df[, !names(df) %in% fine_om_vars]
  }
  
  return(df)
}

env_f_s <- remove_correlated_vars(
  env_f_s,
  site = "FR",
  extra_remove = c("Turbidity", "CDOM", "Chlorophill.A", "Water_Volume", "Conductivity", "Saturation_."),
  fine_om_vars = c("weight_500", "weight_300")
)

env_b_s <- remove_correlated_vars(
  env_b_s,
  site = "BR",
  extra_remove = c("Ammonium", "Water_Volume", "Saturation_."),
  fine_om_vars = c("weight_500", "weight_300")
)

env_i_s <- remove_correlated_vars(
  env_i_s,
  site = "IN",
  extra_remove = c("Chlorid", "Nitrat", "Sulfat", "Water_Volume", "Opening_ar", "Saturation_.", "Conductivity"),
  fine_om_vars = c("weight_500", "weight_300", "Coarse_weight")
)

# Remove unused community data columns------
df_fr <- df_fr[, -1]
df_br <- df_br[, -1]
df_i  <- df_i[, -1]


# Spatial Data Processing and MEM Calculation - France----------

space <- read.csv("CoordinatesFrance.csv", row.names = 1, sep = ";")

# Convert to Cartesian and center
space <- transform(space,
                   y = y * 40e6 / 360,
                   x = x * cos(y * pi / 180) * 40e6 / 360
)
space$x <- space$x - mean(range(space$x))
space$y <- space$y - mean(range(space$y))

# Build Gabriel graph and weights
nbgab <- graph2nb(gabrielneigh(space), sym = TRUE)
dists <- nbdists(nbgab, space)
fdist <- lapply(dists, function(x) 1 - x / max(dist(space)))
lw <- nb2listw(nbgab, glist = fdist)

# Generate MEM
mem.gab <- mem(lw)
MEM_p <- mem.gab[, 1:16]
MEM_n <- mem.gab[, 17:ncol(mem.gab)]

# Visualize selected MEM axes
plot(mem.gab[, c(1, 15, 35)], SpORcoords = space)

# Beta-Diversity Decomposition: France ------------------

# Calculate taxonomic beta-diversity using Podani family (percentage difference)
out.comp.fr <- beta.div.comp(df_fr, coef = "S", quant = TRUE)

# Extract and transform dissimilarity matrices
richness.fr <- out.comp.fr$rich
replacement.fr <- sqrt(out.comp.fr$repl)
totalbeta.fr <- sqrt(out.comp.fr$D)

# LCBD (Local Contribution to Beta Diversity)
LCBD.fr.rich <- LCBD.comp(richness.fr)
LCBD.fr.repl <- LCBD.comp(replacement.fr, sqrt.D = TRUE)
LCBD.fr.beta <- LCBD.comp(totalbeta.fr, sqrt.D = TRUE)

# Summary of partitioning
out.comp.fr$Note
out.comp.fr$part

# Ordinations (PCoA)
pco.fr.rich <- wcmdscale(richness.fr)
pco.fr.repl <- wcmdscale(replacement.fr)
pco.fr.beta <- wcmdscale(totalbeta.fr)

#Variable Selection and Variance Partitioning: France - taxonomic-richness ------------

# Select environmental variables

set.seed(38832)
mod_env <- ordiR2step(dbrda(pco.fr.rich ~ 1, data = env_f_s),
                      dbrda(pco.fr.rich ~ ., data = env_f_s), trace = FALSE,
                      Pin = 0.9, Pout = 1)

# View selected environmental variables


# Select spatial variables
set.seed(373287732)
mod_mem <- ordiR2step(dbrda(pco.fr.rich ~ 1, data = MEM_p),
                      dbrda(pco.fr.rich ~ ., data = MEM_p), trace = TRUE,
                      Pin = 0.9, Pout = 1)

# MEM_n 
mem_test <- ordiR2step(dbrda(pco.fr.rich ~ 1, data = MEM_n),
                       dbrda(pco.fr.rich ~ ., data = MEM_n), trace = TRUE)


# Final selected variables
env_rich_sig <- env_f_s[, c("sum_fine_organic_matter", "Ammonium")]
mems <- mem.gab[, c(1, 20)]

# Variance partitioning and plot
vrichsig <- varpart(pco.fr.rich, env_rich_sig, mems)
plot(vrichsig, main = "France Richness ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# Significance tests
anova.cca(rda(pco.fr.rich, env_rich_sig))
anova.cca(rda(pco.fr.rich, env_rich_sig, cbind(mems)))
anova.cca(rda(pco.fr.rich, env_rich_sig, mems))
anova.cca(rda(pco.fr.rich, mems))
anova.cca(rda(pco.fr.rich, mems, cbind(env_rich_sig)))

# Variable Selection and Variance Partitioning: France (Replacement) ------------

# Select environmental variables
set.seed(137)
mod_env_repl <- ordiR2step(dbrda(pco.fr.repl ~ 1, data = env_f_s),
                           dbrda(pco.fr.repl ~ ., data = env_f_s), trace = TRUE, Pin = 0.9, Pout = 1)


# Select spatial variables
set.seed(1234)
mod_mem_repl <- ordiR2step(dbrda(pco.fr.repl ~ 1, data = MEM_p),
                           dbrda(pco.fr.repl ~ ., data = MEM_p), trace = FALSE, Pin = 0.9, Pout = 1)


# MEM_n selection
mem_test_repl <- ordiR2step(dbrda(pco.fr.repl ~ 1, data = MEM_n),
                            dbrda(pco.fr.repl ~ ., data = MEM_n), trace = TRUE)


# Final variables
env_rep_sig <- env_f_s[, c("Phosphat", "Chlorid")]
mem_rep_sig <- mem.gab[, c(1, 6, 9, 10, 14)]

# Variance partitioning
vreplsig <- varpart(pco.fr.repl, env_rep_sig, mem_rep_sig)
plot(vreplsig, main = "France Replacement ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.fr.repl, env_rep_sig))
anova.cca(rda(pco.fr.repl, env_rep_sig, cbind(mem_rep_sig)))
anova.cca(rda(pco.fr.repl, env_rep_sig, mem_rep_sig))
anova.cca(rda(pco.fr.repl, mem_rep_sig))
anova.cca(rda(pco.fr.repl, mem_rep_sig, cbind(env_rep_sig)))

# Variable Selection and Variance Partitioning: France (Total Beta) ------------
# Select environmental variables
set.seed(123812892)
mod_env_beta <- ordiR2step(dbrda(pco.fr.beta ~ 1, data = env_f_s),
                           dbrda(pco.fr.beta ~ ., data = env_f_s), trace = FALSE, Pin = 0.9, Pout = 1)

# Select spatial variables
set.seed(127)
mod_mem_beta <- ordiR2step(dbrda(pco.fr.beta ~ 1, data = MEM_p),
                           dbrda(pco.fr.beta ~ ., data = MEM_p), trace = FALSE, Pin = 0.9, Pout = 1)

#MEM_n test
mem_test_beta <- ordiR2step(dbrda(pco.fr.beta ~ 1, data = MEM_n),
                            dbrda(pco.fr.beta ~ ., data = MEM_n), trace = TRUE)

# Final variables
env_beta_sig <- env_f_s[, c("Ph", "Coarse_weight", "Ammonium")]
mem_beta_sig <- mem.gab[, c(1, 6, 14, 16)]

# Variance partitioning
vbetasig <- varpart(pco.fr.beta, env_beta_sig, mem_beta_sig)
plot(vbetasig, main = "France Total Beta ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.fr.beta, env_beta_sig))
anova.cca(rda(pco.fr.beta, env_beta_sig,cbind (mem_beta_sig)))
anova.cca(rda(pco.fr.beta, env_beta_sig,mem_beta_sig))
anova.cca(rda(pco.fr.beta, mem_beta_sig))
anova.cca(rda(pco.fr.beta, mem_beta_sig, cbind(env_beta_sig)))

# Functional Beta-Diversity Decomposition: France-----------

traitf<-read.csv("TraitFrance.csv",row.names=1, sep=";")

traitf <- traitf[!(rownames(traitf) %in% c("Brachypalpoids_Lentus")), ]

traitf <- drop_na(traitf)

#For the traits analyis we can't have communitiy with only one species
species_count <- rowSums(df_fr> 0)

# Identify rows with only one species
rows_with_one_species <- which(species_count == 1)
rows_with_one_species

#WTH=14,29,54,10,21,35 

# Remove those rows from the dataframe
df_fr<- df_fr[-rows_with_one_species, ]
env_f_s<- env_f_s[-rows_with_one_species, ]
MEM_n<- MEM_n[-rows_with_one_species, ]
MEM_p<- MEM_p[-rows_with_one_species, ]
mem.gab<- mem.gab[-rows_with_one_species, ]
space<- space[-rows_with_one_species, ]

### I need to remove traits where the variance is equal to 0
traitf <- dplyr::select(traitf, -AS4, -BS1, -DM2, -FG5, -LO7, -RM5)

# Fuzzy coding
set.seed(32)
gaw_dist <- gawdis::gawdis(traitf, groups=c(rep(1,3),rep(2,4),rep(3,1),rep(4,5),rep(5,6),rep(6,4),rep(7,3)),
                           fuzzy=c(1:7), w.type = "optimized")

#the unbalance distribution of my traits is the nature of my communities, that's fine

set.seed(9)
hvf <- BAT::kernel.build(comm = as.matrix(df_fr), trait = gaw_dist, 
                         abund = TRUE, #samples are abundances
                         method = "box", 
                         axes = 3) #conver

##calculate functional beta diversity
system.time(beta <- BAT::kernel.beta(comm = hvf, func="Sorenson") )

LCBD_Beta_complete = LCBD.comp(beta$Btotal, sqrt.D=TRUE) 
LCBD_Beta_complete$beta

LCBD_Beta_repl = LCBD.comp(beta$Brepl, sqrt.D=TRUE)  
LCBD_Beta_repl$beta 

LCBD_Beta_rich = LCBD.comp(beta$Brich, sqrt.D=TRUE)  
LCBD_Beta_rich$beta 

# Ordinations (PCoA) – Functional Beta Components (France)

pco.fr.rich <- wcmdscale(beta$Brich)
pco.fr.repl <- wcmdscale(beta$Brepl)
pco.fr.beta <- wcmdscale(beta$Btotal)

# Variable Selection – Functional Richness (France)

# Environmental selection
set.seed(137)
mod_env_fr_rich <- ordiR2step(dbrda(pco.fr.rich ~ 1, data = env_f_s),
                              dbrda(pco.fr.rich ~ ., data = env_f_s),
                              trace = TRUE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_p)
set.seed(127823)
mod_mem_fr_rich <- ordiR2step(dbrda(pco.fr.rich ~ 1, data = MEM_p),
                              dbrda(pco.fr.rich ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_n)
set.seed(127823)
mod_memn_fr_rich <- ordiR2step(dbrda(pco.fr.rich ~ 1, data = MEM_n),
                               dbrda(pco.fr.rich ~ ., data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)

# Final selected variables
env_rich_sig <- env_f_s[, c("Chlorid")]
mem_rich_sig <- mem.gab[, c(19, 24, 33)]

# Variance partitioning
vrichsig <- varpart(pco.fr.rich, env_rich_sig, mem_rich_sig)
plot(vrichsig, main = "France Functional Richness ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.fr.rich, env_rich_sig))
anova.cca(rda(pco.fr.rich, env_rich_sig, cbind(mem_rich_sig)))
anova.cca(rda(pco.fr.rich, env_rich_sig, mem_rich_sig))
anova.cca(rda(pco.fr.rich, mem_rich_sig))
anova.cca(rda(pco.fr.rich, mem_rich_sig, cbind(env_rich_sig)))

# Variable Selection – Functional Replacement (France)

# Environmental selection
set.seed(137)
mod_env_fr_repl <- ordiR2step(dbrda(pco.fr.repl ~ 1, data = env_f_s),
                              dbrda(pco.fr.repl ~ ., data = env_f_s),
                              trace = TRUE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_p)
set.seed(48)
mod_mem_fr_repl <- ordiR2step(dbrda(pco.fr.repl ~ 1, data = MEM_p),
                              dbrda(pco.fr.repl ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_n)
set.seed(1234)
mod_memn_fr_repl <- ordiR2step(dbrda(pco.fr.repl ~ 1, data = MEM_n),
                               dbrda(pco.fr.repl ~ ., data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)

# Final variables
env_repl_sig <- env_f_s[, c("Ammonium")]
mem_repl_sig <- mem.gab[, c(1,21)]

# Variance partitioning
vreplsig <- varpart(pco.fr.repl, env_repl_sig, mem_repl_sig)
plot(vreplsig, main = "France Functional Replacement ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests

anova.cca(rda(pco.fr.repl, env_repl_sig))
anova.cca(rda(pco.fr.repl, env_repl_sig, cbind (mem_repl_sig)))
anova.cca(rda(pco.fr.repl, env_repl_sig,mem_repl_sig))
anova.cca(rda(pco.fr.repl, mem_repl_sig))
anova.cca(rda(pco.fr.repl, mem_repl_sig,cbind(env_repl_sig)))

# Variable Selection – Functional Total Beta (France)

# Environmental selection
set.seed(123812892)
mod_env_fr_beta <- ordiR2step(dbrda(pco.fr.beta ~ 1, data = env_f_s),
                              dbrda(pco.fr.beta ~ ., data = env_f_s),
                              trace = FALSE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_p)
set.seed(14435)
mod_mem_fr_beta <- ordiR2step(dbrda(pco.fr.beta ~ 1, data = MEM_p),
                              dbrda(pco.fr.beta ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_n)
set.seed(1)
mod_memn_fr_beta <- ordiR2step(dbrda(pco.fr.beta ~ 1, data = MEM_n),
                               dbrda(pco.fr.beta ~ ., data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)


# Final variables
env_beta_sig <- env_f_s[, c("Chlorid")]
mem_beta_sig <- mem.gab[, c(1, 30)]

# Variance partitioning
vbetasig <- varpart(pco.fr.beta, env_beta_sig, mem_beta_sig)
plot(vbetasig, main = "France Functional Total Beta ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.fr.beta, env_beta_sig))
anova.cca(rda(pco.fr.beta, env_beta_sig, cbind(mem_beta_sig)))
anova.cca(rda(pco.fr.beta, mem_beta_sig, env_beta_sig))
anova.cca(rda(pco.fr.beta, mem_beta_sig))
anova.cca(rda(pco.fr.beta, mem_beta_sig, cbind(env_beta_sig)))

#Spatial Data Processing and MEM Calculation - Brazil-----------
spaceb<- read.csv("CoordinatesBrazil.csv", row.names=1, sep=";")

spaceb<-spaceb[, -c(3)]

# Convert to Cartesian and center
spaceb <- transform(
  spaceb,
  y = lat_S* 40e6 / 360,
  x = long_W * cos(lat_S * pi / 180) * 40e6 / 360
)

spaceb$x <- spaceb$x - mean(range(spaceb$x))
spaceb$y<- spaceb$y - mean(range(spaceb$y))

#remove old columns
spaceb<-spaceb[,-c(1,2)]

# Build Gabriel graph and weights
nbgabmy<-graph2nb(gabrielneigh(spaceb), sym=TRUE)
distgab<- nbdists(nbgabmy, spaceb)
fdist <- lapply(distgab, function(x) 1 - x/max(dist(spaceb)))
listwgab <- nb2listw(nbgabmy, glist = fdist)
listwgab

# Generate MEM
mem.gab_b <- mem(listwgab)
mem.gab_b

MEM_p<-mem.gab_b[,c(1:14)]
MEM_n<-mem.gab_b[,c(15:32)]


# Beta-Diversity Decomposition: Brazil-----------
set.seed(23733)

# Calculate taxonomic beta-diversity using Podani family (percentage difference)
out.comp.br = beta.div.comp(df_br, coef="S", quant=TRUE)

# Extract and transform dissimilarity matrices
richness.br <- out.comp.br$rich
replacement.br <- sqrt(out.comp.br$repl)
totalbeta.br <- sqrt(out.comp.br$D)

# LCBD (Local Contribution to Beta Diversity)
LCBD.br.rich <- LCBD.comp(richness.br)
LCBD.br.repl <- LCBD.comp(replacement.br, sqrt.D = TRUE)
LCBD.br.beta <- LCBD.comp(totalbeta.br, sqrt.D = TRUE)

# Summary of partitioning
out.comp.br$Note
out.comp.br$part

# Ordinations (PCoA)
pco.br.rich <- wcmdscale(richness.br)
pco.br.repl <- wcmdscale(replacement.br)
pco.br.beta <- wcmdscale(totalbeta.br)

#Variable Selection and Variance Partitioning: Brazil- taxonomic-richness ------------

# Select environmental variables

set.seed(38832)
mod_env <- ordiR2step(dbrda(pco.br.rich ~ 1, data = env_b_s),
                      dbrda(pco.br.rich ~ ., data = env_b_s), trace = TRUE,
                      Pin = 0.9, Pout = 1)

# Select spatial variables
set.seed(373287732)
mod_mem <- ordiR2step(dbrda(pco.br.rich ~ 1, data = MEM_p),
                      dbrda(pco.br.rich ~ ., data = MEM_p), trace = TRUE,
                      Pin = 0.9, Pout = 1)

# MEM_n 
mem_test <- ordiR2step(dbrda(pco.br.rich ~ 1, data = MEM_n),
                       dbrda(pco.br.rich ~ ., data = MEM_n), trace = TRUE)
mem_test$anova
# Final selected variables
env_rich_sig <- env_b_s[, c("Coarse_weight")]
mems <- mem.gab_b[, c(17,15,27)]

# Variance partitioning and plot
vrichsig <- varpart(pco.br.rich, env_rich_sig, mems)
plot(vrichsig, main = "Brazil Richness ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# Significance tests
anova.cca(rda(pco.br.rich, env_rich_sig))
anova.cca(rda(pco.br.rich, env_rich_sig, cbind(mems)))
anova.cca(rda(pco.br.rich, env_rich_sig, mems))
anova.cca(rda(pco.br.rich, mems))
anova.cca(rda(pco.br.rich, mems, cbind(env_rich_sig)))

# Variable Selection and Variance Partitioning: Brazil (Replacement) ------------

# Select environmental variables
set.seed(137)
mod_env_repl <- ordiR2step(dbrda(pco.br.repl ~ 1, data = env_b_s),
                           dbrda(pco.br.repl ~ ., data = env_b_s), trace = TRUE, Pin = 0.9, Pout = 1)


# Select spatial variables
set.seed(1)
mod_mem_repl <- ordiR2step(dbrda(pco.br.repl ~ 1, data = MEM_p),
                           dbrda(pco.br.repl ~ ., data = MEM_p), trace = TRUE, Pin = 0.9, Pout = 1)


# MEM_n selection
set.seed(532)
mem_test_repl <- ordiR2step(dbrda(pco.br.repl ~ 1, data = MEM_n),
                            dbrda(pco.br.repl ~ MEM26 + MEM30, data = MEM_n), trace = TRUE)


# Final variables
env_rep_sig <- env_b_s[, c("sum_fine_organic_matter")]
mem_rep_sig <- mem.gab_b[, c(7, 26, 30)]

# Variance partitioning
vreplsig <- varpart(pco.br.repl, env_rep_sig, mem_rep_sig)
plot(vreplsig, main = "Brazil Replacement ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.br.repl, env_rep_sig))
anova.cca(rda(pco.br.repl, env_rep_sig, cbind(mem_rep_sig)))
anova.cca(rda(pco.br.repl, env_rep_sig, mem_rep_sig))
anova.cca(rda(pco.br.repl, mem_rep_sig))
anova.cca(rda(pco.br.repl, mem_rep_sig, cbind(env_rep_sig)))

# Variable Selection and Variance Partitioning: Brazil (Total Beta) ------------
# Select environmental variables
set.seed(123812892)
mod_env_beta <- ordiR2step(dbrda(pco.br.beta ~ 1, data = env_b_s),
                           dbrda(pco.br.beta ~ ., data = env_b_s), trace = TRUE, Pin = 0.9, Pout = 1)

# Select spatial variables
set.seed(1223)
mod_mem_beta <- ordiR2step(dbrda(pco.br.beta ~ 1, data = MEM_p),
                           dbrda(pco.br.beta ~ ., data = MEM_p), trace = TRUE, Pin = 0.9, Pout = 1)


#MEM_n test
set.seed(127823)
mem_test_beta <- ordiR2step(dbrda(pco.br.beta ~ 1, data = MEM_n),
                            dbrda(pco.br.beta ~ ., data = MEM_n), trace = TRUE)


# Final variables
env_beta_sig <- env_b_s[, c( "Coarse_weight")]
mem_beta_sig <- mem.gab_b[, c(17,19)]

# Variance partitioning
vbetasig <- varpart(pco.br.beta, env_beta_sig, mem_beta_sig)
plot(vbetasig, main = "Brazil Total Beta ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.br.beta, env_beta_sig))
anova.cca(rda(pco.br.beta, env_beta_sig,cbind (mem_beta_sig)))
anova.cca(rda(pco.br.beta, env_beta_sig,mem_beta_sig))
anova.cca(rda(pco.br.beta, mem_beta_sig))
anova.cca(rda(pco.br.beta, mem_beta_sig, cbind(env_beta_sig)))

## Functional Beta-Diversity Decomposition: Brazil---------------
traitb<-read.csv("TraitBrazil.csv",row.names=1, sep=";")

df_br <- df_br[, !(names(df_br) %in% c("Ostracoda_a", "Nematoda_br"))]
traitb <- traitb[, !(names(traitb) %in% c("Tadpole_a","Hydrophilidae_larva", "Tropisternus_larva",  "Ephydridae","Culicidae","Tanytarsini" ))]

traitb <- drop_na(traitb)
str(traitb)

# Count the number of non-zero entries (species) in each row
species_count <- rowSums(df_br> 0)
rows_with_one_species <- which(species_count == 1)

df_br<- df_br[-rows_with_one_species, ]
df_br <- df_br[!rownames(df_br) %in% c("113", "114"), ]
env_b_s<- env_b_s[!rownames(env_b_s) %in% c("113", "114"), ]
MEM_n<- MEM_n[!rownames(MEM_n) %in% c("30", "11"), ]
MEM_p<- MEM_p[!rownames(MEM_p) %in% c("30", "11"), ]
mem.gab_b<- mem.gab_b[!rownames(mem.gab_b) %in% c("30", "11"), ]
space<-spaceb[-rows_with_one_species, ]

# fuzzy coding
set.seed(57)
gaw_dist <- gawdis::gawdis(traitb, groups=c(rep(1,4),rep(2,5),rep(3,2),rep(4,6),rep(5,7),rep(6,5),rep(7,3)),
                           fuzzy=c(1:7))

hvb <- BAT::kernel.build(comm = as.matrix(df_br), trait = gaw_dist, 
                         abund = TRUE, 
                         method = "box", 
                         axes = 3) 

beta_b <- BAT::kernel.beta(comm = hvb, func="Sorenson") 

LCBD_Beta_complete = LCBD.comp(beta_b$Btotal, sqrt.D=TRUE) 
LCBD_Beta_complete$beta

LCBD_Beta_repl = LCBD.comp(beta_b$Brepl, sqrt.D=TRUE)  
LCBD_Beta_repl$beta

LCBD_Beta_rich = LCBD.comp(beta_b$Brich, sqrt.D=TRUE)  
LCBD_Beta_rich$beta

# Ordinations (PCoA) – Functional Beta Components (Brazil)

pco.br.rich <- wcmdscale(beta_b$Brich)
pco.br.repl <- wcmdscale(beta_b$Brepl)
pco.br.beta <- wcmdscale(beta_b$Btotal)

# Variable Selection – Functional Richness (Brazil)

# Environmental selection
set.seed(123812892)
mod_env_br_rich <- ordiR2step(dbrda(pco.br.rich ~ 1, data = env_b_s),
                              dbrda(pco.br.rich ~ ., data = env_b_s),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_p)
set.seed(1)
mod_mem_br_rich <- ordiR2step(dbrda(pco.br.rich ~ 1, data = MEM_p),
                              dbrda(pco.br.rich ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_n)
set.seed(1)
mod_memn_br_rich <- ordiR2step(dbrda(pco.br.rich ~ 1, data = MEM_n),
                               dbrda(pco.br.rich ~ MEM16 + MEM19, data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)


# Final selected variables
env_rich_sig <- env_b_s[, c("O2_mgL")]
mem_rich_sig <- mem.gab_b[, c(16,19)]

# Variance partitioning
vrichsig <- varpart(pco.br.rich, env_rich_sig, mem_rich_sig)
plot(vrichsig, main = "Brazil Functional Richness ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.br.rich, env_rich_sig))
anova.cca(rda(pco.br.rich, env_rich_sig, cbind(mem_rich_sig)))
anova.cca(rda(pco.br.rich, env_rich_sig, mem_rich_sig))
anova.cca(rda(pco.br.rich, mem_rich_sig))
anova.cca(rda(pco.br.rich, mem_rich_sig, cbind(env_rich_sig)))

# Variable Selection – Functional Replacement (Brazil)

# Environmental selection
set.seed(137)
mod_env_br_repl <- ordiR2step(dbrda(pco.br.repl ~ 1, data = env_b_s),
                              dbrda(pco.br.repl ~., data = env_b_s),
                              trace = TRUE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_p)
set.seed(675843)
mod_mem_br_repl <- ordiR2step(dbrda(pco.br.repl ~ 1, data = MEM_p),
                              dbrda(pco.br.repl ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_n)
set.seed(1)
mod_memn_br_repl <- ordiR2step(dbrda(pco.br.repl ~ 1, data = MEM_n),
                               dbrda(pco.br.repl ~ ., data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)

# Final variables
env_repl_sig <- env_b_s[, c("Nitrat")]
mem_repl_sig <- mem.gab_b[, c(14, 32)]

# Variance partitioning
vreplsig <- varpart(pco.br.repl, env_repl_sig, mem_repl_sig)
plot(vreplsig, main = "Brazil Functional Replacement ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests

anova.cca(rda(pco.br.repl, env_repl_sig))
anova.cca(rda(pco.br.repl, env_repl_sig, cbind (mem_repl_sig)))
anova.cca(rda(pco.br.repl, env_repl_sig,mem_repl_sig))
anova.cca(rda(pco.br.repl, mem_repl_sig))
anova.cca(rda(pco.br.repl, mem_repl_sig,cbind(env_repl_sig)))

# Variable Selection – Functional Total Beta (Brazil)

# Environmental selection
set.seed(9)
mod_env_br_beta <- ordiR2step(dbrda(pco.br.beta ~ 1, data = env_b_s),
                              dbrda(pco.br.beta ~ ., data = env_b_s),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_p)
set.seed(6758439)
mod_mem_br_beta <- ordiR2step(dbrda(pco.br.beta ~ 1, data = MEM_p),
                              dbrda(pco.br.beta ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_n)
set.seed(1)
mod_memn_br_beta <- ordiR2step(dbrda(pco.br.beta ~ 1, data = MEM_n),
                               dbrda(pco.br.beta ~ ., data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)


# Final variables
env_beta_sig <- env_b_s[, c("O2_mgL")]
mem_beta_sig <- mem.gab_b[, c(16)]

# Variance partitioning
vbetasig <- varpart(pco.br.beta, env_beta_sig, mem_beta_sig)
plot(vbetasig, main = "Brazil Functional Total Beta ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.br.beta, env_beta_sig))
anova.cca(rda(pco.br.beta, env_beta_sig, cbind(mem_beta_sig)))
anova.cca(rda(pco.br.beta, mem_beta_sig, env_beta_sig))
anova.cca(rda(pco.br.beta, mem_beta_sig))
anova.cca(rda(pco.br.beta, mem_beta_sig, cbind(env_beta_sig)))


#Spatial Data Processing and MEM Calculation - India-----------
space<- read.csv("CoordinatesIndia.csv", row.names=1, sep=";")

# Convert to Cartesian and center
space <- transform(
  space,
  y = lat* 40e6 / 360,
  x = long * cos(lat * pi / 180) * 40e6 / 360
)

space$x <- space$x - mean(range(space$x))
space$y<- space$y - mean(range(space$y))

space<-space[,-c(1,2)]

# Build Gabriel graph and weights
nbgabmy<-graph2nb(gabrielneigh(space), sym=TRUE)
nb2listw(nbgabmy)
distgab <- nbdists(nbgabmy, space)
fdist <- lapply(distgab, function(x) 1 - x/max(dist(space)))
listwgab <- nb2listw(nbgabmy, glist = fdist)
listwgab

# Generate MEM
mem.gab_i <- mem(listwgab)
MEM_p<-mem.gab_i[,c(1:15)]
MEM_n<-mem.gab_i[,c(16:34)]

# Beta-Diversity Decomposition: India----------
set.seed(1)
# Calculate taxonomic beta-diversity using Podani family (percentage difference)
out.comp.in= beta.div.comp(df_i, coef="S", quant=TRUE)

# Extract and transform dissimilarity matrices
richness.in <- out.comp.in$rich
replacement.in <- sqrt(out.comp.in$repl)
totalbeta.in <- sqrt(out.comp.in$D)

# LCBD (Local Contribution to Beta Diversity)
LCBD.in.rich <- LCBD.comp(richness.in)
LCBD.in.repl <- LCBD.comp(replacement.in, sqrt.D = TRUE)
LCBD.in.beta <- LCBD.comp(totalbeta.in, sqrt.D = TRUE)

# Summary of partitioning
out.comp.in$Note
out.comp.in$part

# Ordinations (PCoA)
pco.in.rich <- wcmdscale(richness.in)
pco.in.repl <- wcmdscale(replacement.in)
pco.in.beta <- wcmdscale(totalbeta.in)


#Variable Selection and Variance Partitioning: India- taxonomic-richness ------------

# Select environmental variables
set.seed(3882)
mod_env <- ordiR2step(dbrda(pco.in.rich ~ 1, data = env_i_s),
                      dbrda(pco.in.rich ~ ., data = env_i_s), trace = TRUE,
                      Pin = 0.9, Pout = 1)

# Select spatial variables
set.seed(373287732)
mod_mem <- ordiR2step(dbrda(pco.in.rich ~ 1, data = MEM_p),
                      dbrda(pco.in.rich ~ ., data = MEM_p), trace = TRUE,
                      Pin = 0.9, Pout = 1)

# MEM_n 
mem_test <- ordiR2step(dbrda(pco.in.rich ~ 1, data = MEM_n),
                       dbrda(pco.in.rich ~ ., data = MEM_n), trace = TRUE)


# Final selected variables
env_rich_sig <- env_i_s[, c("Turbidity", "Ammonium")]
mems <- mem.gab_i[, c(4, 12)]

# Variance partitioning and plot
vrichsig <- varpart(pco.in.rich, env_rich_sig, mems)
plot(vrichsig, main = "India Richness ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# Significance tests
anova.cca(rda(pco.in.rich, env_rich_sig))
anova.cca(rda(pco.in.rich, env_rich_sig, cbind(mems)))
anova.cca(rda(pco.in.rich, env_rich_sig, mems))
anova.cca(rda(pco.in.rich, mems))
anova.cca(rda(pco.in.rich, mems, cbind(env_rich_sig)))

# Variable Selection and Variance Partitioning: Brazil (Replacement) ------------

# Select environmental variables
set.seed(1)
mod_env_repl <- ordiR2step(dbrda(pco.in.repl ~ 1, data = env_i_s),
                           dbrda(pco.in.repl ~ Orientation +Ph, data = env_i_s), trace = TRUE, Pin = 0.9, Pout = 1)

# Select spatial variables
set.seed(12)
mod_mem_repl <- ordiR2step(dbrda(pco.in.repl ~ 1, data = MEM_p),
                           dbrda(pco.in.repl ~ ., data = MEM_p), trace = TRUE, Pin = 0.9, Pout = 1)

# MEM_n selection
set.seed(9)
mem_test_repl <- ordiR2step(dbrda(pco.in.repl ~ 1, data = MEM_n),
                            dbrda(pco.in.repl ~., data = MEM_n), trace = TRUE)


# Final variables
env_rep_sig <- env_i_s[, c("Ph", "Orientation")]
mem_rep_sig <- mem.gab_i[, c(1,21)]

# Variance partitioning
vreplsig <- varpart(pco.in.repl, env_rep_sig, mem_rep_sig)
plot(vreplsig, main = "India Replacement ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.in.repl, env_rep_sig))
anova.cca(rda(pco.in.repl, env_rep_sig, cbind(mem_rep_sig)))
anova.cca(rda(pco.in.repl, env_rep_sig, mem_rep_sig))
anova.cca(rda(pco.in.repl, mem_rep_sig))
anova.cca(rda(pco.in.repl, mem_rep_sig, cbind(env_rep_sig)))

# Variable Selection and Variance Partitioning: Brazil (Total Beta) ------------
# Select environmental variables
set.seed(1)
mod_env_beta <- ordiR2step(dbrda(pco.in.beta ~ 1, data = env_i_s),
                           dbrda(pco.in.beta ~ . , data = env_i_s), trace = TRUE, Pin = 0.9, Pout = 1)

# Select spatial variables
set.seed(127823)
mod_mem_beta <- ordiR2step(dbrda(pco.in.beta ~ 1, data = MEM_p),
                           dbrda(pco.in.beta ~ ., data = MEM_p), trace = TRUE, Pin = 0.9, Pout = 1)

#MEM_n test
set.seed(1273)
mem_test_beta <- ordiR2step(dbrda(pco.in.beta ~ 1, data = MEM_n),
                            dbrda(pco.in.beta ~ ., data = MEM_n), trace = TRUE)


# Final variables
env_beta_sig <- env_i_s[, c( "Turbidity")]
mem_beta_sig <- mem.gab_i[, c(1,26)]

# Variance partitioning
vbetasig <- varpart(pco.in.beta, env_beta_sig, mem_beta_sig)
plot(vbetasig, main = "India Total Beta ~ Env + MEM", digits = 2,
     Xnames = c("Env", "MEM"), cex = 1.5, bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.in.beta, env_beta_sig))
anova.cca(rda(pco.in.beta, env_beta_sig,cbind (mem_beta_sig)))
anova.cca(rda(pco.in.beta, env_beta_sig,mem_beta_sig))
anova.cca(rda(pco.in.beta, mem_beta_sig))
anova.cca(rda(pco.in.beta, mem_beta_sig, cbind(env_beta_sig)))

### Functional Beta-Diversity Decomposition: India-------------------------------------
traiti<-read.csv("TraitIndia.csv",row.names=1, sep=";")

traiti <- drop_na(traiti)

df_i <- df_i[, !(names(df_i) %in% c("Acaria", "Anisoptera", "Copepod", "Daphnia", "Diptera_A", "Muscidae_A", "Nematoda_I", "location"))]
traiti <- drop_na(traiti)

#### fuzzy coding----------
set.seed(4)
gaw_dist <- gawdis::gawdis(traiti, groups=c(rep(1,4),rep(2,5),rep(3,2),rep(4,6),rep(5,7),rep(6,5),rep(7,3)),
                           fuzzy=c(1:7))

hvi <- BAT::kernel.build(comm = as.matrix(df_i), trait = gaw_dist,
                         abund = TRUE,
                         method = "gaussian", 
                         axes = 3) 

beta_in <- BAT::kernel.beta(comm = hvi, func="Sorenson") 

LCBD_Beta_complete = LCBD.comp(beta_in$Btotal, sqrt.D=TRUE) 
LCBD_Beta_complete$beta

LCBD_Beta_repl = LCBD.comp(beta_in$Brepl, sqrt.D=TRUE)  
LCBD_Beta_repl$beta

LCBD_Beta_rich = LCBD.comp(beta_in$Brich, sqrt.D=TRUE)  
LCBD_Beta_rich$beta

# Ordinations (PCoA) – Functional Beta Components (India)
pco.in.rich <- wcmdscale(beta_in$Brich)
pco.in.repl <- wcmdscale(beta_in$Brepl)
pco.in.beta <- wcmdscale(beta_in$Btotal)

# Variable Selection – Functional Richness (India)

# Environmental selection
set.seed(123812892)
mod_env_in_rich <- ordiR2step(dbrda(pco.in.rich ~ 1, data = env_i_s),
                              dbrda(pco.in.rich ~ ., data = env_i_s),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_p)
set.seed(127823)
mod_mem_in_rich <- ordiR2step(dbrda(pco.in.rich ~ 1, data = MEM_p),
                              dbrda(pco.in.rich ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_n)
mod_memn_in_rich <- ordiR2step(dbrda(pco.in.rich ~ 1, data = MEM_n),
                               dbrda(pco.in.rich ~ ., data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)


# Final selected variables
env_rich_sig <- env_i_s[, c("Turbidity")]
mem_rich_sig <- mem.gab_i[, c(19,21, 24, 33)]

# Variance partitioning
vrichsig <- varpart(pco.in.rich, env_rich_sig, mem_rich_sig)
plot(vrichsig, main = "India Functional Richness ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.in.rich, env_rich_sig))
anova.cca(rda(pco.in.rich, env_rich_sig, cbind(mem_rich_sig)))
anova.cca(rda(pco.in.rich, env_rich_sig, mem_rich_sig))
anova.cca(rda(pco.in.rich, mem_rich_sig))
anova.cca(rda(pco.in.rich, mem_rich_sig, cbind(env_rich_sig)))

# Variable Selection – Functional Replacement (India)

# Environmental selection
set.seed(123812892)
mod_env_in_repl <- ordiR2step(dbrda(pco.in.repl ~ 1, data = env_i_s),
                              dbrda(pco.in.repl ~., data = env_i_s),
                              trace = TRUE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_p)
set.seed(137)
mod_mem_in_repl <- ordiR2step(dbrda(pco.in.repl ~ 1, data = MEM_p),
                              dbrda(pco.in.repl ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)

# Spatial selection (MEM_n)
set.seed(1)
mod_memn_in_repl <- ordiR2step(dbrda(pco.in.repl ~ 1, data = MEM_n),
                               dbrda(pco.in.repl ~ MEM33, data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)

# Final variables
env_repl_sig <- env_i_s[, c("Ph")]
mem_repl_sig <- mem.gab_i[, c(1,10,33)]

# Variance partitioning
vreplsig <- varpart(pco.in.repl, env_repl_sig, mem_repl_sig)
plot(vreplsig, main = "inazil Functional Replacement ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.in.repl, env_repl_sig))
anova.cca(rda(pco.in.repl, env_repl_sig, cbind (mem_repl_sig)))
anova.cca(rda(pco.in.repl, env_repl_sig,mem_repl_sig))
anova.cca(rda(pco.in.repl, mem_repl_sig))
anova.cca(rda(pco.in.repl, mem_repl_sig,cbind(env_repl_sig)))

# Variable Selection – Functional Total Beta (India)

# Environmental selection
set.seed(123812892)
mod_env_in_beta <- ordiR2step(dbrda(pco.in.beta ~ 1, data = env_i_s),
                              dbrda(pco.in.beta ~ ., data = env_i_s),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_p)
set.seed(127823)
mod_mem_in_beta <- ordiR2step(dbrda(pco.in.beta ~ 1, data = MEM_p),
                              dbrda(pco.in.beta ~ ., data = MEM_p),
                              trace = TRUE, Pin = 0.9, Pout = 1)


# Spatial selection (MEM_n)
set.seed(4)
mod_memn_in_beta <- ordiR2step(dbrda(pco.in.beta ~ 1, data = MEM_n),
                               dbrda(pco.in.beta ~ ., data = MEM_n),
                               trace = TRUE, Pin = 0.9, Pout = 1)


# Final variables
env_beta_sig <- env_i_s[, c("Ph")]
mem_beta_sig <- mem.gab_i[, c(1,21)]

# Variance partitioning
vbetasig <- varpart(pco.in.beta, env_beta_sig, mem_beta_sig)
plot(vbetasig, main = "India Functional Total Beta ~ Env + MEM",
     digits = 2, Xnames = c("Env", "MEM"), cex = 1.5,
     bg = c("#5ab4ac", "#d8b365"), alpha = 80)

# ANOVA tests
anova.cca(rda(pco.in.beta, env_beta_sig))
anova.cca(rda(pco.in.beta, env_beta_sig, cbind(mem_beta_sig)))
anova.cca(rda(pco.in.beta, mem_beta_sig, env_beta_sig))
anova.cca(rda(pco.in.beta, mem_beta_sig))
anova.cca(rda(pco.in.beta, mem_beta_sig, cbind(env_beta_sig)))


#Figure 1--------
df<-read.csv("Fig_1_Data.csv", sep=";")

colnames(df) <- c("Value", "Forest", "Component", "Dimension")

df$Type <- recode(df$Dimension,
                  "Taxonomic" = "A) Taxonomic",
                  "Functional" = "B) Functional")

df$Forest <- factor(df$Forest, levels = c("Temperate-Mediterranean", "Neotropical", "Palaeotropical"))
df$Component <- factor(df$Component, levels = c("Richness differences", "Turnover"))  # controls stack order

my_colors <- c("Richness differences" = "#4d9221", "Turnover" = "#c51b7d")

ggplot(df, aes(x = Forest, y = Value, fill = Component)) +
  geom_bar(stat = "identity", position = "stack", color = "white", linewidth = 1.1) +
  facet_grid(. ~ Type) +
  scale_fill_manual(values = my_colors) +
  labs(
    x = "Biogeographical region",
    y = expression("Total Beta"),
    fill = "Component"
  ) +
  theme_classic(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(color = "gray30", fill = NA, linewidth = 0.5),
    strip.text.x = element_text(size = 14, face = "bold", color = "gray30"),
    strip.background = element_rect(fill = "gray80"),
    axis.text.x = element_text(angle = 50, hjust = 1),
    legend.title = element_blank(),
    legend.position = "top",
    legend.justification = c("left", "top"),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.key = element_blank()
  )

ggsave("Fig_1.pdf", width = 7, height = 8)

#Figure 2--------

data <-read.csv("Fig_2_data.csv", sep=";")

# Ensure factors are properly ordered for plotting
data$Location <- factor(data$Location, levels = c("Temperate-Mediterranean", "Neotropical", "Palaeotropical"))
data$Component <- factor(data$Diversity_matrix, levels = c("β-diversity", "Richness differences", "Turnover"))
data$dimension <- factor(data$Diversity_dimension, levels = c("Taxonomic", "Functional"))
data$Type <- factor(data$Type, levels = c("Pure_env", "env_space", "Pure_space"))
str(data$Type)


ggplot(data, aes(x = Location, y = R_adjusted, fill = Type)) +
  geom_bar(stat = "identity", width = 0.7, color = "white", linewidth = 0.8) +
  facet_grid(Component ~ dimension) +
  scale_fill_manual(
    values = c("Pure_env" = "#E4D6A7", "Pure_space" = "#4d9221", "env_space" = "#c51b7d"),
    labels = c("Environmental factors","Environmental*Spatial", "Spatial factors")
  ) +
  theme_minimal(base_size = 13) +
  theme(
    strip.text.x = element_text(size = 13, face = "bold", color = "gray30"),
    strip.text.y = element_text(size = 13, face = "bold", color = "gray30"),
    strip.background = element_rect(fill = "gray80"),
    axis.title.y = element_text(margin = margin(r = 10)),
    axis.text.x = element_text(angle = 50, hjust = 1),
    panel.border = element_rect(colour = "gray30", fill = NA, linewidth = 0.5),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    legend.title = element_blank(),
    legend.position = "top", 
    legend.justification = c("left", "top"),
    legend.direction = "horizontal"  
  ) +
  ylab("adjusted R²") +
  xlab("Biogeographical region")

ggsave("Fig_2.pdf", width = 7, height = 9)