---
title: "Second_workflow"
author: "Bingrui chen"
date: "2025-4-13"
output:
  html_document:
    toc: true
    toc_depth: 4
    toc_float: true
    anchor_sections: true
    code_download: true
    highlight: textmate
    theme: journal
  pdf_document:
    toc: true
    toc_depth: '4'
---

```{r setup, include=FALSE}

knitr::opts_knit$set(
  root.dir = "D:/Doctoral_stage/data/Micr_data/result_finnal/zendo_database"#pathway of your files
)

```

# 1 Load R package
 
```{r}

library(ggh4x)
library(magrittr)
library(tidyverse)
library(ggplot2)
library(RColorBrewer)
library(psych)
library(vegan)
library(rdacca.hp)
library(dplyr)
library(Hmisc)
library(reshape2)
library(rstatix)
library(ggpubr)
library(FactoMineR)
library(factoextra)
library(ggbreak)
library(patchwork)
library(ggdist)
library(corrplot)
library(ggcorrplot2)
library(linkET)
library(lme4)
library(lmerTest)
library(piecewiseSEM)
library(bestNormalize)
library(performance)
library(ggeffects)
library(caret)
library(randomForest)
library(rfPermute)
library(glmm.hp)
library(multifunc)
library(missForest)
library(ggsci)
library(scales)
library(ggpmisc)
library(microeco)
library(car)
library(MuMIn)
library(broom.mixed)

```

# 2 Load data

```{r}

sample <- read.csv("./envir_V1.0.csv", header = TRUE, row.names = 1)

total_data <- sample %>% mutate(N_treat = case_when(
    Group == "CK" ~ 0.1,
    Group == "N50" ~ 0.2,
    Group == "N100" ~ 0.3))

colnames(total_data)[colnames(total_data) == "Tree_RGR"] <- "TRGR"
colnames(total_data)[colnames(total_data) == "litter_production"] <- "TLP"

soil_env_data <- total_data

soil_env_data$Site <- factor(soil_env_data$Site, levels = c("GH","DLS","GNJ","JFL"))
soil_env_data$Group <- factor(soil_env_data$Group, levels = c("CK","N50","N100"))

soil_env_data <- soil_env_data %>%
  group_by(Site) %>%
  mutate(across(where(is.numeric), ~ as.numeric(scale(.x, center = T)))) %>%
  ungroup()

```

```{r}

lmms_diff_num <- function(y_var, data = data) {
  
  model <- lmer(as.formula(paste(y_var, "~", "N_treat", "+ (1|Block)")), 
                data = data, na.action = na.omit)
  
  model_res <- summary(model)
  model_res <- as.data.frame(model_res$coefficients) %>% .[-1, ]
  model_res$'Pr(>|t|)' <- round(model_res$'Pr(>|t|)', 4)
  model_res$sig <- ifelse(model_res$'Pr(>|t|)' >= 0.05, '', 
                          ifelse(model_res$'Pr(>|t|)' > 0.01, '*', 
                                 ifelse(model_res$'Pr(>|t|)' > 0.001, '**', '***')))
 model_res$sig_position <- dplyr::case_when(
  model_res$Estimate < 0 ~  0.03,
  model_res$Estimate > 0 ~ -0.04,
  TRUE                  ~  0)
  normality_p <- check_residuals(model) 
  if (!is.numeric(normality_p)) {
    
    normality_p <- unlist(normality_p)  
  }
  model_res$Normality_P <- round(as.numeric(normality_p), 4)
  
  return(list(model_res = model_res))
}


lmms_diff_class <- function(y_var, data = data) {
  model <- lmer(as.formula(paste(y_var, "~", "Group", "+ (1|Block)")), 
                data = data, na.action = na.omit)
  model_res <- summary(model)
  model_res <- as.data.frame(model_res$coefficients) %>% .[-1, ]
  model_res$'Pr(>|t|)' <- round(model_res$'Pr(>|t|)', 4)
  model_res$sig <- ifelse(model_res$'Pr(>|t|)' >= 0.05, '', 
                          ifelse(model_res$'Pr(>|t|)' > 0.01, '*', 
                                 ifelse(model_res$'Pr(>|t|)' > 0.001, '**', '***')))
 model_res$sig_position <- dplyr::case_when(
  model_res$Estimate < 0 ~  0.05,
  model_res$Estimate > 0 ~ -0.05,
  TRUE                  ~  0)
  normality_p <- check_residuals(model) 
  if (!is.numeric(normality_p)) {
    normality_p <- unlist(normality_p)  
  }
  model_res$Normality_P <- round(as.numeric(normality_p), 4)
  
  return(list(model_res = model_res))
}


```

# 3 Changes in soil properties

## 3.1 Trend

```{r}

env_data_split <- split(soil_env_data, soil_env_data$Site)

env_y_vars <- c("C_P", "C_N", "N_P","pH", "SC", "SN", "SP") 

env_fun_diff_num_GH <- map_dfr(env_y_vars, ~ {
  res <- lmms_diff_num(.x, data = env_data_split$GH)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

env_fun_diff_num_DLS <- map_dfr(env_y_vars, ~ {
  res <- lmms_diff_num(.x, data = env_data_split$DLS)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

env_fun_diff_num_GNJ <- map_dfr(env_y_vars, ~ {
  res <- lmms_diff_num(.x, data = env_data_split$GNJ)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

env_fun_diff_num_JFL <- map_dfr(env_y_vars, ~ {
  res <- lmms_diff_num(.x, data = env_data_split$JFL)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

env_fun_diff_num_GH$Site <- "Boreal forest"
env_fun_diff_num_DLS$Site <- "Temperate forest"
env_fun_diff_num_GNJ$Site <- "Subtropical forest"
env_fun_diff_num_JFL$Site <- "Tropical forest"

env_fun_diff_num_mul <- rbind(env_fun_diff_num_GH, env_fun_diff_num_DLS, 
                              env_fun_diff_num_GNJ, env_fun_diff_num_JFL)

env_fun_diff_num_mul$Site <- factor(env_fun_diff_num_mul$Site, 
                                    levels = c("Boreal forest","Temperate forest",
                                               "Subtropical forest","Tropical forest"))

env_fun_diff_num_mul$y_var <- factor(env_fun_diff_num_mul$y_var, 
                                     levels = c("pH","SC","SN","SP","C_N","C_P","N_P"))

env_fun_diff_num_mul

```

```{r}

soil_var_labels <- data.frame(
  y_var = c("pH", "SC", "SN", "SP", "C_N", "C_P", "N_P"),
  label = c("(A)", "(B)", "(C)", "(D)", "(E)", "(F)", "(G)"))

soil_var_labels$y_var <- factor(soil_var_labels$y_var, levels = levels(env_fun_diff_num_mul$y_var))

env_fun_diff_estimate_plot <- ggplot(data = env_fun_diff_num_mul, 
                                     aes(x = Site, y = Estimate, fill = Site)) +
    geom_text(data = soil_var_labels, aes(x = -Inf, y = Inf, label = label), 
            inherit.aes = FALSE, hjust = -0.3, vjust = 1.5, size = 4, 
            fontface = "bold", family = "serif") +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7) + 
  geom_errorbar(
    aes(ymin = Estimate - `Std. Error`, 
        ymax = Estimate + `Std. Error`,
        color = Site), 
    position = position_dodge(width = 0.8),  
    width = 0.1, 
    lwd = 0.4) +
  geom_hline(yintercept = 0, color = "Black", linetype = "dashed", lwd = 0.3) +
  scale_fill_manual(values = c("Boreal forest" = "#08306B", "Temperate forest" = "#2171B5", 
                               "Subtropical forest" = "#6BAED6", "Tropical forest" = "#C6DBEF")) +
  scale_color_manual(values = c("Boreal forest" = "#08306B", "Temperate forest" = "#2171B5", 
                                "Subtropical forest" = "#6BAED6", "Tropical forest" = "#C6DBEF")) +
  facet_manual(vars(y_var), 
               design = "ABCD
                         #EFG", 
               scales = "fixed",
               axes = "all",
               remove_labels = "all",
               labeller = labeller(y_var = c(
                 pH  = "Soil pH", SC = "Soil C content", SN = "Soil N content", SP = "Soil P content",
                 C_N = "Soil C:N ratios", C_P = "Soil C:P ratios", N_P = "Soil N:P ratios"
               ))) +
  scale_x_discrete(limits = c("Boreal forest","Temperate forest",
                              "Subtropical forest","Tropical forest"),
                   labels = c(" "," "," "," "), 
                   position = "bottom") +
  theme_classic() + 
  geom_text(data = env_fun_diff_num_mul, aes(x = Site, y = sig_position, label = sig),
            position = position_dodge(width = 0.8),  size = 4, color = "black") +
  labs(title = NULL, x = NULL, y = "Effect size") +
  #guides(color = guide_legend(title = NULL))+
  theme(legend.position = "inside",
        legend.position.inside = c(0.11, 0.25), # 
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 17),
        axis.text.x = element_text(size = 15, color = "black",angle = 30, hjust = 1),
        axis.text.y = element_text(size = 15, color = "black"),
        axis.ticks = element_line(color = "black", size = 0.3),
        #axis.ticks = element_blank(),
        axis.ticks.length.x = unit(-0.10, "cm"), # 
        strip.background.y = element_rect(colour = NA), 
        axis.line = element_blank(),
        axis.title.x = element_text(size = 16, color = "black"),
        axis.title.y = element_text(size = 16, color = "black"),
        text = element_text(family = "serif"),
        strip.text = element_text(color = "black", size = 17),
        strip.background = element_rect(colour = NA, fill = "#e9ecef", size = 0.5),
        panel.border = element_rect(color = "black", fill = NA, size = 0.5)) 
  
env_fun_diff_estimate_plot

#ggsave("./Fig_2.pdf", env_fun_diff_estimate_plot, width = 10, height = 6)

```

## 3.2 Absolute variation

