---
title: "Rhizosphere Priming Effects on Ancient Peat Carbon"
subtitle: "PRIMETIME Analysis - EGUSPHERE"
author: "Louis A. Mielke, Tom C. Parker, Lorna E. Street, et al."
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_depth: 2
    toc_float: true
    number_sections: true
    theme: united
    highlight: tango
---

```{r setup}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
# Set seed for reproducibility
set.seed(2021)
```

#Objectives - The primary goal of this analysis is to evaluate the distributions of measured and modeled carbon fluxes. Specifically:

Total respiration and raw 14C02 signatures from 2mm mesh (R14C) and 1um (H14C) bags (Figure 2) with and without root and rhizosphere processes, respectively.

Compare modeled Ancient carbon flux based on atmospheric 14C signatures from 2017 and 2021 across different vegetation types (Birch vs. Willow) and treatments (Control vs. Girdled) (see Figure 3).

Quantify the Rhizosphere Priming Effect (RPE) and compare across vegetation and treatments (Figure 5).

Supplementary Information

#Script Setup & Data Loading

First, we load the necessary libraries and import the dataset.

```{r load-packages}
# Load required packages
library(dplyr);library(tidyverse);library(ggplot2);library(ggridges)
library(ggExtra);library(gridExtra);library(ggpubr);library(grDevices)
library(nlme);library(openxlsx);library(knitr);library(kableExtra);library(ggpmisc)
```

#Naming Conventions:

just a note that %modern, PMC, and 14C signature are used interexchangebly when annotating

For df_long:
-Tsoil_in is soil temp (C) 5cm inside the core;
-Tsoil_out is soil temp (C) 30cm outside the core;
-M_mv_in is soil moisture 5cm inside the core;
-M_mv_out is soil moisture 30cm outside the core;

```{r load-data}
# Note: Ensure these files are in your working directory
df <- read.csv("PRIMETIME_IngrowthCores2.csv", header=T, row.names=1) #CO2 data
df_long <- read.csv("PRIMETIME_IngrowthCores_long2.csv", header=T) #CO2 data long
doc_long <- read.csv("PRIMETIME_IngrowthCores_DO14C_long.csv", header=T) #DOC Data

# Basic cleaning
df <- df %>% rownames_to_column(var = "Plot") %>% na.omit()

# Log transformations for Flux data
df$logHflux = log(df$Hflux)
df$logRflux = log(df$Rflux)

# Vectorizing estimates for atmospheric (air) and peat 14CO2 signatures

df$air14C2021 = rnorm(mean = 100.12, sd = .46, n = nrow(df))
df$air14C2018 = rnorm(mean = 101.13, sd = .46, n = nrow(df))
df$air14C2017 = rnorm(mean = 101.80, sd = .46, n = nrow(df))
df$udd_14C    = rnorm(mean = 71.41,  sd = .33, n = nrow(df)) #bulk peat
df$udd_14CO2    = rnorm(mean = 78.65,  sd = .36, n = nrow(df)) #CO2 from peat

```

# Calculate Modeled Fluxes


For the dataset 'df:'

-H14C  is the %modern of the respiration in 1um mesh bags (Exclusion); 
-R14C  is the %modern with root ingrowth (Ingrowth);
-air14C represents the modelled %modern of the ancient carbon released after removing autotrophic inputs into the heterotroph cores (Exclusion); 
-air14C20XX is the %modern of the atmosphere (followed by the year 20XX)
-frecent  is the fraction of recent sources (atmosphere through belowground allocation or litter leaching) to respiration
-fancient   is the fraction of respiration from the ancient peat
-Hflux21R   is the modelled 'heterophic' flux of ancient peat from 2mm soil mesh cores
-Hflux_adj is the modelled heterophic flux of ancient peat from soil mesh cores in the 1um mesh cores
-RecentfluxH  is from the ingrowth core with 1um mesh (Excluson);
-RecentfluxR  is from the ingrowth core with 2 mm mesh allowing roots, fungi and other organisms in (Ingrowth);
-Ancientflux is heterophic respiration flux from ancient peat in 1um or 2mm soil mesh core mod_df_long;
-Recentflux is recent respiration flux from atmospheric inputs in 1um or 2mm soil mesh core for mod_df_long

Method:

We calculated the fractional contribution of recent respiration (frecent) in to the CO2 efflux in the “Ingrowth” and “Excluson” soil cores using equation 2:

(2a) fancient = (∆sample-∆air)/(∆u-∆air)	or in this case fancient = (resp14C - air14CO2) / (udd14CO2 - air14CO2)). 

(2b) fancient = the ancient carbon released taking out recent inputs (frecent = 1-fancient)