```{r}

env_data_clas_split <- split(sample, sample$Site)

env_y_vars <- c("C_P", "C_N", "N_P","pH", "SC", "SN", "SP") 

env_fun_diff_clas_GH <- map_dfr(env_y_vars, ~ {
  res <- lmms_diff_class(.x, data = env_data_clas_split$GH)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

env_fun_diff_clas_DLS <- map_dfr(env_y_vars, ~ {
  res <- lmms_diff_class(.x, data = env_data_clas_split$DLS)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

env_fun_diff_clas_GNJ <- map_dfr(env_y_vars, ~ {
  res <- lmms_diff_class(.x, data = env_data_clas_split$GNJ)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

env_fun_diff_clas_JFL <- map_dfr(env_y_vars, ~ {
  res <- lmms_diff_class(.x, data = env_data_clas_split$JFL)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

env_fun_diff_clas_GH$Site <- "Boreal forest"
env_fun_diff_clas_DLS$Site <- "Temperate forest"
env_fun_diff_clas_GNJ$Site <- "Subtropical forest"
env_fun_diff_clas_JFL$Site <- "Tropical forest"

env_fun_diff_clas_mul <- rbind(env_fun_diff_clas_GH, env_fun_diff_clas_DLS, 
                              env_fun_diff_clas_GNJ, env_fun_diff_clas_JFL)

env_fun_diff_clas_mul$Site <- factor(env_fun_diff_clas_mul$Site, 
                                    levels = c("Boreal forest","Temperate forest",
                                               "Subtropical forest","Tropical forest"))

env_fun_diff_clas_mul$N_treat[env_fun_diff_clas_mul$N_treat == "GroupN50"] <- "N50"
env_fun_diff_clas_mul$N_treat[env_fun_diff_clas_mul$N_treat == "GroupN100"] <- "N100"
env_fun_diff_clas_mul$N_treat <- factor(env_fun_diff_clas_mul$N_treat, levels = c("N50","N100"))

env_fun_diff_clas_mul$y_var <- factor(env_fun_diff_clas_mul$y_var, 
                                     levels = c("pH","SC","SN","SP","C_N","C_P","N_P"))

env_fun_diff_clas_mul

```

```{r}

soil_env_est_plot <- ggplot(data = env_fun_diff_clas_mul, 
                            aes(x = Site, y = Estimate, fill = N_treat)) +
    geom_text(data = soil_var_labels, aes(x = -Inf, y = Inf, label = label), 
              inherit.aes = FALSE, hjust = -0.3, vjust = 1.5, size = 3.5, 
              fontface = "bold", family = "serif") +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7) + 
  geom_errorbar(
    aes(ymin = Estimate - `Std. Error`, 
        ymax = Estimate + `Std. Error`,
        color = N_treat), 
    position = position_dodge(width = 0.8),  
    width = 0.1, 
    lwd = 0.4) +
  geom_hline(yintercept = 0, color = "Black", linetype = "dashed", lwd = 0.3) +
  scale_fill_manual(values = c("N50" = "#08306B", "N100" = "#6BAED6")) +
  scale_color_manual(values = c("N50" = "#08306B", "N100" = "#6BAED6")) +
  facet_manual(vars(y_var), 
               design = "ABCD
                         #EFG", 
               scales = "free_y",
               axes = "all",
               remove_labels = "all",
               labeller = labeller(y_var = c(
                 pH  = "Soil pH", SC = "Soil C content", SN = "Soil N content", SP = "Soil P content",
                 C_N = "Soil C:N ratios", C_P = "Soil C:P ratios", N_P = "Soil N:P ratios"
               ))) +
  scale_x_discrete(limits = c("Boreal forest","Temperate forest",
                              "Subtropical forest","Tropical forest"),
                   labels = c("Boreal","Temperate","Subtropical","Tropical"), 
                   position = "bottom") +
  theme_classic() + 
  geom_text(data = env_fun_diff_clas_mul, aes(x = Site, y = sig_position, label = sig),
            position = position_dodge(width = 0.8), size = 4, color = "black") +
  labs(title = NULL, x = NULL, y = "Effect size") +
  theme(legend.position = "inside",
        legend.position.inside = c(0.11, 0.25), # 
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 17),
        axis.text.x = element_text(size = 12, color = "black",angle = 30, hjust = 1),
        axis.text.y = element_text(size = 12, color = "black"),
        axis.ticks = element_line(color = "black", size = 0.3),
        #axis.ticks = element_blank(),
        axis.ticks.length.x = unit(-0.10, "cm"), # 
        strip.background.y = element_rect(colour = NA), 
        axis.line = element_blank(),
        axis.title.x = element_text(size = 14, color = "black"),
        axis.title.y = element_text(size = 14, color = "black"),
        text = element_text(family = "serif"),
        strip.text = element_text(color = "black", size = 14),
        strip.background = element_rect(colour = NA, fill = "#e9ecef", size = 0.5),
        panel.border = element_rect(color = "black", fill = NA, size = 0.5))

soil_env_est_plot

#ggsave("./Fig_S1.pdf", plot = soil_env_est_plot, width = 10, height = 6)

```

# 4 Diversity of microbial communities

```{r}

#Import OTU data
b_otu <- read.csv("./b_OTU_V1.0.csv", header = TRUE, row.names = 1)
f_otu <- read.csv("./f_OTU_V1.0.csv", header = TRUE, row.names = 1)

#Import environmental data
sample <- read.csv("./envir_V1.0.csv", header = TRUE, row.names = 1)
sample$Group <- factor(sample$Group, levels = c("CK", "N50", "N100"))
sample$Site <- factor(sample$Site, levels = c("GH", "DLS", "GNJ", "JFL"))

#Import species annotation data
b_tax <- read.csv("./b_tax_V1.0.csv", header = TRUE, row.names = 1)
f_tax <- read.csv("./f_tax_V1.0.csv", header = TRUE, row.names = 1)

#Import Fasta data
b_fasta <- seqinr::read.fasta("./b_fasta_V1.0.fasta")

```

```{r}

#Build a microtable and store it in a microtable object.
#bacteria
b_dataset <- microtable$new(sample_table = sample, 
                            otu_table = b_otu, 
                            tax_table = b_tax,
                            rep_fasta = b_fasta)
#fungi
f_dataset <- microtable$new(sample_table = sample, 
                            otu_table = f_otu, 
                            tax_table = f_tax)

b_dataset$tidy_dataset()
f_dataset$tidy_dataset()

b_dataset
f_dataset

```

```{r}

sites <- c("GH", "DLS", "GNJ", "JFL")
groups <- c("CK", "N50", "N100")

b_result_list <- list()
f_result_list <- list()


for (site in sites) {
  for (group in groups) {

    dataset_name <- paste0(tolower(site), "_", tolower(group), "_b_dataset")
    
    temp_dataset <- clone(b_dataset)
    temp_dataset$sample_table %<>% subset(Site == site & Group == group)
    temp_dataset$tidy_dataset()
    temp_dataset$filter_taxa(rel_abund = 0, freq = 0.75)
    
    b_result_list[[dataset_name]] <- temp_dataset
    
    message(paste("Processing completed：Site =", site, "| Group =", group))
  }
}

for (site in sites) {
  for (group in groups) {

    dataset_name <- paste0(tolower(site), "_", tolower(group), "_f_dataset")
    
    temp_dataset <- clone(f_dataset)
    temp_dataset$sample_table %<>% subset(Site == site & Group == group)
    temp_dataset$tidy_dataset()
    temp_dataset$filter_taxa(rel_abund = 0, freq = 0.5)
    
    f_result_list[[dataset_name]] <- temp_dataset
    
    message(paste("Processing completed：Site =", site, "| Group =", group))
  }
}

```

```{r}

merge_otu_tables <- function(
    otu_list,       
    fill_na = 0     
) {
  require(tidyverse)
  
  merged_long <- map_df(otu_list, ~ {
    .x %>%
      as.data.frame() %>%
      tibble::rownames_to_column("OTU_ID") %>%
      pivot_longer(cols = -OTU_ID, 
                   names_to = "SampleID", 
                   values_to = "Abundance")
  })
  
  merged_wide <- merged_long %>%
    pivot_wider(
      id_cols = OTU_ID,
      names_from = SampleID,
      values_from = Abundance,
      values_fill = list(Abundance = fill_na)
    ) %>%
    as.data.frame() %>%
    tibble::column_to_rownames("OTU_ID")
  
  return(merged_wide)
}

b_otu_list <- list(
  otu1 = b_result_list$gh_ck_b_dataset$otu_table,
  otu2 = b_result_list$gh_n50_b_dataset$otu_table,
  otu3 = b_result_list$gh_n100_b_dataset$otu_table,
  otu4 = b_result_list$dls_ck_b_dataset$otu_table,
  otu5 = b_result_list$dls_n50_b_dataset$otu_table,
  otu6 = b_result_list$dls_n100_b_dataset$otu_table,
  otu7 = b_result_list$gnj_ck_b_dataset$otu_table,
  otu8 = b_result_list$gnj_n50_b_dataset$otu_table,
  otu9 = b_result_list$gnj_n100_b_dataset$otu_table,
  otu10 = b_result_list$jfl_ck_b_dataset$otu_table,
  otu11 = b_result_list$jfl_n50_b_dataset$otu_table,
  otu12 = b_result_list$jfl_n100_b_dataset$otu_table
)

f_otu_list <- list(
  otu1 = f_result_list$gh_ck_f_dataset$otu_table,
  otu2 = f_result_list$gh_n50_f_dataset$otu_table,
  otu3 = f_result_list$gh_n100_f_dataset$otu_table,
  otu4 = f_result_list$dls_ck_f_dataset$otu_table,
  otu5 = f_result_list$dls_n50_f_dataset$otu_table,
  otu6 = f_result_list$dls_n100_f_dataset$otu_table,
  otu7 = f_result_list$gnj_ck_f_dataset$otu_table,
  otu8 = f_result_list$gnj_n50_f_dataset$otu_table,
  otu9 = f_result_list$gnj_n100_f_dataset$otu_table,
  otu10 = f_result_list$jfl_ck_f_dataset$otu_table,
  otu11 = f_result_list$jfl_n50_f_dataset$otu_table,
  otu12 = f_result_list$jfl_n100_f_dataset$otu_table
)

b_otu_merged <- merge_otu_tables(b_otu_list)
f_otu_merged <- merge_otu_tables(f_otu_list)

```

```{r}

#bacteria
b_dataset_filter <- microtable$new(sample_table = sample, 
                            otu_table = b_otu_merged, 
                            tax_table = b_tax,
                            rep_fasta = b_fasta)
#fungi
f_dataset_filter <- microtable$new(sample_table = sample, 
                            otu_table = f_otu_merged, 
                            tax_table = f_tax)

b_dataset_filter$tidy_dataset()
f_dataset_filter$tidy_dataset()

b_dataset_filter
f_dataset_filter

```

## 4.1 Community diversity

```{r}

#bacteria
b_alpha <- trans_alpha$new(dataset = b_dataset_filter, group = "Group", by_group = "Site")
#fungi
f_alpha <- trans_alpha$new(dataset = f_dataset_filter, group = "Group", by_group = "Site")

```

```{r}

#分析结果整理与保存
names(b_alpha$data_alpha)
names(f_alpha$data_alpha)

b_alpha_result <- b_alpha$data_alpha %>% select(SampleID = SampleID, 
                                        Measure = Measure, 
                                        Value = Value, 
                                        Group = Group, 
                                        Site = Site)

f_alpha_result <- f_alpha$data_alpha %>% select(SampleID = SampleID, 
                                        Measure = Measure, 
                                        Value = Value, 
                                        Group = Group, 
                                        Site = Site)

b_alpha_result$Kindom <- "bacteria"
f_alpha_result$Kindom <- "fungi"

b_alpha_result %<>% pivot_wider(names_from = Measure, values_from = Value)
f_alpha_result %<>% pivot_wider(names_from = Measure, values_from = Value)

alpha_b_pca_split <- split(b_alpha_result, b_alpha_result$Site)
alpha_f_pca_split <- split(f_alpha_result, f_alpha_result$Site)

```

```{r}

names(b_alpha_result)[names(b_alpha_result) == "Chao1"] <- "Chao1_b"
names(f_alpha_result)[names(f_alpha_result) == "Chao1"] <- "Chao1_f"

names(b_alpha_result)[names(b_alpha_result) == "Simpson"] <- "Simpson_b"
names(f_alpha_result)[names(f_alpha_result) == "Simpson"] <- "Simpson_f"

names(b_alpha_result)[names(b_alpha_result) == "Pielou"] <- "Pielou_b"
names(f_alpha_result)[names(f_alpha_result) == "Pielou"] <- "Pielou_f"

total_data <- merge(total_data, b_alpha_result[, c("SampleID", "Chao1_b","Simpson_b","Pielou_b")], 
                    by = "SampleID")

total_data <- merge(total_data, f_alpha_result[, c("SampleID", "Chao1_f","Simpson_f","Pielou_f")], 
                    by = "SampleID")

```

## 4.2 Network complexity

```{r}

np_b_GH <- read.csv("./GH_b_network_properties_SE.csv", header = TRUE, row.names = 1)
np_b_DLS <- read.csv("./DLS_b_network_properties_SE.csv", header = TRUE, row.names = 1)
np_b_GNJ <- read.csv("./GNJ_b_network_properties_SE.csv", header = TRUE, row.names = 1)
np_b_JFL <- read.csv("./JFL_b_network_properties_SE.csv", header = TRUE, row.names = 1)
np_b <- bind_rows(np_b_GH, np_b_DLS, np_b_GNJ, np_b_JFL)

np_f_GH <- read.csv("./GH_f_network_properties_SE.csv", header = TRUE, row.names = 1)
np_f_DLS <- read.csv("./DLS_f_network_properties_SE.csv", header = TRUE, row.names = 1)
np_f_GNJ <- read.csv("./GNJ_f_network_properties_SE.csv", header = TRUE, row.names = 1)
np_f_JFL <- read.csv("./JFL_f_network_properties_SE.csv", header = TRUE, row.names = 1)
np_f <- bind_rows(np_f_GH, np_f_DLS, np_f_GNJ, np_f_JFL)

names(np_b)[names(np_b) == "Density"] <- "Connectance"
names(np_f)[names(np_f) == "Density"] <- "Connectance"

np_b <- merge(np_b, total_data[, c("SampleID", "Block","Site","N_treat")], by = "SampleID")
np_f <- merge(np_f, total_data[, c("SampleID", "Block","Site","N_treat")], by = "SampleID")

np_b$Inv_average_path_length <- -1 * np_b$Average_path_length
np_f$Inv_average_path_length <- -1 * np_f$Average_path_length


```


```{r}

total_data <- merge(total_data, 
                    np_b[, c("SampleID","Vertex","Edge","Inv_average_path_length",
                                 "Average_degree","Average_path_length",
                                 "Modularity","Connectance")], 
                    by = "SampleID")

names(total_data)[names(total_data) == "Vertex"] <- "V_b"
names(total_data)[names(total_data) == "Edge"] <- "E_b"
names(total_data)[names(total_data) == "Average_degree"] <- "AD_b"
names(total_data)[names(total_data) == "Average_path_length"] <- "APL_b"
names(total_data)[names(total_data) == "Inv_average_path_length"] <- "IAPL_b"
names(total_data)[names(total_data) == "Modularity"] <- "Mod_b"
names(total_data)[names(total_data) == "Connectance"] <- "Con_b"

```

```{r}

total_data <- merge(total_data, 
                    np_f[, c("SampleID","Vertex","Edge","Inv_average_path_length",
                                 "Average_degree","Average_path_length",
                                 "Modularity","Connectance")], 
                    by = "SampleID")

names(total_data)[names(total_data) == "Vertex"] <- "V_f"
names(total_data)[names(total_data) == "Edge"] <- "E_f"
names(total_data)[names(total_data) == "Average_degree"] <- "AD_f"
names(total_data)[names(total_data) == "Average_path_length"] <- "APL_f"
names(total_data)[names(total_data) == "Inv_average_path_length"] <- "IAPL_f"
names(total_data)[names(total_data) == "Modularity"] <- "Mod_f"
names(total_data)[names(total_data) == "Connectance"] <- "Con_f"

```


```{r}

ad_nc_data <- total_data %>%
  group_by(Site) %>%
  mutate(across(where(is.numeric), ~ as.numeric(scale(.x, center = T)))) %>%
  ungroup()

colnames(ad_nc_data)

```

```{r}

ad_nc_data_split <- split(ad_nc_data, ad_nc_data$Site)

mic_y_vars <- c("Chao1_b","Simpson_b","Pielou_b",
                "Chao1_f","Simpson_f","Pielou_f",
                "V_b","E_b","AD_b", "IAPL_b","Con_b",
                "V_f","E_f","AD_f", "IAPL_f","Con_f"
                )

ad_nc_diff_GH <- map_dfr(mic_y_vars, ~ {
  res <- lmms_diff_num(.x, data = ad_nc_data_split$GH)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

ad_nc_diff_DLS <- map_dfr(mic_y_vars, ~ {
  res <- lmms_diff_num(.x, data = ad_nc_data_split$DLS)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

ad_nc_diff_GNJ <- map_dfr(mic_y_vars, ~ {
  res <- lmms_diff_num(.x, data = ad_nc_data_split$GNJ)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

ad_nc_diff_JFL <- map_dfr(mic_y_vars, ~ {
  res <- lmms_diff_num(.x, data = ad_nc_data_split$JFL)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

ad_nc_diff_GH$Site <- "Boreal forest"
ad_nc_diff_DLS$Site <- "Temperate forest"
ad_nc_diff_GNJ$Site <- "Subtropical forest"
ad_nc_diff_JFL$Site <- "Tropical forest"

ad_nc_diff_res <- rbind(ad_nc_diff_GH, ad_nc_diff_DLS, ad_nc_diff_GNJ, ad_nc_diff_JFL)

ad_nc_diff_res$Site <- factor(ad_nc_diff_res$Site, levels = c("Boreal forest","Temperate forest",
                                                              "Subtropical forest","Tropical forest"))

ad_nc_diff_res <- ad_nc_diff_res %>%
    mutate(Tax = case_when(
      grepl("_b", ad_nc_diff_res$y_var) ~ "Bacteria",
      grepl("_f", ad_nc_diff_res$y_var) ~ "Fungi",
      TRUE ~ NA_character_
    ))

ad_nc_diff_res$y_var_clear <- gsub("_b|_f", "", ad_nc_diff_res$y_var)
ad_nc_diff_res$y_var_clear <- factor(ad_nc_diff_res$y_var_clear, 
                                    levels = c("Chao1","Simpson","Pielou",
                                               "V","E","AD", "IAPL","Con","Mod"))

ad_nc_diff_res <- ad_nc_diff_res %>%
  mutate(Type = if_else(y_var_clear %in% c("Chao1", "Simpson", "Pielou"), 
                        "alpha_diversity", 
                        "network_complexity"))

ad_nc_diff_res$facet <- paste(ad_nc_diff_res$Tax, ad_nc_diff_res$Type, sep = "_")

#ad_nc_diff_res$facet <- factor(ad_nc_diff_res$facet, 
                               #levels = c("Bacteria_GH","Bacteria_DLS","Bacteria_GNJ","Bacteria_JFL",
                                          #"Fungi_GH","Fungi_DLS","Fungi_GNJ","Fungi_JFL"
                                          #))

ad_nc_diff_res

```


```{r}

#ad_nc_diff_res_split <- split(ad_nc_diff_res, ad_nc_diff_res$Type)

nc_var_labels <- data.frame(
  facet = c("Bacteria_alpha_diversity","Bacteria_network_complexity", 
           "Fungi_alpha_diversity","Fungi_network_complexity"),
  label = c("(A)","(B)","(C)","(D)"))

#ad_nc_diff_res$facet <- factor(ad_nc_diff_res$facet, 
                                    #levels = c("Bacteria_network_complexity", "Bacteria_alpha_diversity",
                                               #"Fungi_network_complexity", "Fungi_alpha_diversity"))

nc_diff_estimate_plot <- ggplot(data = ad_nc_diff_res, 
                                     aes(x = y_var_clear, y = Estimate, fill = Site)) +
  geom_text(data = nc_var_labels, aes(x = -Inf, y = Inf, label = label), 
            inherit.aes = FALSE, hjust = -0.3, vjust = 1.5, size = 3.5, 
            fontface = "bold", family = "serif") +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7) + 
  geom_errorbar(
    aes(ymin = Estimate - `Std. Error`, 
        ymax = Estimate + `Std. Error`,
        color = Site), 
    position = position_dodge(width = 0.8),  
    width = 0.1, 
    lwd = 0.4) +
  geom_hline(yintercept = 0, color = "Black", linetype = "dashed", lwd = 0.3) +
  scale_fill_manual(values = c("Boreal forest" = "#08306B", "Temperate forest" = "#2171B5", 
                               "Subtropical forest" = "#6BAED6", "Tropical forest" = "#C6DBEF")) +
  scale_color_manual(values = c("Boreal forest" = "#08306B", "Temperate forest" = "#2171B5", 
                                "Subtropical forest" = "#6BAED6", "Tropical forest" = "#C6DBEF")) +
  facet_manual(vars(facet), 
               design = "AB
                         CD", 
               scales = "free",
               axes = "all",
               widths  = c(1, 1.8, 1, 1.8),
               remove_labels = "all",
               labeller = labeller(facet = c(
                 Bacteria_network_complexity  = "Bacterial network complexity", 
                 Bacteria_alpha_diversity = "Bacterial alpha diversity", 
                 Fungi_network_complexity = "Fungal network complexity",
                 Fungi_alpha_diversity = "Fungal alpha diversity"))
               ) +
  scale_x_discrete(labels = c("V" = "Vertex", 
                              "E" = "Edge",
                              "AD" = "Average degree",
                              "IAPL" = "Inv APL",
                                "Con" = "Connectance"), 
                   position = "bottom") +
  theme_classic() + 
  geom_text(data = ad_nc_diff_res, aes(x = y_var_clear, y = sig_position, label = sig),
            position = position_dodge(width = 0.8), size = 3, color = "black") +
  labs(title = NULL, x = NULL, y = "Effect size") +
  theme(legend.position = "bottom",
        #legend.position.inside = c(0.11, 0.25), 
        legend.text = element_text(size = 10),
        legend.title = element_blank(),
        legend.key.size = unit(0.4, "cm"),
        axis.text.x = element_text(size = 15, color = "black"),
        axis.text.y = element_text(size = 15, color = "black"),
        axis.ticks = element_line(color = "black", size = 0.3),
        #axis.ticks.length = unit(0.15, "cm"), 
        strip.background.y = element_rect(colour = NA), 
        axis.line = element_blank(),
        axis.title.x = element_text(size = 16, color = "black"),
        axis.title.y = element_text(size = 16, color = "black"),
        text = element_text(family = "serif"),
        strip.text = element_text(color = "black", size = 17),
        strip.background = element_rect(colour = NA, fill = "#e9ecef", size = 0.5),
        panel.border = element_rect(color = "black", fill = NA, size = 0.5))

print(nc_diff_estimate_plot)

#ggsave("./Fig_3.pdf", plot = nc_diff_estimate_plot, width = 10, height = 7)

```