where ∆sample or 'R14C' is the 14CO2 content of respiration from the Ingrowth (2mm) mesh cores  and 'H14C' is the 14CO2 content of respiration from the Exclusion (1um) mesh cores and ∆air 'air14C2021' is the 14C content of  recently photosynthesised carbon determined by the 14C signature of the atmosphere in 2021 (100.12 ± 0.49% modern (± analytical error).

fancient21 was determined to be the product of the flux of the +RRP peat core (Ft) and the estimated fraction of the heterotrophic CO2 efflux based on photosynthates from the 2021 atmosphere.

Hflux_adj = totalFlux*fancient

```{r calculations}
mod_df <- df %>%
  group_by(Vegetation, Treatment, Plot) %>%
  mutate(
    # 2021 Partitioning ancient (peat) and recent respiration in 2mm soil mesh cores
    fancient21R = (R14C - air14C2021) / (udd_14CO2 - air14C2021),
    frecent21R = 1 - fancient21R,
    Hflux21R = fancient21R*Rflux,
    RecentRFlux21 = frecent21R*Rflux,
    # 2021 Partitioning in 1um soil mesh cores (adjusted for recent inputs)
    fancient21H = (H14C - air14C2021) / (udd_14CO2 - air14C2021) ,
    frecent21H = 1 - fancient21H,
    Hflux21_adj = fancient21H*Hflux,
    RecentHFlux21 = frecent21H*Hflux,
    # 2017 Partitioning ancient (peat) and recent respiration in 2mm soil mesh cores
    fancient17R = (R14C - air14C2017) / (udd_14CO2 - air14C2017),
    frecent17R = 1 - fancient17R,
    Hflux17R = fancient17R*Rflux,
    RecentRFlux17 = frecent17R*Rflux,
    # 2017 Partitioning in 1um soil mesh cores (adjusted for recent inputs)
    fancient17H = (H14C - air14C2017) / (udd_14CO2 - air14C2017) ,
    frecent17H = 1 - fancient17H,
    Hflux17_adj = fancient17H*Hflux,
    RecentHFlux17 = frecent17H*Hflux,
    # Partitioning in 2mm cores (adjusted for bulk Uddjaure)
    fancient21Rb = (R14C - air14C2021) / (udd_14C - air14C2021),
    frecent21Rb = 1 - fancient21Rb,
    Hflux21Rb = fancient21Rb*Rflux,
    RecentRFluxb = frecent21Rb*Rflux,
    #Partitioning in 1um soil mesh cores (adjusted for bulk Uddjaure)
    fancient21Hb = (H14C - air14C2021) / (udd_14C - air14C2021) ,
    frecent21Hb = 1 - fancient21Hb,
    Hflux21b_adj = fancient21Hb*Hflux,
    RecentHFluxb = frecent21Hb*Hflux,
    # 2017 Partitioning ancient (bulk peat) in 2mm soil mesh cores
    fancient17Rb = (R14C - air14C2017) / (udd_14C - air14C2017),
    frecent17Rb = 1 - fancient17Rb,
    Hflux17Rb = fancient17Rb*Rflux,
    RecentRFlux17b = frecent17Rb*Rflux,
    # 2017 Partitioning in 1um soil mesh cores (adjusted for bulk peat)
    fancient17Hb = (H14C - air14C2017) / (udd_14C - air14C2017) ,
    frecent17Hb = 1 - fancient17Hb,
    Hflux17b_adj = fancient17Hb*Hflux,
    RecentHFlux17b = frecent17Hb*Hflux,
  ) %>%
  ungroup()

# ------------- identify factors instead of characters
df_long$Mesh <- as.factor(df_long$Mesh)
df_long$Treatment <- as.factor(df_long$Treatment)
df_long$Treatment <- as.factor(df_long$Treatment)
df_long$Plot <- as.factor(df_long$Plot)

df_long$atm14C2021 = rnorm(mean = 100.12, sd = .46,n = 42) #2021 distribution of %modern of atmospheric carbon estimate vectorized
df_long$atm14C2018 = rnorm(mean = 101.13, sd = .46,n = 42) #2018 %modern of atmospheric carbon estimate vectorized
df_long$atm14C2017 = rnorm(mean = 101.80, sd = .46,n = 42) #2017 %modern of atmospheric carbon estimate vectorized
df_long$udd_14CO2 = rnorm(mean = 78.65, sd = .36,n = 42) # %modern of ancient soil organic carbon respired (Uddjuare peat) estimate vectorized
df_long$udd_14C = rnorm(mean = 71.41, sd = .33,n = 42) # %modern of ancient soil organic carbon (Uddjuare peat) estimate vectorized

#and again with the long pivot dataset
mod_df_long <- df_long %>%
  group_by(Vegetation, Treatment) %>% #group by veg (birch / willow) and treatment (ctrl / girdled)
  mutate(fancient17 = (X14CO2 - atm14C2017) / (udd_14CO2 - atm14C2017)) %>% #estimated fraction of autotrophic respiration in heterotroph only mesh bags from a mixing model (2017) based on respiration %modern of bulk peat
  mutate(frecent17 = 1-fancient17) %>% #fraction of uddjaure respiration based on the mixing model (2017) based on respiration %modern of bulk peat
  mutate(AncientFlux17 = fancient17*Flux) %>% #estimated flux of uddjaure respiration in 2mm ingrowth core (2017) based on respiration %modern of bulk peat
  mutate(RecentFlux17 = frecent17*Flux) %>%
  mutate(fancient21 = (X14CO2 - atm14C2021) / (udd_14CO2- atm14C2021)) %>% #estimated fraction of autotrophic respiration from a mixing model (2021) based on respiration %modern of bulk peat
  mutate(AncientFlux21 = fancient21*Flux) %>% #estimated flux of uddjaure respiration in 2mm ingrowth core (2021) based on respiration %modern of bulk peat
  mutate(frecent21 = 1-fancient21) %>% #fraction of uddjaure respiration based on the mixing model (2017) and respiration %modern
  mutate(RecentFlux21 = frecent21*Flux) %>%
  mutate(fancient17b = (X14CO2 - atm14C2017) / (udd_14C - atm14C2017)) %>% #estimated fraction of autotrophic respiration in heterotroph only mesh bags from bulk peat in a mixing model (2017)
  mutate(frecent17b = 1-fancient17b) %>% #fraction of uddjaure respiration from bulk peat based on the mixing model (2017)
  mutate(AncientFlux17b = fancient17b*Flux) %>% #estimated flux of uddjaure respiration from bulk peat in 2mm ingrowth core (2017)
  mutate(RecentFlux17b = frecent17b*Flux) %>%
  mutate(fancient21b = (X14CO2 - atm14C2021) / (udd_14C- atm14C2021)) %>% #estimated fraction of autotrophic respiration from a mixing model (2021)
  mutate(AncientFlux21b = fancient21b*Flux) %>% #estimated flux of uddjaure respiration from bulk peat in 2mm ingrowth core (2021)
  mutate(frecent21b = 1-fancient21b) %>% #fraction of uddjaure respiration based on the mixing model (2017) and bulk peat
  mutate(RecentFlux21b = frecent21b*Flux) %>%
  ungroup()
mod_df_long

birch <- mod_df_long[mod_df_long$Vegetation == "Birch",]
willow <- mod_df_long[mod_df_long$Vegetation == "Willow",]

```
 
# Figure 2 14CO2 and total CO2 efflux 
``` {r plot CO2}

# Plotting Birch 14CO2

b14CO2 <- ggplot(data = birch, aes(x = Mesh, y = X14CO2, fill = Vegetation)) +
  geom_abline(intercept = 100.12, slope = 0, color = 'lightblue1',size=1)+
  geom_abline(intercept = 78.65, slope = 0, color = 'red3',size =1)+
  geom_jitter(width = 0.05, color = "#006600", size = 4, alpha = 0.5) +  # Jittered points
  facet_wrap(~interaction(Treatment), ncol=2)+
  scale_fill_manual(values=c("#006600"))+ # blue "#000066"
  ylab("% modern")+
  scale_y_continuous(limits = c(77, 103), expand = c(0, 0)) +
  xlab(" ")+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  theme_classic()+
  scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))+
  theme(plot.title = element_text(hjust = 0.6),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=14, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        strip.text = element_text(size = 14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"))


#14Co2 (%modern) willow this is 13C corrected for air contamination
w14CO2 <- ggplot(data = willow, aes(x = Mesh, y = X14CO2, fill = Vegetation)) +
  geom_abline(intercept = 100.12, slope = 0, color = 'lightblue1',size=1)+
  geom_abline(intercept = 78.65, slope = 0, color = 'red3',size =1)+
  geom_jitter(width = 0.05, color = "#000066", size = 4, alpha = 0.6) +  # Jittered points
  facet_wrap(~interaction(Treatment), ncol=2) +
  scale_fill_manual(values=c("#000066")) + # blue
  ylab("% modern")+
  xlab(" ") +
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  scale_y_continuous(limits = c(77, 103), expand = c(0, 0)) +
   scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))+
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=14, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        strip.text = element_text(size = 14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"))


#total Co2 efflux birch 
  bCO2 <- ggplot(data = birch, aes(x = Mesh, y = Flux, fill = Vegetation)) +
  geom_jitter(width = 0.05, color = "#006600", size = 4, alpha = 0.6) +  # Jittered points
  facet_wrap(~interaction(Treatment), ncol=2)+
  scale_fill_manual(values=c("#006600"))+ # blue
  ylab(expression(CO[2] ~ "efflux" ~ (mu*mol ~ m^{-2} ~ s^{-1})))+
  xlab(" ")+
  theme_classic()+
  scale_y_continuous(limits = c(0, 8.5))+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))+
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=14, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        strip.text = element_text(size = 14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"))

#total CO2 efflux willow
wCO2 <- ggplot(data = willow, aes(x = Mesh, y = Flux, fill = Vegetation)) +
  geom_jitter(width = 0.05, color = "#000066", size = 4, alpha = 0.6)+
  facet_wrap(~interaction(Treatment), ncol=2)+
  scale_fill_manual(values=c("#000066"))+ #blue
  ylab(expression(CO[2] ~ "efflux" ~ (mu*mol ~ m^{-2} ~ s^{-1})))+
  xlab(" ")+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=14, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        strip.text = element_text(size = 14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"))+
  scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))+
  scale_y_continuous(limits = c(0, 4.5))

#combine plots
allco2 <- ggarrange(bCO2, b14CO2, wCO2, w14CO2,
                    labels = c("a)", "b)", "c)","d)"),
                    ncol = 2,nrow=2,
                    common.legend = TRUE)
allco2
ggsave(allco2, filename = "Rplot_CO2flux_14CO2_CtrlGirdled_BirchWillow_Points.svg", device = "svg", height = 8, width = 12, units = "in") 

```

##  CO2 Statistics & Sensitivity Analysis
```{r co2 statistics}

#Birch - comparing heterotrophic and root ingrowth differences in 14CO2 and Conc based on gridling and mesh

#Total respiration linear mixed model (Package nlme version 3.1-166)
birchCO2_lme <-lme(log(Flux) ~ Mesh*Treatment,
                   data = birch,
                   random = ~ 1|Plot)
plot(resid(birchCO2_lme))
qqnorm(resid(birchCO2_lme))
birchCO2_effects <- summary(birchCO2_lme)$tTable
birchCO2_df <- round(as.data.frame(birchCO2_effects),3)
colnames(birchCO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(birchCO2_df,"birchCO2_lme")

#14CO2 signature linear mixed model (Package nlme version 3.1-166)
birch14CO2_lme <- lme(log(X14CO2) ~ Mesh*Treatment,
                      data = birch,
                      random = ~ 1|Plot)
plot(resid(birch14CO2_lme))
qqnorm(resid(birch14CO2_lme))
birch14CO2_effects <- summary(birch14CO2_lme)$tTable
birch14CO2_df <- round(as.data.frame(birch14CO2_effects),3)
colnames(birch14CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(birch14CO2_df,"birch14CO2_lme")

#ancient peat-derived uddjaure respiration (after partitioning based on 2021 atmosphere)
birch_anc_CO2_lme <-lme(log(AncientFlux21) ~ Mesh*Treatment,
                        data = birch,
                        random = ~ 1|Plot)
plot(resid(birch_anc_CO2_lme))
qqnorm(resid(birch_anc_CO2_lme))
birch_anc_CO2_effects <- summary(birch_anc_CO2_lme)$tTable
birch_anc_CO2_df <- round(as.data.frame(birch_anc_CO2_effects),3)
colnames(birch_anc_CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(birch_anc_CO2_df,"birch_anc_CO2_lme")


#ancient peat-derived uddjaure respiration (after partitioning based on 2017 atmosphere)
birch17_anc_CO2_lme <-lme(log(AncientFlux17) ~ Mesh*Treatment,
                        data = birch,
                        random = ~ 1|Plot)
plot(resid(birch17_anc_CO2_lme))
qqnorm(resid(birch17_anc_CO2_lme))
birch17_anc_CO2_effects <- summary(birch17_anc_CO2_lme)$tTable
birch17_anc_CO2_df <- round(as.data.frame(birch17_anc_CO2_effects),3)
colnames(birch17_anc_CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(birch17_anc_CO2_df,"birch17_anc_CO2_lme")

#recent plant-derived respiration (after partitioning) 2021
birch_mod_CO2_lme <-lme(log(RecentFlux21) ~ Mesh*Treatment,
                        data = birch,
                        random = ~ 1|Plot)
plot(resid(birch_mod_CO2_lme))
qqnorm(resid(birch_mod_CO2_lme))
birch_mod_CO2_effects <- summary(birch_mod_CO2_lme)$tTable
birch_mod_CO2_df <- round(as.data.frame(birch_mod_CO2_effects),3)
colnames(birch_mod_CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(birch_mod_CO2_df,"birch_mod_CO2_lme")

#modern respiration (after partitioning based on 2017 atmosphere)
birch17_mod_CO2_lme <-lme(log(RecentFlux17) ~ Mesh*Treatment,
                        data = birch,
                        random = ~ 1|Plot)
plot(resid(birch17_mod_CO2_lme))
qqnorm(resid(birch17_mod_CO2_lme))
birch17_mod_CO2_effects <- summary(birch17_mod_CO2_lme)$tTable
birch17_mod_CO2_df <- round(as.data.frame(birch17_mod_CO2_effects),3)
colnames(birch17_mod_CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(birch17_mod_CO2_df,"birch17_mod_CO2_lme")

# Willow - CO2 Concentration linear mixed model (Package nlme version 3.1-166)
willowCO2_lme <-lme(log(Flux) ~ Mesh*Treatment,
                    data = willow,
                    random = ~ 1|Plot)
plot(resid(willowCO2_lme))
qqnorm(resid(willowCO2_lme))
willowCO2_effects <- summary(willowCO2_lme)$tTable
willowCO2_df <- round(as.data.frame(willowCO2_effects),3)
colnames(willowCO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(willowCO2_df,"willowCO2_lme")

#14CO2 signature linear mixed model (Package nlme version 3.1-166)
willow14CO2_lme <- lme(log(X14CO2) ~ Mesh*Treatment,
                       data = willow,
                       random = ~ 1|Plot)
summary(willow14CO2_lme)
plot(resid(willow14CO2_lme))
qqnorm(resid(willow14CO2_lme))
willow14CO2_effects <- summary(willow14CO2_lme)$tTable
willow14CO2_df <- round(as.data.frame(willow14CO2_effects),3)
colnames(willow14CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(willow14CO2_df, "willow14CO2_lme" )

#ancient uddjaure respiration (after partitioning) 2021 atmosphere
willow_anc_CO2_lme <-lme(log(AncientFlux21) ~ Mesh*Treatment,
                        data = willow,
                        random = ~ 1|Plot)
plot(resid(willow_anc_CO2_lme))
qqnorm(resid(willow_anc_CO2_lme))
willow_anc_CO2_effects <- summary(willow_anc_CO2_lme)$tTable
willow_anc_CO2_df <- round(as.data.frame(willow_anc_CO2_effects),3)
colnames(willow_anc_CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(willow_anc_CO2_df,"willow_ancient_CO2_lme")

#ancient uddjaure respiration (after partitioning based on 2017)
willow17_anc_CO2_lme <-lme(log(AncientFlux17) ~ Mesh*Treatment,
                         data = willow,
                         random = ~ 1|Plot)
plot(resid(willow17_anc_CO2_lme))
qqnorm(resid(willow17_anc_CO2_lme))
willow17_anc_CO2_effects <- summary(willow17_anc_CO2_lme)$tTable
willow17_anc_CO2_df <- round(as.data.frame(willow17_anc_CO2_effects),3)
colnames(willow17_anc_CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(willow17_anc_CO2_df,"willow_ancient17_CO2_lme")


#modern respiration (after partitioning)
willow_mod_CO2_lme <-lme(log(RecentFlux21) ~ Mesh*Treatment,
                        data = willow,
                        random = ~ 1|Plot)
plot(resid(willow_mod_CO2_lme))
qqnorm(resid(willow_mod_CO2_lme))
willow_mod_CO2_effects <- summary(willow_mod_CO2_lme)$tTable
willow_mod_CO2_df <- round(as.data.frame(willow_mod_CO2_effects),3)
colnames(willow_mod_CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(willow_mod_CO2_df,"willow_recent21_CO2_lme")


#modern respiration (after partitioning based on 2017 atm)
willow17_mod_CO2_lme <-lme(log(RecentFlux17) ~ Mesh*Treatment,
                         data = willow,
                         random = ~ 1|Plot)
plot(resid(willow17_mod_CO2_lme))
qqnorm(resid(willow17_mod_CO2_lme))
willow17_mod_CO2_effects <- summary(willow17_mod_CO2_lme)$tTable
willow17_mod_CO2_df <- round(as.data.frame(willow17_mod_CO2_effects),3)
colnames(willow17_mod_CO2_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(willow17_mod_CO2_df,"willow_recent17_CO2_lme")
```

# Figure 3 DO14C and DOC Conc 
```{r doc plot}

names(doc_long) #check data

# ------------- identify factors instead of characters
doc_long$Mesh <- as.factor(doc_long$Mesh)
doc_long$Treatment <- as.factor(doc_long$Treatment)
doc_long$Plot <- as.factor(doc_long$Plot)

birchdoc <- doc_long[doc_long$Vegetation == "Birch",]
willowdoc <- doc_long[doc_long$Vegetation == "Willow",]


#comparing heterotrophic and root ingrowth differences in DO14C and Conc based on gridling and mesh
#then print

bDO14C <- ggplot(data = birchdoc, aes(x = Mesh, y = DO14C, fill = Vegetation)) +
  geom_jitter(width = 0.05, color = "#006600", size = 4, alpha = 0.5) +  # Jittered points
  facet_wrap(~interaction(Treatment))+
  scale_fill_manual(values=c("#006600")) +
  ylim(70, 90)+
  xlab(" ")+
  ylab("DO14C (‰)")+
  theme_classic()+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  theme(plot.title = element_text(hjust = 0.6),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=12, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"))+
  scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))

wDO14C <- ggplot(data = willowdoc, aes(x = Mesh, y = DO14C, fill = Vegetation)) +
  geom_jitter(width = 0.05, color = "#000066", size = 4, alpha = 0.5) +  # Jittered points
  facet_wrap(~interaction(Treatment))+
  scale_fill_manual(values=c("#000066")) +
  ylim(70, 90)+
  ylab("DO14C (‰)")+
  theme_classic()+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  theme(plot.title = element_text(hjust = 0.6),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=12, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"))+
    scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))


bDOC<- ggplot(data = birchdoc, aes(x = Mesh, y = DOC_Conc, fill = Vegetation)) +
  geom_jitter(width = 0.05, color = "#006600", size = 4, alpha = 0.5) +  # Jittered points
  facet_wrap(~interaction(Treatment))+
  scale_fill_manual(values=c("#006600"))+
  ylim(0, 25)+
  ylab("DOC Concentration (mg L-1)")+
  xlab("Treatment")+
  theme_classic()+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  theme(plot.title = element_text(hjust = 0.6),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=12, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"))+
  scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))

wDOC<- ggplot(data = willowdoc, aes(x = Mesh, y = DOC_Conc, fill = Vegetation)) +
  geom_jitter(width = 0.05, color = "#000066", size = 4, alpha = 0.5) +  # Jittered points
  facet_wrap(~interaction(Treatment))+
  scale_fill_manual(values=c("#000066"))+
  ylim(0, 25)+
  ylab("DOC Concentration (mg L-1)")+
  xlab("Treatment")+
  theme_classic()+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  theme(plot.title = element_text(hjust = 0.6),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=12, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"))+
  scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))


#combine plots
alldoc <- ggarrange(bDOC, bDO14C, wDOC, wDO14C,
                    labels = c("a)", "b)", "c)","d)"),
                    ncol = 2,nrow=2)
alldoc
ggsave(alldoc, filename = "Rplot_DOCconc_DO14C_CtrlGirdled_BirchWillow_Points.svg", device = "svg",height = 8, width = 12, units = "in") 

```

## DOC Statistics
```{r doc stats}

#Birch
#DOC Concentration linear mixed model (Package nlme version 3.1-166)

birchDOC_lme <-lme(log(DOC_Conc) ~ Mesh*Treatment,
                   data = birchdoc,
                   random = ~ 1|Plot)
plot(resid(birchDOC_lme))
qqnorm(resid(birchDOC_lme))
birchDOC_effects <- summary(birchDOC_lme)$tTable
birchDOC_doc <- round(as.data.frame(birchDOC_effects),3)
colnames(birchDOC_doc) <- c("Estimate", "Std. Error", "doc", "t-value", "p-value")
rownames_to_column(birchDOC_doc,"birchDOC_lme")

#DO14C signature linear mixed model (Package nlme version 3.1-166)
birchDO14C_lme <- lme(log(DO14C) ~ Mesh*Treatment,
                      data = birchdoc,
                      random = ~ 1|Plot)
plot(resid(birchDO14C_lme))
qqnorm(resid(birchDO14C_lme))
birchDO14C_effects <- summary(birchDO14C_lme)$tTable
birchDO14C_doc <- round(as.data.frame(birchDO14C_effects),3)
colnames(birchDO14C_doc) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(birchDO14C_doc,"birchDO14C_lme")

#Willow 
#DOC Concentration linear mixed model (Package nlme version 3.1-166)
willowDOC_lme <-lme(log(DOC_Conc) ~ Mesh*Treatment,
                    data = willowdoc,
                    random = ~ 1|Plot)
plot(resid(willowDOC_lme))
qqnorm(resid(willowDOC_lme))
willowDOC_effects <- summary(willowDOC_lme)$tTable
willowDOC_doc <- round(as.data.frame(willowDOC_effects),3)
colnames(willowDOC_doc) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(willowDOC_doc,"willowDOC_lme")

#DO14C signature linear mixed model (Package nlme version 3.1-166)
willowDO14C_lme <- lme(log(DO14C) ~ Mesh*Treatment,
                       data = willowdoc,
                       random = ~ 1|Plot)
summary(willowDO14C_lme)
plot(resid(willowDO14C_lme))
qqnorm(resid(willowDO14C_lme))
willowDO14C_effects <- summary(willowDO14C_lme)$tTable
willowDO14C_doc <- round(as.data.frame(willowDO14C_effects),3)
colnames(willowDO14C_doc) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
rownames_to_column(willowDO14C_doc, "willowDO14C_lme" )

#PRINT TABLE OUTPUT

# Write the combined table to the CSV, appending the results of each model

#xl_lst <- list('birchDOC_doc' = birchDOC_doc, 'birchDO14C_doc' = birchDO14C_doc, 'willowDOC_doc' = willowDOC_doc,'willowDO14C_doc' = willowDO14C_doc)

# Write the data list in a new workbook
#write.xlsx(xl_lst, file = "DOC_lme_workbook.xlsx", rowNames = TRUE)

```

# Figure 4 Respiration Stacked Bargraph 
```{r stacked barplot}
df2 <- mod_df %>%
  group_by(Treatment, Vegetation) %>%
  summarize_at(c("H14C","R14C","Hflux","Rflux","Hflux21_adj","Hflux21R","frecent21R","RecentHFlux21","RecentRFlux21"), funs(mean = mean,n(),sd = sd, se = sd(.)/sqrt(n())))%>%
  ungroup()
df2

summ <- mod_df %>%
  summarize_at(c("H14C","R14C","Hflux","Rflux","Hflux21_adj","Hflux21R","frecent21R","RecentHFlux21","RecentRFlux21"), funs(mean = mean,n(),sd = sd, se = sd(.)/sqrt(n())))%>%
  ungroup()
summ

df2_fluxmean_long <- df2 %>%
  dplyr::select(.,"Vegetation", "Treatment", "RecentRFlux21_mean", "Hflux21R_mean","RecentHFlux21_mean","Hflux21_adj_mean") %>%
  pivot_longer(.,cols = c("RecentRFlux21_mean","Hflux21R_mean","Hflux21_adj_mean","RecentHFlux21_mean"), names_to = "FluxType",values_to = "Flux_mean")

#average total flux over vegetation, treatments and mesh for error bars
flux_total <- df2 %>%
  dplyr::select(.,"Vegetation", "Treatment", "Rflux_mean", "Hflux21R_mean", "Hflux_mean", "Hflux21_adj_mean") %>%
  pivot_longer(.,cols = c("Rflux_mean","Hflux21R_mean", "Hflux_mean","Hflux21_adj_mean"), names_to = "FluxType",values_to = "Flux_mean")

df2_fluxse_long <- df2 %>%
  dplyr::select(.,"Rflux_se","Hflux21R_se", "Hflux_se","Hflux21_adj_se") %>%
  pivot_longer(.,cols = c("Rflux_se","Hflux21R_se", "Hflux_se","Hflux21_adj_se"), names_to = "Flux_seType",values_to = "Flux_se")

df2_long <- cbind(df2_fluxmean_long,df2_fluxse_long)
df2_long$Mesh <- ifelse(df2_long$FluxType == "RecentRFlux21_mean", "Ingrowth", ifelse(df2_long$FluxType == "Hflux21R_mean","Ingrowth", "Exclusion"))
df2_long$TrtMesh <- paste(df2_long$Treatment,df2_long$Mesh)
df2_long$TrtMesh <- factor(df2_long$TrtMesh, c("Control Ingrowth", "Control Exclusion","Girdled Ingrowth", "Girdled Exclusion"))

df2_total_long <- cbind(flux_total,df2_fluxse_long)
df2_total_long$Mesh <- ifelse(df2_total_long$FluxType == "Rflux_mean", "Ingrowth", ifelse(df2_total_long$FluxType == "Hflux21R_mean","Ingrowth", "Exclusion"))
df2_total_long$TrtMesh <- paste(df2_total_long$Treatment,df2_total_long$Mesh)
df2_total_long$TrtMesh <- as.factor(df2_total_long$TrtMesh)
df2_total_long$TrtMesh <- fct_relevel(df2_total_long$TrtMesh, "Control Ingrowth", "Control Exclusion", "Girdled Ingrowth", "Girdled Exclusion")

# Stacked Barplot
StackResp <- ggplot(data = df2_long, aes(x = TrtMesh, y = Flux_mean, fill = FluxType)) +
  facet_grid(.~ Vegetation)+
  geom_bar(position = position_stack(reverse = TRUE), stat= "identity") +
  geom_errorbar(data = df2_total_long, aes(x = TrtMesh, ymin = Flux_mean-Flux_se, ymax = Flux_mean+Flux_se),width=0.3,size=0.2)+
  scale_fill_manual(values=c("grey","grey","grey","black","black","black"),labels = c("Rflux_mean" = "","Hflux21_adj_mean" = "Ancient Flux","Hflux21R_mean" = "Ancient Flux","Hflux21R_mean" = "Ancient Flux", "RecentRFlux21_mean" = "Recent Flux", "Hflux_mean" = "", "RecentHFlux21_mean" = "Recent Flux"))+ # blue
  #scale_color_manual(values=c("black"))
  ylab("umol CO2 hr-1") +
  xlab(" ") +
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.5),
        #legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=20, color = "black"),
        axis.text.x= element_text(size=14, color = "black",angle =70, hjust = 1),
        axis.text.y = element_text(size = 16, color = "black"))+
  scale_y_continuous(limits = c(0, 5), expand = c(0,0))

StackResp
 #ggsave(StackResp, filename = "Rplot_StackedResp_CtrlGirdled_Birch&Willow.svg", #device = "svg",
    # height = 8, width = 10, units = "in") 

```

# Figure 5 Rhizosphere Priming Effects 
```{r priming effects}

 #Compiling the ratios of rooted vs unrooted mesh bags to get a priming effect ratio (positive priming > 1)
 #Statistics should be analyzed with the natural log, since mean values will be overestimated otherwise
 
 priming_effect <- mod_df %>%
   mutate(pe = Hflux21R / Hflux21_adj)%>%
  mutate(pe_b = Hflux21Rb / Hflux21b_adj)%>%
  mutate(pe17 = Hflux17R / Hflux17_adj)%>%
  mutate(pe17_b = Hflux17Rb / Hflux17b_adj)%>%
 as.data.frame()
 
 priming_effect$Treatment <- as.factor(priming_effect$Treatment)
 priming_effect$Vegetation <- as.factor(priming_effect$Vegetation)
 priming_effect$Plot <- as.factor(priming_effect$Plot)
priming_effect_noBC6 <- priming_effect[priming_effect$Plot != "BC6",]
 
 birch_pe <- priming_effect[priming_effect$Vegetation == "Birch",]
 willow_pe <- priming_effect[priming_effect$Vegetation == "Willow",]
 
 bcpe <- priming_effect[priming_effect$Vegetation == "Birch" & priming_effect$Treatment == "Control" ,]
 wcpe <- priming_effect[priming_effect$Vegetation == "Willow" & priming_effect$Treatment == "Control" ,]
 bgpe <- priming_effect[priming_effect$Vegetation == "Birch" & priming_effect$Treatment == "Girdled" ,]
 wgpe <- priming_effect[priming_effect$Vegetation == "Willow" & priming_effect$Treatment == "Girdled" ,]
 
 #overall mean
 exp(mean(log(priming_effect$pe))) #mean with 14CO2 estimation and 2021 air
  exp(mean(log(priming_effect$pe_b))) #mean with bulk 14C estimation and 2021 air
   exp(mean(log(priming_effect$pe17))) #mean with 14CO2 estimation and 2017 air
   exp(mean(log(priming_effect$pe17_b))) #mean with bulk 14C estimation and 2017 air
 
 # setting up means for the graph for RPE
 exp(mean(log(bcpe$pe))) #mean birch control
 exp(mean(log(wcpe$pe))) #mean willow control
 exp(mean(log(bgpe$pe))) #mean birch girdled
 exp(mean(log(wgpe$pe))) #mean willow girdled
 
 RPE_points <- ggplot(data = priming_effect, aes(x = interaction(Vegetation,Treatment), y = pe, fill = Vegetation, color = Vegetation)) +
   geom_jitter(shape = 21,width = 0.05, size = 4, alpha = 0.5, stroke = 1.1)+  # Jittered points
   scale_fill_manual(values=c("#006600","#000066")) +
   scale_color_manual(values=c("#006600","#000066")) +
   ylab("Rhizosphere Priming Effect Ratio")+
   xlab(" ")+
   geom_hline(yintercept = 1, color = "grey", linetype = "dashed", linewidth = 1.0)+
   theme_classic()+
   theme(plot.title = element_text(hjust = 0.6),
         legend.position = ("none"),
         axis.line = element_line(size = .5, linetype = "solid"), 
         axis.title=element_text(size=16, color = "black"),
         axis.text.x= element_text(size=14, color = "black"),
         axis.text.y = element_text(size = 16, color = "black"))+
   scale_x_discrete(labels=c("Birch.Control" = "Birch Control",
                             "Birch.Girdled" = "Birch Girdled",
                             "Willow.Control" = "Willow Control", 
                             "Willow.Girdled" = "Willow Girdled"))+
   stat_summary(fun = function(x) exp(mean(log(x))), 
                geom = "crossbar", 
                aes(ymin = after_stat(y), ymax = after_stat(y)), 
                width = 0.4, color = "black")+
   scale_y_continuous(breaks = seq(0, 6, by = 1))
  RPE_points
  
# ggsave(RPE_points, filename = "Figure 5_RPE_points_raw.svg", device = "svg",
  #       height = 5, width = 7, units = "in")

  
 # TESTING RHIZOSPHERE PRIMING EFFECTS

 # test significance of RPEs - with random effects
 #root_rpe <-lme(log(pe) ~ Vegetation*Treatment,
 #               data = priming_effect,
 #               random = ~ 1|Plot)
 # summary(root_rpe) # sd of random effects tiny compared to sd of residuals so drop 
 # intervals(root_rpe)
 # plot(resid(root_rpe))
 # qqnorm(resid(root_rpe))
 # rpe_root_effects <- summary(root_rpe)$tTable
 # rpe_root_df <- round(as.data.frame(rpe_root_effects),3)
 # colnames(rpe_root_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
 # rownames_to_column(rpe_root_df,"rpe_root_lme")
 
 # Fit model again with fixed effects only:
 # change contrasts
 options(contrasts = c("contr.sum", "contr.poly"))
 # change back to default
 #options(contrasts = c("contr.treatment", "contr.poly"))
 root_rpe <-lm(log(pe) ~ Vegetation*Treatment,
                data = priming_effect)
 rpe_root_effects <- summary(root_rpe)$coefficients[, 1:4]
 rpe_root_df <- round(as.data.frame(rpe_root_effects),3)
 colnames(rpe_root_df) <- c("Estimate", "Std. Error", "t-value", "p-value")
 rownames_to_column(rpe_root_df,"rpe_root_lm")
 
 # confidence intervals
 confint(root_rpe)["(Intercept)", ]
 exp(confint(root_rpe)["(Intercept)", ])
 
root_rpe17 <-lm(log(pe17) ~ Vegetation*Treatment,
                data = priming_effect)
 rpe_root17_effects <- summary(root_rpe17)$coefficients[, 1:4]
 rpe_root17_df <- round(as.data.frame(rpe_root17_effects),3)
 colnames(rpe_root17_df) <- c("Estimate", "Std. Error", "t-value", "p-value")
 rownames_to_column(rpe_root17_df,"rpe_root17_lm")
 
 # confidence intervals
 confint(root_rpe)["(Intercept)", ]
 exp(confint(root_rpe)["(Intercept)", ])
 
 #test intercept without predictors
 summary(lm(log(pe) ~ 1,
            data = priming_effect))
 
 # without BC6 - random effects
 root_rpe_noB6C <-lme(log(pe) ~ Vegetation*Treatment,
                      data = priming_effect_noBC6,
                      random = ~ 1|Plot)
 summary(root_rpe_noB6C ) # sd of random effects tiny compared to sd of residuals so drop 
 # intervals(root_rpe_noB6C)
 # plot(resid(root_rpe_noB6C ))
 # qqnorm(resid(root_rpe_noB6C ))
 # rpe_root_effects_noB6C <- summary(root_rpe_noB6C )$tTable
 # rpe_root_noB6C_df <- round(as.data.frame(rpe_root_effects_noB6C),3)
 # colnames(rpe_root_noB6C_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
 # rownames_to_column(rpe_root_noB6C_df,"rpe_root_noB6C_lme")
 
 # fixed effects only - no BC6
 root_rpe_noB6C <-lm(log(pe) ~ Vegetation*Treatment,
                     data = priming_effect_noBC6)
 rpe_root_effects_noB6C <- summary(root_rpe_noB6C)$coefficients[, 1:4]
 rpe_root_df_noB6C <- round(as.data.frame(rpe_root_effects_noB6C),3)
 colnames(rpe_root_df_noB6C) <- c("Estimate", "Std. Error", "t-value", "p-value")
 rownames_to_column(rpe_root_df_noB6C,"rpe_root_lm_noB6C")
 
 # confidence intervals
 confint(root_rpe_noB6C)["(Intercept)", ]
 exp(confint(root_rpe_noB6C)["(Intercept)", ])
 
 #test intercept without predictors
 summary(lm(log(pe) ~ 1,
            data = priming_effect_noBC6))
 
 
# xl_lst <- list('rpe_root_df' = rpe_root_df, 'rpe_root_df_noB6C' = rpe_root_df_noB6C)
 
 # Write the data list in a new workbook
# write.xlsx(xl_lst, file = "rpe_lme_workbook.xlsx", rowNames = TRUE)
 
```

# Supplementary Information
```{r supp info}

## Table A1 and Figure A1 Root Respiration Sensitivity Test
library(plotrix)
library(segmented)

box<-read.table("root box expt.csv",header=TRUE,sep=",")
box$plot=as.factor(box$plot)
#turn slope into flux (um s-1)
box$flux=box$slope*((101.3*(184/1000000))/(0.008*283))

plot(box$flux~box$fineroots)#clear relationship between fine root biomass and flux
plot(box$flux~box$courseroots)#NO relationship between coarse root biomass and flux

#remove boxes with dead roots
box=box[box$plot!= "1",]
box=box[box$plot!= "13",]
box=box[box$plot!= "18",]
box=box[box$plot!= "10",]
box=box[box$plot!= "12",]
box$plot=droplevels(box$plot)
#turn slope into flux (um s-1)
box$flux=box$slope*((101.3*(184/1000000))/(0.008*283))
#vol 2*8*11.5
#normalise by g fine root
box$fluxFR=box$flux/box$fineroots

#normalise by the measurement 30 mins before cut
aveFR=data.frame(plot=levels(box$plot),aveFRflux=as.numeric(tapply(box$fluxFR[box$hour==-0.5],box$plot[box$hour==0.5],mean,na.rm=T)))

#merge datasets in order to get normalised rate
box=merge(aveFR,box,by='plot')
box$fluxFRnorm=box$fluxFR/box$aveFRflux

#make line graph of normalised fluxes 
#remove first flux measurement
cflux=tapply(box$fluxFRnorm[box$treatment=='c'&box$hour>-4],box$hour[box$treatment=='c'&box$hour>-4],mean,na.rm=T)
sflux=tapply(box$fluxFRnorm[box$treatment=='s'&box$hour>-4],box$hour[box$treatment=='s'&box$hour>-4],mean,na.rm=T)


cfluxSE=tapply(box$fluxFRnorm[box$treatment=='c'&box$hour>-4],box$hour[box$treatment=='c'&box$hour>-4],std.error,na.rm=T)
sfluxSE=tapply(box$fluxFRnorm[box$treatment=='s'&box$hour>-4],box$hour[box$treatment=='s'&box$hour>-4],std.error,na.rm=T)

hour=as.numeric(levels(as.factor(box$hour[box$hour>-4])))
hour_comp=c(hour[1:11],9,11,13,15,17,19)

#plot compressed x axis
par(mfrow=c(1,1),mar=c(5,5,1,1))
plotCI(hour_comp,y=cflux,cfluxSE,xlim=c(-4,20),ylim=c(0.4,1.8),sfrac=0.003,xaxt='n',yaxt='n',xlab='',ylab='',pch=16,las=2,cex=2,col='black')
par(new=T)
plotCI(hour_comp+0.05,y=sflux,sfluxSE,xlim=c(-4,20),ylim=c(0.4,1.8),sfrac=0.003,xaxt='n',yaxt='n',xlab='',ylab='',pch=1,las=2,cex=2,col='red')
abline(v=0,lty=2,col='red')
abline(h=1,lty=2,col='black')
axis(side=1,at=hour_comp,labels=hour,cex.axis=2)
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8),labels=c(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8),cex.axis=2,las=2)

#calculate s/c from average values
plot((sflux/cflux)~hour)

#segment regression analysis to find inflection point
severed.mod=lm(fluxFRnorm~hour,data=box[box$treatment=='s'&box$hour>-0.5,])
severed.seg=segmented(severed.mod,seg.z=hour,psi=48)
summary(severed.seg)
anova(severed.seg)
draw.history(severed.seg)
severed.seg

control.mod=lm(fluxFRnorm~hour,data=box[box$treatment=='c'&box$hour>-0.5,])
control.seg=segmented(control.mod,seg.z=hour,psi=6)
summary(control.seg)
anova(control.seg)
control.seg

#uncompressed x axis
par(mfrow=c(1,1),mar=c(7,7,1,1))
plotCI(hour,y=cflux,cfluxSE,xlim=c(-4,170),ylim=c(0.4,1.8),sfrac=0.003,xaxt='n',yaxt='n',xlab='',ylab='',pch=16,las=2,cex=2,col='black')
par(new=T)
plotCI(hour+0.05,y=sflux,sfluxSE,xlim=c(-4,170),ylim=c(0.4,1.8),sfrac=0.003,xaxt='n',yaxt='n',xlab='',ylab='',pch=1,las=2,cex=2,col='red')
abline(v=0,lty=2,col='red')
abline(h=1,lty=2,col='black')
axis(side=1,at=hour,labels=hour,cex.axis=2)
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8),labels=c(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8),cex.axis=2,las=2)
mtext(side=1,line=4,adj=0.5, 'Time from severing (h)',cex=2.5)
mtext(side=2,line=5,adj=0.5, 'Normalised root respiration',cex=2.5)
lines(severed.seg)

#ggsave(., filename = ".svg", device = "svg",
#       height = 5, width = 7, units = "in")

## Figure B1 Keeling Plot with Regression

df_long$recip.CO2_volume_ml <- 1/df_long$CO2_volume_ml

mesh.2mm <- subset(df_long, Mesh == "2mm")
mesh.1um <- subset(df_long, Mesh == "1um")

coef.2mm <- coef(lm(d13C.VPDB ~ recip.CO2_volume_ml, data = mesh.2mm))
R2.2mm <- summary(lm(d13C.VPDB ~ recip.CO2_volume_ml, data = mesh.2mm))$r.squared
coef.1um <- coef(lm(d13C.VPDB ~ recip.CO2_volume_ml, data = mesh.1um))
R2.1um <- summary(lm(d13C.VPDB ~ recip.CO2_volume_ml, data = mesh.1um))$r.squared

co2vol_d13co2_mesh <- ggplot(data = df_long, aes(x = recip.CO2_volume_ml, y = d13C.VPDB, color = Mesh)) +
  geom_point(size = 3) +
  stat_smooth(method = "lm", se = FALSE, fullrange = TRUE) +
  scale_color_manual(values=c("grey","black")) + 
  ylab(expression(phantom()^{13}*CO[2]*" (‰)")) +
  xlim(c(0,0.4))+
  ylim(c(-32,-12))+
  xlab(expression(1/CO[2]*"volume (ml)")) +
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.6),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=20, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        axis.text.y = element_text(size = 16, color = "black"))+
  annotate("text", x=0.1, y=-13, label="y = 12.4x - 25.5, R2 = 0.56", size = 5, color = "grey")+
  annotate("text", x=0.1, y=-15, label="y = 26.3x - 28.8, R2 = 0.41", size = 5, color = "black")
co2vol_d13co2_mesh

ggsave(co2vol_d13co2_mesh, filename = "Figure_B1_KeelingPlot.svg", device = "svg",
      height = 4, width = 6, units = "in")


## Figure B2 Root Biomass in Birch and Willow 

mod_df_long$root_Cg_m2 <- mod_df_long$Root.biomass_g*.5/(64*pi)*10000 
root.biomass <- mod_df_long %>%
 group_by(Mesh,Vegetation,Treatment) %>%
  dplyr::summarise(average_value = mean(root_Cg_m2, na.rm = TRUE), .groups = "drop")

birch_girdled <- birch %>% 
  filter(Treatment == "Girdled" & Mesh == "2mm")

broot <- ggplot(data = birch, aes(x = Mesh, y = Root.biomass_g, fill = Vegetation)) +
  #geom_boxplot() +
  geom_jitter(width = 0.05, color = "#006600", size = 4, alpha = 0.6) +  # Jittered points
  facet_wrap(~interaction(Treatment), ncol=2)+
  scale_fill_manual(values=c("#006600"))+ # blue
  ylab("Root biomass (g)") +
  xlab(" ")+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  theme_classic()+
  scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))+
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=20, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        axis.text.y = element_text(size = 16, color = "black"))+
  scale_y_continuous(limits = c(0, 8))

mean(willow$Root.biomass_g)
sd(willow$Root.biomass_g)/sqrt(length(willow$Root.biomass_g))

willow_control <- willow %>% 
  filter(Treatment == "Control" & Mesh == "2mm")

willow_girdled <- willow %>% 
  filter(Treatment == "Girdled" & Mesh == "2mm")

#willow root biomass g C / m2
(mean(willow_girdled$Root.biomass_g)/(2*64*pi))*10000
(sd(willow_girdled$Root.biomass_g)/sqrt(length(willow_girdled$Root.biomass_g)))/(2*64*pi)*10000

(mean(willow_control$Root.biomass_g)/(2*64*pi))*10000
(sd(willow_control$Root.biomass_g)/sqrt(length(willow_control$Root.biomass_g)))/(2*64*pi)*10000

wroot <- ggplot(data = willow, aes(x = Mesh, y = Root.biomass_g, fill = Vegetation)) +
  #geom_boxplot() +
  geom_jitter(width = 0.05, color = "#000066", size = 4, alpha = 0.6)+
  facet_wrap(~interaction(Treatment), ncol=2)+
  scale_fill_manual(values=c("#000066"))+ #blue
  ylab("Root biomass (g)")+
  xlab(" ")+
  stat_summary(fun = mean, geom = "crossbar", width = 0.5, color = "black")+
  theme_classic()+
  scale_x_discrete(limits = c("2mm", "1um"), labels=c("1um" = "Exclusion (1-μm)", "2mm" = "Ingrowth (2-mm)"))+
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=20, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        axis.text.y = element_text(size = 16, color = "black"))+
  scale_y_continuous(limits = c(0, 8))

#save file to edit in powerpoint first change y axis title

all_root <-  ggarrange(broot, wroot,
                       labels = c("a)", "b)"),
                       ncol = 1, nrow = 2)
all_root

ggsave(all_root, filename = "Figure_B2_rootbiomass.svg", device = "svg",
      height = 4, width = 6, units = "in")


## Root biomass STATS linear mixed model
#(Package nlme version 3.1-166)
root_lme <-lme(log(Root.biomass_g+1) ~ Mesh*Treatment,
                     data = mod_df_long,
                     random = ~ 1|Plot)
 plot(resid(root_lme))
 qqnorm(resid(root_lme))
root_effects <- summary(root_lme)$tTable
root_df <- round(as.data.frame(root_effects),3)
 colnames(root_df) <- c("Estimate", "Std. Error", "DF", "t-value", "p-value")
 rownames_to_column(root_df,"root_lme")

##  pMC and CO2 flux regression

df_long$group <- paste(df_long$Vegetation, df_long$Mesh, sep = "_")
Flux_14CO2 <- ggplot(data = df_long, aes(x = Flux, y = X14CO2, fill = group, color = Vegetation)) +
  #geom_boxplot()+
  geom_point(shape = 21, size = 4, alpha = 0.5, stroke = 1.1)+  # Jittered points
  scale_fill_manual(values=c("Birch_2mm" = "#006600",
                             "Birch_1um" = "white",
                             "Willow_2mm" = "#000066", 
                             "Willow_1um" = "white"))+ 
  scale_color_manual(values=c("Birch" = "#006600",
                              "Willow" = "#000066")) + 
  ylab(expression(phantom()^{14}*CO[2]*" (%modern)")) +
  xlab(expression(CO[2]*" flux "*"(umol "* m^-2 * s^-1*")")) +
  ylim(82,100)+
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.6),
        legend.position = ("none"),
        axis.line = element_line(size = .5, linetype = "solid"), 
        axis.title=element_text(size=20, color = "black"),
        axis.text.x= element_text(size=14, color = "black"),
        axis.text.y = element_text(size = 16, color = "black"))
Flux_14CO2
#ggsave(Flux_14CO2, filename = "Figure A3_raw.svg", device = "svg",
#       height = 5, width = 7, units = "in")


 
```