# 5 Response of ecosystem functions to N fertilization

```{r}

emf_data <- total_data %>%
  group_by(Site) %>%
  mutate(across(
    .cols = c(TRGR,TLP,N_treat), 
    .fns = ~ as.numeric(scale(.x, center = TRUE))
  )) %>%
  ungroup()


```


```{r}

emf_data_split <- split(emf_data,emf_data$Site)

emf_y_vars <- c("TRGR","TLP") 

emf_diff_num_GH <- map_dfr(emf_y_vars, ~ {
  res <- lmms_diff_num(.x, data = emf_data_split$GH)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

emf_diff_num_DLS <- map_dfr(emf_y_vars, ~ {
  res <- lmms_diff_num(.x, data = emf_data_split$DLS)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

emf_diff_num_GNJ <- map_dfr(emf_y_vars, ~ {
  res <- lmms_diff_num(.x, data = emf_data_split$GNJ)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

emf_diff_num_JFL <- map_dfr(emf_y_vars, ~ {
  res <- lmms_diff_num(.x, data = emf_data_split$JFL)
  res$model_res %>% 
    mutate(y_var = .x) %>% 
    rownames_to_column("N_treat")
})

emf_diff_num_GH$Site <- "Boreal forest"
emf_diff_num_DLS$Site <- "Temperate forest"
emf_diff_num_GNJ$Site <- "Subtropical forest"
emf_diff_num_JFL$Site <- "Tropical forest"

emf_diff_num_mul <- rbind(emf_diff_num_GH, emf_diff_num_DLS, emf_diff_num_GNJ, emf_diff_num_JFL)

emf_diff_num_mul$y_var <- factor(emf_diff_num_mul$y_var, 
                                 levels = c("TRGR","TLP"))

emf_diff_num_mul$Site <- factor(emf_diff_num_mul$Site, 
                               levels = c("Boreal forest","Temperate forest", 
                                          "Subtropical forest","Tropical forest"))
emf_diff_num_mul

```

```{r}

ef_var_labels <- data.frame(
  y_var = c("TRGR", "TLP"),
  label = c("(A)", "(B)"))

ef_var_labels$y_var <- factor(ef_var_labels$y_var, levels = levels(emf_diff_num_mul$y_var))

emf_diff_estimate_plot <- ggplot(data = emf_diff_num_mul, 
                                 aes(x = Site, y = Estimate, fill = Site))+
      geom_text(data = ef_var_labels, aes(x = -Inf, y = Inf, label = label), 
            inherit.aes = FALSE, hjust = -0.3, vjust = 1.5, size = 5, 
            fontface = "bold", family = "serif") +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8),  width = 0.7) +  
  geom_errorbar(
    aes(ymin = Estimate - `Std. Error`, 
        ymax = Estimate + `Std. Error`,
        color = Site),  
    position = position_dodge(width = 0.8),  
    width = 0.2, 
    lwd = 0.4) +
    facet_manual(vars(y_var), 
               design = "AB", 
               scales = "free_y",
               axes = "all",
               remove_labels = "all",
               labeller = labeller(y_var = c(
                 TRGR  = "Tree relative growth rate",
                 TLP = "Total litterfall production"
               ))) +
  geom_hline(yintercept = 0, color = "Black", linetype = "dashed", lwd = 0.3) +
  scale_fill_manual(values = c("Boreal forest" = "#08306B", "Temperate forest" = "#2171B5", 
                               "Subtropical forest" = "#6BAED6", "Tropical forest" = "#C6DBEF")) +
  scale_color_manual(values = c("Boreal forest" = "#08306B", "Temperate forest" = "#2171B5", 
                                "Subtropical forest" = "#6BAED6", "Tropical forest" = "#C6DBEF")) +
    scale_x_discrete(limits = c("Boreal forest","Temperate forest",
                              "Subtropical forest","Tropical forest"),
                   labels = c(" "," "," "," "), 
                   position = "bottom") +
  theme_classic() + 
  geom_text(data = emf_diff_num_mul, aes(y = sig_position, x = Site, label = sig),
            position = position_dodge(width = 0.8), vjust = 0.5,  size = 3, color = "black") +
  labs(title = NULL,x = NULL, y = "Effect size")+
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.text = element_text(size = 7),
        legend.title = element_blank(),
         legend.key.size = unit(0.35, "cm"),
        axis.ticks = element_line(color = "black", size = 0.3),
        axis.ticks.length.x = unit(-0.10, "cm"), # 
        #axis.ticks = element_blank(),
        axis.text.x = element_text(size = 10, color = "black",angle = 30, hjust = 1),
        axis.text.y = element_text(size = 10, color = "black"),
        axis.line = element_blank(),
        axis.title.x = element_text(size = 12, color = "black"),
        axis.title.y = element_text(size = 12, color = "black"),
        text = element_text(family = "serif"),
        strip.text = element_text(color = "black", size = 13),
        strip.background = element_rect(colour = NA, fill = "#e9ecef",size = 0.3),
        panel.border = element_rect(color = "black", fill = NA, size = 0.3))

emf_diff_estimate_plot

#ggsave("./Fig_4.pdf", plot = emf_diff_estimate_plot, width = 7, height = 4.5)

```

# 6 Connecting soil env, mic and prod

## 6.1 Correlation analysis

```{r}

plot_half_corr <- function(data, var_names, method = "pearson", title = title) {
  library(dplyr)
  library(reshape2)
  library(ggplot2)
  library(psych)

  testRes <- corr.test(data[, var_names], method = method, adjust = "none")

  cor_res <- melt(testRes$r) %>%
    mutate(
      p = melt(testRes$p, as.is = TRUE)[[3]],
      p_signif = symnum(
        p,
        corr = FALSE,
        na = FALSE,
        cutpoints = c(0, 0.001, 0.01, 0.05, 1),
        symbols = c("***", "**", "*", "")
      )
    ) %>%
    setNames(c("env", "diversity", "r", "p", "p_signif"))

  cor_res_filtered <- cor_res %>%
    mutate(
      env = factor(env, levels = var_names),
      diversity = factor(diversity, levels = var_names)
    ) %>%
    filter(as.numeric(env) > as.numeric(diversity))

  p_out <- ggplot(cor_res_filtered, aes(diversity, env)) +
    geom_tile(color = "black", fill = "white", linewidth = 0.3) +
    geom_point(aes(size = abs(r), color = r), shape = 15) +
    geom_text(aes(label = p_signif), size = 4, color = "black", hjust = 0.5, vjust = 0.7) +
    scale_x_discrete(position = "top") + 
    scale_color_gradientn(
      colours = RColorBrewer::brewer.pal(5, "RdBu"),
      limits = c(-1, 1),
      name = "Pearson\nCorrelation"
    ) +
    labs(title = title) +
    scale_size(range = c(1, 8), guide = NULL) +
    theme_minimal() +
    theme(
      axis.text.x.top = element_text(angle = 45, vjust = 0, hjust = 0),
      panel.grid = element_blank(),
      axis.title = element_blank(),
      axis.ticks = element_blank()
    )

  return(p_out)
}

```

```{r}

var_names <- c("SC","SN", "SP","C_N","C_P","N_P","pH",
               "Chao1_b","Simpson_b","Pielou_b",
               "Chao1_f","Simpson_f","Pielou_f",
               "V_b","E_b", "AD_b","IAPL_b", "Con_b",
               "V_f","E_f", "AD_f", "IAPL_f", "Con_f")

GH_cor_plot <- plot_half_corr(emf_data_split$GH, var_names, title = "Boreal forest")
DLS_cor_plot <- plot_half_corr(emf_data_split$DLS, var_names, title = "Temperate forest")
GNJ_cor_plot <- plot_half_corr(emf_data_split$GNJ, var_names, title = "Subtropical forest")
JFL_cor_plot <- plot_half_corr(emf_data_split$JFL, var_names, title = "Tropical forest")

```

```{r}

fig_S4 <- ggarrange(GH_cor_plot, DLS_cor_plot, 
                    GNJ_cor_plot, JFL_cor_plot,
                    labels=c("(A)","(B)","(C)","(D)"),
                    align = "hv", nrow = 2,ncol = 2)

fig_S4
#ggsave("fig_S4.pdf", plot = fig_S4, width = 14, height = 11)

```

```{r}

mic_env_data_split <- split(emf_data, emf_data$Site)

mic_env_data_split$GH <- mic_env_data_split$GH %>%
  mutate(across(
    .cols = -c(SampleID, Site_Group, Group ,Site, Block),
    .fns = ~ as.numeric(scale(.x, center = TRUE))
  ))

mic_env_data_split$DLS <- mic_env_data_split$DLS %>%
  mutate(across(
    .cols = -c(SampleID, Site_Group,Group ,Site, Block), 
    .fns = ~ as.numeric(scale(.x, center = TRUE))
  ))

mic_env_data_split$GNJ <- mic_env_data_split$GNJ %>%
  mutate(across(
    .cols = -c(SampleID, Site_Group,Group ,Site, Block), 
    .fns = ~ as.numeric(scale(.x, center = TRUE))
  ))

mic_env_data_split$JFL <- mic_env_data_split$JFL %>%
  mutate(across(
    .cols = -c(SampleID,Site_Group, Group ,Site, Block), 
    .fns = ~ as.numeric(scale(.x, center = TRUE))
  ))

```


```{r}

library(lmerTest)

lmer_multi_summary <- function(data, y_var, x_vars, site_name) {

  res_list <- lapply(x_vars, function(x) {

    form <- as.formula(paste0(y_var, " ~ ", x ))
    mod  <- lm(form, data = data)

    sum_mod  <- summary(mod)
    coef_tab <- as.data.frame(sum_mod$coefficients)

    slope <- coef_tab[x, ]

    r2_vals <- MuMIn::r.squaredGLMM(mod)

    p_val <- slope[["Pr(>|t|)"]]

    sig_label <- dplyr::case_when(
      p_val < 0.001 ~ "***",
      p_val < 0.01  ~ "**",
      p_val < 0.05  ~ "*",
      TRUE ~ " "
    )

    data.frame(
      predictor = x,
      estimate  = slope[["Estimate"]],
      std_error = slope[["Std. Error"]],
      t_value   = slope[["t value"]],
      p_value   = p_val,
      signif    = sig_label,
      R2m = r2_vals[1],
      R2c = r2_vals[2],
      Site = site_name
    )
  })

  dplyr::bind_rows(res_list)
}


```

```{r}

DLS_EMF_pca <- PCA(mic_env_data_split$DLS[,c("TRGR","TLP")], scale.unit = F, graph = T)
DLS_EMF_pca_coord <- as.data.frame(DLS_EMF_pca$ind$coord)  
DLS_EMF_pca$var

mic_env_data_split$DLS$EMF_PC1 <- DLS_EMF_pca_coord$Dim.1
#mic_env_data_split$DLS$soil_nutrients_PC2 <- DLS_soil_pca_coord$Dim.2

```

```{r}

GH_TRGR_soil_lmer_res <- lmer_multi_summary(
  data = mic_env_data_split$GH,
  site_name = "Boreal forest",
  y_var = "TRGR",
  x_vars = c("SN", "C_N","N_P","Chao1_b","Chao1_f",
              "V_b","E_b", "AD_b","IAPL_b", "Con_b",
              "V_f","E_f", "AD_f", "IAPL_f")#,
  #random = "Block"#,
  #p_threshold = 0.05
)

DLS_TRGR_soil_lmer_res <- lmer_multi_summary(
  data = mic_env_data_split$DLS,
  site_name = "Temperate forest",
  y_var = "EMF_PC1",
  x_vars = c("SN","C_N","Chao1_b","Chao1_f",
              "Con_b","V_f","E_f", "AD_f", "IAPL_f","Con_f")#,
  #random = "Block"#,
  #p_threshold = 0.05
)

GNJ_TRGR_soil_lmer_res <- lmer_multi_summary(
  data = mic_env_data_split$GNJ,
  site_name = "Subtropical forest",
  y_var = "TLP",
  x_vars = c("SC","SP","C_N","C_P","N_P","Chao1_b","Chao1_f",
              "V_b","E_b", "AD_b","IAPL_b", "Con_b",
             "V_f", "AD_f", "IAPL_f")#,
  #random = "Block"#,
  #p_threshold = 0.05
)

EF_lmer_res <- bind_rows(GH_TRGR_soil_lmer_res,DLS_TRGR_soil_lmer_res,GNJ_TRGR_soil_lmer_res)
EF_lmer_res <- EF_lmer_res %>%
  mutate(
    Category = case_when(
      grepl("_b$", predictor) ~ "Bacterial community",
      grepl("_f$", predictor) ~ "Fungal community",
      TRUE ~ "Soil properties"
    )
  )

EF_lmer_res$Site <- factor(EF_lmer_res$Site, 
                               levels = c("Boreal forest",
                                          "Temperate forest",
                                          "Subtropical forest"))

EF_lmer_res$Category <- factor(EF_lmer_res$Category, 
                               levels = c("Soil properties",
                                          "Bacterial community",
                                          "Fungal community"))

EF_lmer_res$predictor <- factor(EF_lmer_res$predictor, 
                               levels = c(
               "Chao1_b","V_b","E_b", "AD_b","IAPL_b", "Con_b",
               "Chao1_f","V_f","E_f", "AD_f", "IAPL_f", "Con_f",
               "N_P","C_P","C_N","SP","SN", "SC"))


```


```{r}

EF_lmer_labels <- data.frame(
  Site = c("Boreal forest", "Temperate forest","Subtropical forest"),
  label = c("(A)", "(B)", "(C)"))
EF_lmer_labels$Site <- factor(EF_lmer_labels$Site, 
                               levels = c("Boreal forest",
                                          "Temperate forest",
                                          "Subtropical forest"))

Fig_5 <- ggplot(EF_lmer_res, aes(x = estimate, y = predictor, colour = Category)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey40") +
  geom_errorbarh(aes(xmin = estimate - std_error,
                     xmax = estimate + std_error),
                 height = 0.2,
                 size = 0.8,
                 show.legend = FALSE) +
  geom_text(data = EF_lmer_labels, aes(x = -Inf, y = Inf, label = label), 
            inherit.aes = FALSE, hjust = -0.3, vjust = 1.5, size = 5, 
            fontface = "bold", family = "serif") +
  geom_point(size = 4) +
  geom_text(aes(label = signif,x = estimate + 1.5*std_error),size = 6,
            hjust = 0.5,show.legend = FALSE, color = "black") +
  scale_y_discrete(labels = c("Chao1_b" = "Chao1",
                              "V_b" = "V",
                              "E_b" = "E", 
                              "AD_b" = "AD",
                              "IAPL_b" = "IAPL", 
                              "Con_b" = "Con",
                              "Chao1_f" = "Chao1",
                              "V_f" = "V",
                              "E_f" = "E", 
                              "AD_f" = "AD", 
                              "IAPL_f" = "IAPL", 
                              "Con_f" = "Con",
                              "SC" = "SC",
                              "SN" = "SN", 
                              "SP" = "SP",
                              "C_N" = "C:N",
                              "C_P" = "C:P",
                              "N_P" = "N:P")) +
  facet_wrap(~ Site, scales = "free_y") +
    scale_colour_manual(values = c(
      "Soil properties" = "#C6DBEF",
      "Fungal community" = "#6BAED6",
      "Bacterial community" = "#08306B")) +
    labs(x = "Estimate", y = NULL)+
  theme_bw() +
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.text = element_text(size = 10),
        #legend.title = element_text(size = 7),
        axis.text.x = element_text(size = 15, color = "black"),
        axis.text.y = element_text(size = 15, color = "black"),
        axis.ticks = element_line(color = "black", size = 0.3),
        #axis.ticks.length = unit(0.15, "cm"), #
        strip.background.y = element_rect(colour = NA), 
        axis.line = element_blank(),
        axis.title.x = element_text(size = 16, color = "black"),
        axis.title.y = element_text(size = 16, color = "black"),
        text = element_text(family = "serif"),
        strip.text = element_text(color = "black", size = 17),
        strip.background = element_rect(colour = NA, fill = "#e9ecef", size = 0.5),
        panel.border = element_rect(color = "black", fill = NA, size = 0.5),
        panel.grid = element_blank()
        ) 

Fig_5

#ggsave("./Fig_5.pdf", Fig_5, width = 12, height = 7)

```

```{r}

JFL_TRGR_soil_lmer_res <- lmer_multi_summary(
  data = mic_env_data_split$JFL,
  site_name = "Tropical forest",
  y_var = c("TRGR"),
  x_vars = c("SC","C_N","pH","Chao1_b","Chao1_f",
              "V_b","E_b", "AD_b","IAPL_b", "Con_b",
             "V_f","E_f", "AD_f","IAPL_f", "Con_f")#,
  #random = "Block"#,
  #p_threshold = 0.05
)

JFL_TLP_soil_lmer_res <- lmer_multi_summary(
  data = mic_env_data_split$JFL,
  site_name = "Tropical forest",
  y_var = c("TLP"),
  x_vars = c("SC","C_N","pH","Chao1_b","Chao1_f",
              "V_b","E_b", "AD_b","IAPL_b", "Con_b",
             "V_f","E_f", "AD_f","IAPL_f", "Con_f")#,
  #random = "Block"#,
  #p_threshold = 0.05
)

JFL_TRGR_soil_lmer_res$Yar_name <- "TRGR"
JFL_TLP_soil_lmer_res$Yar_name <- "TLP"

JFL_cor_res <- bind_rows(JFL_TRGR_soil_lmer_res, JFL_TLP_soil_lmer_res)

JFL_cor_res <- JFL_cor_res %>%
  mutate(
    Category = case_when(
      grepl("_b$", predictor) ~ "Bacterial community",
      grepl("_f$", predictor) ~ "Fungal community",
      TRUE ~ "Soil properties"
    )
  )

JFL_cor_res$predictor <- factor(JFL_cor_res$predictor, 
                               levels = c(
               "Chao1_b","V_b","E_b", "AD_b","IAPL_b", "Con_b",
               "Chao1_f","V_f","E_f", "AD_f", "IAPL_f", "Con_f",
               "N_P","C_P","C_N","SP","SN", "SC","pH"))

JFL_cor_res$Yar_name <- factor(JFL_cor_res$Yar_name, levels = c("TRGR","TLP"))

```


```{r}

Fig_S5_labels <- data.frame(
  Yar_name = c("TRGR", "TLP"),
  label = c("(A)", "(B)"))
Fig_S5_labels$Yar_name <- factor(Fig_S5_labels$Yar_name, levels = c("TRGR","TLP"))

Fig_S5 <- ggplot(JFL_cor_res, aes(x = estimate, y = predictor, colour = Category)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey40") +
  geom_errorbarh(aes(xmin = estimate - std_error,
                     xmax = estimate + std_error),
                 height = 0.2,
                 size = 0.6,
                 show.legend = FALSE) +
  geom_text(data = Fig_S5_labels, aes(x = -Inf, y = Inf, label = label), 
            inherit.aes = FALSE, hjust = -0.3, vjust = 1.5, size = 5, 
            fontface = "bold", family = "serif") +
  geom_point(size = 3) +
  geom_text(aes(label = signif,x = estimate + 1.5*std_error),size = 4,
            hjust = 0.5,show.legend = FALSE, color = "black") +
  scale_y_discrete(labels = c("Chao1_b" = "Chao1",
                              "V_b" = "V",
                              "E_b" = "E", 
                              "AD_b" = "AD",
                              "IAPL_b" = "IAPL", 
                              "Con_b" = "Con",
                              "Chao1_f" = "Chao1",
                              "V_f" = "V",
                              "E_f" = "E", 
                              "AD_f" = "AD", 
                              "IAPL_f" = "IAPL", 
                              "Con_f" = "Con",
                              "pH" = "pH",
                              "SC" = "SC",
                              "SN" = "SN", 
                              "SP" = "SP",
                              "C_N" = "C:N",
                              "C_P" = "C:P",
                              "N_P" = "N:P")) +
  facet_wrap(~ Yar_name, scales = "free_y") +
    scale_colour_manual(values = c(
      "Soil properties" = "#C6DBEF",
      "Fungal community" = "#6BAED6",
      "Bacterial community" = "#08306B")) +
    labs(x = "Estimate", y = NULL)+
  theme_bw() +
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.text = element_text(size = 10),
        #legend.title = element_text(size = 7),
        axis.text.x = element_text(size = 15, color = "black"),
        axis.text.y = element_text(size = 15, color = "black"),
        axis.ticks = element_line(color = "black", size = 0.3),
        #axis.ticks.length = unit(0.15, "cm"), #
        strip.background.y = element_rect(colour = NA), 
        axis.line = element_blank(),
        axis.title.x = element_text(size = 16, color = "black"),
        axis.title.y = element_text(size = 16, color = "black"),
        text = element_text(family = "serif"),
        strip.text = element_text(color = "black", size = 17),
        strip.background = element_rect(colour = NA, fill = "#e9ecef", size = 0.5),
        panel.border = element_rect(color = "black", fill = NA, size = 0.5),
        panel.grid = element_blank()
        ) 

Fig_S5

#ggsave("Fig_S5.pdf", plot = Fig_S5, width = 10, height = 6)

```

## 6.2 Glmm.hp of measurement func

```{r}

#bacteria
GH_np_b_pca <- PCA(mic_env_data_split$GH[,c("V_b","E_b","AD_b","IAPL_b")], 
                   scale.unit = F, graph = F)

GNJ_np_b_pca <- PCA(mic_env_data_split$GNJ[,c("V_b","E_b","AD_b","IAPL_b","Con_b")], 
                    scale.unit = F, graph = F)


```

```{r}

#fungi
GH_np_f_pca <- PCA(mic_env_data_split$GH[,c("E_f","AD_f")], 
                   scale.unit = F, graph = F)

DLS_np_f_pca <- PCA(mic_env_data_split$DLS[,c("AD_f","IAPL_f","Con_f")], 
                    scale.unit = F, graph = F)

GNJ_np_f_pca <- PCA(mic_env_data_split$GNJ[,c("V_f","AD_f","IAPL_f")], 
                    scale.unit = F, graph = F)

```

```{r}

pca_plot <- function(pca_res = pca_res, data_frame = data_frame, title = title){
  data_frame$Group <- factor(data_frame$Group, levels = c("CK","N50","N100"))
  var_explained <- as.data.frame(pca_res$eig[,c("percentage of variance")])
  pca_plot <- fviz_pca_biplot(pca_res,
                              col.ind = data_frame$Group,
                              palette = c( "#6BAED6", "#2171B5", "#08306B"),#C6DBEF
                              addEllipses = TRUE,
                              alpha = 1,
                              label = "var",
                              col.var = "#595959", 
                              repel = T,
                              legend.title = "N_treat",
                              title = title)+
    labs(x=paste0("PC1 (",round(var_explained[1,]),"%)"),
         y=paste0("PC2 ( ",round(var_explained[2,]),"%)"))+
    scale_y_continuous(expand=c(0,0.5))+
    scale_x_continuous(expand=c(0,0.5))+
    theme_bw() + 
    theme(legend.position = "bottom",
          legend.direction = "horizontal",
          legend.text = element_text(size = 11),
          legend.title = element_blank(),
          panel.grid = element_blank(),
          text = element_text(family = "serif"),
          axis.text.x = element_text(size = 13, color = "black"),
          axis.text.y = element_text(size = 13, color = "black"),
          axis.title.x = element_text(size = 16), 
          axis.title.y = element_text(size = 16))+
    guides(color = guide_legend(override.aes = list(shape = NA)))
  
  return(pca_plot)
  
}

```


```{r}

#bacteria
GH_np_b_pca_plot <- pca_plot(pca_res = GH_np_b_pca, data_frame = mic_env_data_split$GH, 
                             title = "Boreal forest's BNC")
GNJ_np_b_pca_plot <- pca_plot(pca_res = GNJ_np_b_pca, data_frame = mic_env_data_split$GNJ, 
                              title = "Subtropical forest's BNC")
#fungi
GH_np_f_pca_plot <- pca_plot(pca_res = GH_np_f_pca, data_frame = mic_env_data_split$GH, 
                             title = "Boreal forest's FNC")
DLS_np_f_pca_plot <- pca_plot(pca_res = DLS_np_f_pca, data_frame = mic_env_data_split$DLS, 
                              title = "Temperate forest's FNC")
GNJ_np_f_pca_plot <- pca_plot(pca_res = GNJ_np_f_pca, data_frame = mic_env_data_split$GNJ, 
                              title = "Subtropical forest's FNC")

```

```{r}

fig_S6 <- ggarrange(GH_np_b_pca_plot,GH_np_f_pca_plot,DLS_np_f_pca_plot,
                    GNJ_np_b_pca_plot,GNJ_np_f_pca_plot,
                    labels=c("(A)","(B)","(C)","(D)","(E)"),
                    align = "hv", nrow = 2, ncol = 3,
                    common.legend = T, legend = "bottom")
fig_S6

#ggsave("Fig_S6.pdf", plot = fig_S6, width = 11, height = 7)

```

```{r}

mic_env_data_split$GH$BNC_PC1 <- as.data.frame(GH_np_b_pca$ind$coord)$Dim.1
mic_env_data_split$GH$FNC_PC1 <- as.data.frame(GH_np_f_pca$ind$coord)$Dim.1

mic_env_data_split$DLS$FNC_PC1 <- as.data.frame(DLS_np_f_pca$ind$coord)$Dim.1
mic_env_data_split$DLS$FNC_PC2 <- as.data.frame(DLS_np_f_pca$ind$coord)$Dim.2

mic_env_data_split$GNJ$BNC_PC1 <- as.data.frame(GNJ_np_b_pca$ind$coord)$Dim.1
mic_env_data_split$GNJ$BNC_PC2 <- as.data.frame(GNJ_np_b_pca$ind$coord)$Dim.2

mic_env_data_split$GNJ$FNC_PC1 <- as.data.frame(GNJ_np_f_pca$ind$coord)$Dim.1

```

```{r}

if ("package:lmerTest" %in% search()) {
  detach("package:lmerTest", unload = TRUE, character.only = TRUE)
}

GH_TRGR <- glmm.hp(lmer(TRGR ~ BNC_PC1 + FNC_PC1 + Chao1_b + Chao1_f + 
                                  (1|Block), data = mic_env_data_split$GH), 
                       iv = list(Alpha_diversity = c("Chao1_b","Chao1_f"), 
                                 Network_complexity = c("BNC_PC1","FNC_PC1"))
                   )

DLS_EMF <- glmm.hp(lmer(EMF_PC1 ~ Con_b + FNC_PC1 + FNC_PC2 + Chao1_b + Chao1_f + 
                                  (1|Block), data = mic_env_data_split$DLS, na.action = na.omit), 
                       iv = list(Alpha_diversity = c("Chao1_b","Chao1_f"), 
                                 Network_complexity = c("Con_b","FNC_PC1","FNC_PC2"))
                    )

DLS_TRGR <- glmm.hp(lmer(TRGR ~ Con_b + FNC_PC1 + FNC_PC2 + Chao1_b + Chao1_f + 
                                  (1|Block), data = mic_env_data_split$DLS,na.action = na.omit), 
                       iv = list(Alpha_diversity = c("Chao1_b","Chao1_f"), 
                                 Network_complexity = c("Con_b","FNC_PC1","FNC_PC2"))
                    )

DLS_TLP <- glmm.hp(lmer(TLP ~ Con_b + FNC_PC1 + FNC_PC2 + Chao1_b + Chao1_f + 
                                  (1|Block), data = mic_env_data_split$DLS,na.action = na.omit), 
                       iv = list(Alpha_diversity = c("Chao1_b","Chao1_f"), 
                                 Network_complexity = c("Con_b","FNC_PC1","FNC_PC2"))
                    )

GNJ_TLP <- glmm.hp(lmer(TLP ~ BNC_PC1 + FNC_PC1 +  Chao1_b + Chao1_f + 
                                  (1|Block), data = mic_env_data_split$GNJ), 
                                 iv = list(Alpha_diversity = c("Chao1_b","Chao1_f"), 
                                  Network_complexity = c("BNC_PC1","FNC_PC1"))
                   )

```

```{r}

r.squaredGLMM(lmer(TRGR ~ BNC_PC1 + FNC_PC1 + Chao1_b + Chao1_f + (1|Block), 
                   data = mic_env_data_split$GH))

r.squaredGLMM(lmer(TRGR ~ Con_b + FNC_PC1 + FNC_PC2 + Chao1_b + Chao1_f + (1|Block), 
                   data = mic_env_data_split$DLS, na.action = na.omit))

r.squaredGLMM(lmer(TLP ~ Con_b + FNC_PC1 + FNC_PC2 + Chao1_b + Chao1_f + (1|Block), 
                   data = mic_env_data_split$DLS,na.action = na.omit))

r.squaredGLMM(lmer(TLP ~ BNC_PC1 + FNC_PC1 +  Chao1_b + Chao1_f + (1|Block), 
                   data = mic_env_data_split$GNJ))

```

```{r}

glmm_hp_analysis <- function(glmm_hp_res, y_name, group_name) {
  data = glmm_hp_res$hierarchical.partitioning
  result <- as.data.frame(data)
  result$Type <- rownames(result)
  result$Group <- group_name
  result$Var_name <- y_name
  return(result)
}

```

```{r}

GH_TRGR_res <- glmm_hp_analysis(GH_TRGR, "TRGR", "GH")

DLS_EMF_res <- glmm_hp_analysis(DLS_EMF, "EMF", "DLS")

DLS_TRGR_res <- glmm_hp_analysis(DLS_TRGR, "TRGR", "DLS")
DLS_TLP_res <- glmm_hp_analysis(DLS_TLP, "TLP", "DLS")


GNJ_TLP_res <- glmm_hp_analysis(GNJ_TLP, "TLP", "GNJ")

```


```{r}

glmm_hp_res_rf <- rbind(GH_TRGR_res, DLS_EMF_res, DLS_TRGR_res, DLS_TLP_res, GNJ_TLP_res)

glmm_hp_res_rf$y_name <- paste(glmm_hp_res_rf$Var_name, glmm_hp_res_rf$Group, sep = "_")
glmm_hp_res_rf$Var_name <- factor(glmm_hp_res_rf$Var_name, 
                                levels = c("TRGR","TLP"))

glmm_hp_res_rf <- glmm_hp_res_rf %>%
    mutate(Site = case_when(
      grepl("GH", glmm_hp_res_rf$Group) ~ "Boreal forest",
      grepl("DLS", glmm_hp_res_rf$Group) ~ "Temperate forest",
      grepl("GNJ", glmm_hp_res_rf$Group) ~ "Subtropical forest",
      grepl("JFL", glmm_hp_res_rf$Group) ~ "Tropical forest",
      TRUE ~ NA_character_
    ))

glmm_hp_res_rf <- glmm_hp_res_rf %>% rename(I_perc = `I.perc(%)`)
glmm_hp_res_rf

```

## 6.3 Best predict factor for mic

```{r}

options(na.action = "na.fail")

GH_b_net_dredge <- dredge(lmer(BNC_PC1 ~ SN + C_N + N_P + (1|Block),
                                     mic_env_data_split$GH))
GH_f_net_dredge <- dredge(lmer(FNC_PC1 ~ SN + C_N + N_P + (1|Block),
                                     mic_env_data_split$GH))

```

```{r}

DLS_b_net_dredge <- dredge(lmer(Con_b ~ SN + C_N + (1|Block),
                                mic_env_data_split$DLS))

DLS_f_net_pc1_dredge <- dredge(lmer(FNC_PC1 ~ SN + C_N + (1|Block),
                                    mic_env_data_split$DLS))

DLS_f_net_pc2_dredge <- dredge(lmer(FNC_PC2 ~ SN + C_N + (1|Block),
                                   mic_env_data_split$DLS))

```

```{r}

GNJ_b_net_dredge <- dredge(lmer(BNC_PC1 ~ SC + SP + C_N + C_P + N_P + (1|Block),
                                     mic_env_data_split$GNJ))

GNJ_f_net_dredge <- dredge(lmer(FNC_PC1 ~ SC + SP + C_N + C_P + N_P + (1|Block),
                                     mic_env_data_split$GNJ))

```

```{r}

#JFL_b_net_pc1_dredge <- dredge(lmer(net_b_PC1 ~ SC +  C_N + pH + (1|Block),
                                     #mic_env_data_split$JFL))


```


```{r}

plot_dredge_importance <- function(dredge_obj,
                                   low_col  = "#DCEAF7",
                                   high_col = "#1F5AA6",
                                   title = title) {

  #imp <- sw(subset(model.sel(dredge_obj), delta <= 4))
  imp <- sw(dredge_obj)

  imp_df <- data.frame(
    Variable = names(imp),
    Importance = as.numeric(imp)
  )

  imp_df$Variable <- factor(
    imp_df$Variable,
    levels = imp_df$Variable[order(imp_df$Importance)]
  )

  p <- ggplot(imp_df, aes(x = Importance,
                     y = Variable,
                     fill = Importance)) +
    geom_col(width = 0.7) +
    scale_fill_gradient(
      low  = low_col,
      high = high_col
    ) +
    theme_classic() +
    theme(
      legend.position = "none"
    ) +
    labs(
      x = "Relative importance",
      y = NULL,
      title = title
    ) + 
      theme(
        axis.ticks = element_line(color = "black", size = 0.3),
        axis.ticks.length.x = unit(-0.10, "cm"), # 
        #axis.ticks = element_blank(),
        axis.text.x = element_text(size = 10, color = "black"),
        axis.text.y = element_text(size = 10, color = "black"),
        axis.line = element_blank(),
        axis.title.x = element_text(size = 12, color = "black"),
        axis.title.y = element_text(size = 12, color = "black"),
        axis.title = element_text(size = 12, color = "black"),
        text = element_text(family = "serif"),
        strip.text = element_text(color = "black", size = 13),
        strip.background = element_rect(colour = NA, fill = "#e9ecef",size = 0.3),
        panel.border = element_rect(color = "black", fill = NA, size = 0.3))
  return(p)
  
}

```

```{r}

GH_b_net_dredge_impo <- plot_dredge_importance(GH_b_net_dredge, 
                                                   title = "Boreal forest's soil BNC (PC1)")
GH_f_net_dredge_impo <- plot_dredge_importance(GH_f_net_dredge, 
                                                   title = "Boreal forest's soil FNC (PC1)")

DLS_b_net_dredge_impo <- plot_dredge_importance(DLS_b_net_dredge, 
                                                title = "Temperate forest's soil BNC (PC1)")
DLS_f_net_pc1_dredge_impo <- plot_dredge_importance(DLS_f_net_pc1_dredge, 
                                                title = "Temperate forest's soil FNC (PC1)")
DLS_f_net_pc2_dredge_impo <- plot_dredge_importance(DLS_f_net_pc2_dredge, 
                                                title = "Temperate forest soil FNC (PC2)")

GNJ_b_net_dredge_impo <- plot_dredge_importance(GNJ_b_net_dredge, 
                                              title = "Subtropical forest's soil BNC (PC1)")
GNJ_f_net_dredge_impo <- plot_dredge_importance(GNJ_f_net_dredge, 
                                              title = "Subtropical forest's soil FNC (PC1)")

#JFL_b_net_pc1_dredge_impo <- plot_dredge_importance(JFL_b_net_pc1_dredge, 
                                                #title = "Tropical forest soil BNC (PC1)")

```

```{r}

fig_S7 <- ggarrange(GH_b_net_dredge_impo, GH_f_net_dredge_impo, DLS_b_net_dredge_impo,
                   DLS_f_net_pc1_dredge_impo, DLS_f_net_pc2_dredge_impo, GNJ_b_net_dredge_impo,
                   GNJ_f_net_dredge_impo, #JFL_b_net_pc1_dredge_impo,
                    labels=c("(A)","(B)","(C)","(D)",
                             "(E)","(F)","(G)"),
                    align = "hv", nrow = 3, ncol = 3)

fig_S7

#ggsave("./Fig_S7.pdf", plot = fig_S7, width = 10, height = 8)

```

# 7 SEM

## 7.1 Model screening

### 7.1.1 GH

```{r}

library(lmerTest)

GH_pSEM_model <- psem(
  lmer(SN ~ N_treat + (1|Block), data = mic_env_data_split$GH),
  lmer(pH ~ N_treat + (1|Block), data = mic_env_data_split$GH),
  lmer(BNC_PC1 ~ SN  + pH + (1|Block), data = mic_env_data_split$GH),
  lmer(FNC_PC1 ~ SN  + pH + (1|Block), data = mic_env_data_split$GH),
  lmer(TRGR ~  BNC_PC1 + FNC_PC1 + SN  + pH + (1|Block), data = mic_env_data_split$GH),
  SN%~~%pH,
  BNC_PC1%~~%N_treat,
  FNC_PC1%~~%N_treat,
  BNC_PC1%~~%FNC_PC1#,
  #data = mic_env_data_split$GH
)

summary(GH_pSEM_model, .progressBar = F)
plot(GH_pSEM_model)

```

### 7.1.2 DLS

```{r}

DLS_pSEM_model_1 <- psem(
  lmer(C_N ~ N_treat + (1|Block), data = mic_env_data_split$DLS),
  lmer(pH ~ N_treat + (1|Block), data = mic_env_data_split$DLS),
  lmer(Con_b ~ C_N + pH + (1|Block), data = mic_env_data_split$DLS),
  lmer(FNC_PC1 ~ C_N + pH + (1|Block), data = mic_env_data_split$DLS),
  lmer(EMF_PC1 ~ Con_b + FNC_PC1 + C_N + pH + (1|Block), data = mic_env_data_split$DLS),
  Con_b%~~%N_treat,
  FNC_PC1%~~%N_treat,
  #Con_b%~~%Con_f,
  EMF_PC1%~~%N_treat#,
  #data = DLS_TRGR_noNA_semdf
)

summary(DLS_pSEM_model_1, .progressBar = F)
plot(DLS_pSEM_model_1)

```


```{r}

#DLS_TRGR_noNA_semdf <- mic_env_data_split$DLS %>% filter(!is.na(TRGR))

DLS_pSEM_model_2 <- psem(
  lmer(C_N ~ N_treat + (1|Block), data = mic_env_data_split$DLS),
  lmer(pH ~ N_treat + (1|Block), data = mic_env_data_split$DLS),
  lmer(Con_b ~ C_N + pH + (1|Block), data = mic_env_data_split$DLS),
  lmer(FNC_PC2 ~ C_N + pH + (1|Block), data = mic_env_data_split$DLS),
  lmer(EMF_PC1 ~ Con_b + FNC_PC2 + C_N + pH + (1|Block), data = mic_env_data_split$DLS),
  Con_b%~~%N_treat,
  FNC_PC2%~~%N_treat,
  Con_b%~~%FNC_PC2,
  EMF_PC1%~~%N_treat#,
  #data = DLS_TRGR_noNA_semdf
)

summary(DLS_pSEM_model_2, .progressBar = F)
plot(DLS_pSEM_model_2)

```

```{r}

#DLS_TRGR_noNA_semdf <- mic_env_data_split$DLS %>% filter(!is.na(TRGR))

DLS_pSEM_model_3 <- psem(
  lmer(SN ~ N_treat + (1|Block), data = mic_env_data_split$DLS),
  lmer(pH ~ N_treat + (1|Block), data = mic_env_data_split$DLS),
  lmer(Con_b ~ SN + pH + (1|Block), data = mic_env_data_split$DLS),
  lmer(FNC_PC1 ~ SN + pH + (1|Block), data = mic_env_data_split$DLS),
  lmer(EMF_PC1 ~ Con_b + FNC_PC1 + SN + pH + (1|Block), data = mic_env_data_split$DLS),
  Con_b%~~%N_treat,
  FNC_PC1%~~%N_treat,
  Con_b%~~%FNC_PC1,
  EMF_PC1%~~%N_treat#,
  #data = DLS_TRGR_noNA_semdf
)

summary(DLS_pSEM_model_3, .progressBar = F)
plot(DLS_pSEM_model_3)

```

```{r}

DLS_pSEM_model_4 <- psem(
  lmer(SN ~ N_treat + (1|Block), data = mic_env_data_split$DLS),
  lmer(pH ~ N_treat + (1|Block), data = mic_env_data_split$DLS),
  lmer(Con_b ~ SN + pH + (1|Block), data = mic_env_data_split$DLS),
  lmer(FNC_PC2 ~ SN + pH + (1|Block), data = mic_env_data_split$DLS),
  lmer(EMF_PC1 ~ Con_b + FNC_PC2 + SN + pH + (1|Block), data = mic_env_data_split$DLS),
  Con_b%~~%N_treat,
  FNC_PC2%~~%N_treat,
  Con_b%~~%FNC_PC2,
  EMF_PC1%~~%N_treat#,
  #data = DLS_TRGR_noNA_semdf
)

summary(DLS_pSEM_model_4, .progressBar = F)
plot(DLS_pSEM_model_4)

```

### 7.1.3 GNJ

```{r}

GNJ_pSEM_model <- psem(
  lmer(C_P ~ N_treat + (1|Block), data = mic_env_data_split$GNJ),
  lmer(pH ~ N_treat + (1|Block), data = mic_env_data_split$GNJ),
  lmer(BNC_PC1 ~ C_P + pH + (1|Block),  data = mic_env_data_split$GNJ),
  lmer(FNC_PC1 ~ C_P + pH + (1|Block), data = mic_env_data_split$GNJ),
  lmer(TLP ~ BNC_PC1 + FNC_PC1 + C_P + pH + (1|Block), 
       data = mic_env_data_split$GNJ),
  C_P%~~%pH,
  BNC_PC1%~~%N_treat,
  FNC_PC1%~~%N_treat,
  BNC_PC1%~~%FNC_PC1#,
  #data = mic_env_data_split$GNJ
)

summary(GNJ_pSEM_model, .progressBar = F)
plot(GNJ_pSEM_model)

```

## 7.2 Direct and indirect effects analysis

```{r}

sem_effects <- function(pSEM_model) {
  coefs_df <- coefs(pSEM_model)
  coefs_df$effect <- coefs_df$Std.Estimate

  build_graph <- function(coefs_df) {
    graph <- list()
    for (i in 1:nrow(coefs_df)) {
      from <- coefs_df$Predictor[i]
      to <- coefs_df$Response[i]
      effect <- coefs_df$effect[i]
      
      if (!is.null(graph[[from]])) {
        graph[[from]] <- rbind(
          graph[[from]],
          data.frame(to = to, effect = effect, stringsAsFactors = FALSE)
        )
      } else {
        graph[[from]] <- data.frame(to = to, effect = effect, stringsAsFactors = FALSE)
      }
    }
    return(graph)
  }
  
  find_all_paths <- function(graph, start, end, visited = character()) {
    if (start == end) return(list(c(end)))
    if (!start %in% names(graph)) return(list())
    
    visited <- c(visited, start)
    paths <- list()
    for (i in 1:nrow(graph[[start]])) {
      next_node <- graph[[start]]$to[i]
      if (!(next_node %in% visited)) {
        sub_paths <- find_all_paths(graph, next_node, end, visited)
        for (sp in sub_paths) {
          paths <- c(paths, list(c(start, sp)))
        }
      }
    }
    return(paths)
  }
  
  path_effect <- function(path, graph) {
    eff <- 1
    for (i in 1:(length(path) - 1)) {
      from <- path[i]
      to <- path[i + 1]
      edge <- graph[[from]]
      eff <- eff * edge$effect[edge$to == to]
    }
    return(eff)
  }

  compute_all_effects <- function(graph, variables) {
    results <- data.frame(
      predictor = character(),
      outcome = character(),
      direct_effect = numeric(),
      indirect_effect = numeric(),
      total_effect = numeric(),
      stringsAsFactors = FALSE
    )
    
    for (pred in variables) {
      for (outc in variables) {
        if (pred != outc) {
          paths <- find_all_paths(graph, pred, outc)
          if (length(paths) > 0) {
            direct_effect <- 0
            indirect_effect <- 0
            
            for (p in paths) {
              eff <- path_effect(p, graph)
              if (length(p) == 2) {
                direct_effect <- direct_effect + eff
              } else {
                indirect_effect <- indirect_effect + eff
              }
            }
            
            total_effect <- direct_effect + indirect_effect
            
            results <- rbind(results, data.frame(
              predictor = pred,
              outcome = outc,
              direct_effect = direct_effect,
              indirect_effect = indirect_effect,
              total_effect = total_effect,
              stringsAsFactors = FALSE
            ))
          }
        }
      }
    }
    return(results)
  }

  graph <- build_graph(coefs_df)
  all_vars <- unique(c(coefs_df$Response, coefs_df$Predictor))
  effects_df <- compute_all_effects(graph, all_vars)
  
  return(effects_df)
}


```

```{r}

GH_SEM_effects_df <- sem_effects(GH_pSEM_model)
DLS_SEM_effects_df <- sem_effects(DLS_pSEM_model_2)
GNJ_SEM_effects_df <- sem_effects(GNJ_pSEM_model)

```


```{r}

GH_SEM_effects_res <- GH_SEM_effects_df %>% 
  filter(outcome == "TRGR") %>%
  pivot_longer(
    cols = c(direct_effect, indirect_effect, total_effect),
    names_to = "effect_type",
    values_to = "effect_size"
  ) %>% filter(predictor != "N_treat")

DLS_SEM_effects_res <- DLS_SEM_effects_df %>% 
  filter(outcome == "EMF_PC1") %>%
  pivot_longer(
    cols = c(direct_effect, indirect_effect, total_effect),
    names_to = "effect_type",
    values_to = "effect_size"
  ) %>% filter(predictor != "N_treat")

GNJ_SEM_effects_res <- GNJ_SEM_effects_df %>% 
  filter(outcome == "TLP") %>%
  pivot_longer(
    cols = c(direct_effect, indirect_effect, total_effect),
    names_to = "effect_type",
    values_to = "effect_size"
  ) %>% filter(predictor != "N_treat")


```

```{r}

TRGR_GH_sem_effects_res <- GH_SEM_effects_res %>% filter(effect_type %in% c("direct_effect"))%>% 
                           mutate(Site = "Boreal forest") 

EMF_DLS_sem_effects_res <- DLS_SEM_effects_res %>% filter(effect_type %in% c("direct_effect"))%>% 
                           mutate(Site = "Temperate forest") 

TLP_GNJ_sem_effects_res <- GNJ_SEM_effects_res %>% filter(effect_type %in% c("direct_effect"))%>% 
                           mutate(Site = "Subtropical forest") 


```

```{r}

sem_effects_res_rf <- rbind(TRGR_GH_sem_effects_res, EMF_DLS_sem_effects_res, 
                            TLP_GNJ_sem_effects_res)

sem_effects_res_rf <- sem_effects_res_rf %>%
    mutate(Var_name = case_when(
  predictor %in% c("SN", "C_N", "C_P") ~ "Soil nutrient",
  predictor == "pH" ~ "Soil pH",
  predictor %in% c("Con_b", "BNC_PC1") ~ "BNC",
  predictor %in% c("FNC_PC1", "FNC_PC2") ~ "FNC",
  TRUE ~ NA_character_
))

sem_effects_res_rf$effect_sign <- ifelse(sem_effects_res_rf$effect_size > 0, "+", 
                         ifelse(sem_effects_res_rf$effect_size < 0, "-", "0"))

sem_effects_res_rf$effect_size_abs <- abs(sem_effects_res_rf$effect_size)

#sem_effects_res_rf$outcome <- factor(sem_effects_res_rf$outcome, levels = c("TRGR","TLP","MBC","MBN"))
#sem_effects_res_rf$Var_name <- factor(sem_effects_res_rf$Var_name, levels = c("BNC","FNC","SN (PC2)","SN (PC1)"))

sem_effects_res_rf

```

```{r}

sem_effects_labels <- data.frame(
  Site = c("Boreal forest","Temperate forest","Subtropical forest"),
  label = c("(A)", "(B)", "(C)"))

sem_effects_res_rf$Site <- factor(sem_effects_res_rf$Site, 
                                  levels = c("Boreal forest","Temperate forest","Subtropical forest"))

sem_effects_labels$Site <- factor(sem_effects_labels$Site, 
                                 levels = c("Boreal forest","Temperate forest","Subtropical forest"))

sem_effects_res_plot <- ggplot(data = sem_effects_res_rf, 
                    aes(x = Var_name, y = effect_size_abs, fill = Var_name)) +
      geom_text(data = sem_effects_labels, aes(x = -Inf, y = Inf, label = label), 
              inherit.aes = FALSE, hjust = -0.3, vjust = 1.5, size = 5, 
              fontface = "bold", family = "serif") +
  geom_bar(stat = "identity", position = "stack", width = 0.7) +
    geom_text(aes(label = effect_sign, y = effect_size_abs), 
            position = position_stack(vjust = 0.5),  #
            size = 3.5, color = "#404040", fontface = "bold") +
  facet_wrap(.~ Site, scales = "free_x", ncol = 3, strip.position = "top") +
  scale_fill_manual(values = c(
    "Soil nutrient" = alpha("#C6DBEF", 0.7),
    "Soil pH" = alpha("#6BAED6", 0.7),
    "BNC" = alpha("#08306B", 0.7),
    "FNC" = alpha("#2171B5", 0.7)
    )) +
  theme_classic() + 
  labs(title = NULL, x = NULL, y = "Direct effect") +
  theme(legend.position = "none",
        #legend.direction = "horizontal",
        #legend.text = element_text(size = 13),
        #legend.title = element_blank(),
        #plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"),
        axis.text.x = element_text(size = 12, color = "black",angle = 30, hjust = 1),
        axis.text.y = element_text(size = 12, color = "black"),
        axis.line = element_blank(), #
        axis.title.x = element_text(size = 14, color = "black"),
        axis.title.y = element_text(size = 14, color = "black"),
        text = element_text(family = "serif"),
        strip.background = element_rect(colour = NA, fill = "#e9ecef", size = 0.5),
        strip.text = element_text(angle = 0, size = 12), #
        panel.border = element_rect(color = "black", fill = NA, size = 0.5))

sem_effects_res_plot

#ggsave("./Fig_S8.pdf", sem_effects_res_plot, width = 10, height = 5)

```



















