### Load packages ####
library(ggplot2)
library(readxl)
library(nlme)
library(ggpubr)
library(dplyr)
library(vegan)
library(tidyr)
library(emmeans)

# set preferences
theme_set(theme_bw())

### load data ####
DataBFS <- read.delim("CarryOverNeffects_PlotData.txt")
DataBFS$CCage <- DataBFS$SamplingYear - DataBFS$year_cut
DataBFS$Plot <- as.factor(DataBFS$Plot)
DataBFS$Clearcut <- as.factor(DataBFS$Clearcut)
DataBFS$Pair <- as.factor(DataBFS$Pair)
DataBFS$Area <- as.factor(DataBFS$Area)
DataBFS$SamplingYear <- as.factor(DataBFS$SamplingYear)
DataBFS$CN <- DataBFS$TotC_./DataBFS$TotN_.
DataBFS$Fertilized <- factor(DataBFS$Fertilized, labels = c("Unfertilized", "Fertilized"))
DataBFS$SoilCstock <- as.numeric(DataBFS$SoilCstock)
DataBFS$SoilNstock <- as.numeric(DataBFS$SoilNstock)
DataBFS$SOM <- as.numeric(DataBFS$SOM)
DataBFS$Moisture <- as.numeric(DataBFS$Moisture)

SoilRespBFS <- read.delim("CarryOverNeffects_SoilRespirationData.txt")
SoilRespBFS$Plot <- as.factor(SoilRespBFS$Plot)
SoilRespBFS$SamplingPoint <- as.factor(SoilRespBFS$SamplingPoint)
SoilRespBFS <- left_join(SoilRespBFS, DataBFS[,c(2,8,9,13,15:27)], by = "Plot")
SoilRespBFS$OrganicLayer <- as.numeric(SoilRespBFS$OrganicLayer)
SoilRespBFS$Moisture <- as.numeric(SoilRespBFS$Moisture)
SoilRespBFS$SOM <- as.numeric(SoilRespBFS$SOM)

TreeGrowth <- read.delim("CarryOverNeffects_TreeGrowthRateData.txt")
TreeGrowth$GrowthRate <- as.numeric(TreeGrowth$GrowthRate)*2
TreeGrowth$Plot <- as.factor(TreeGrowth$Plot)
TreeHeight <- read.delim("CarryOverNeffects_TreeHeightData.txt")
TreeHeight$Plot <- as.factor(TreeHeight$Plot)
Vegetation <- read.delim("CarryOverNeffects_VegetationPlotData.txt")
Vegetation$Plot <- as.factor(Vegetation$Plot)

OrgLayer <- SoilRespBFS %>% group_by(Plot) %>% summarise(OrgLayer=mean(OrganicLayer))
DataBFS <- left_join(DataBFS,OrgLayer,by="Plot")

Data1yo <- DataBFS %>% filter(SamplingYear == 2023)
Data4to13yo <- DataBFS %>% filter(SamplingYear == 2022)
Vegetation <- left_join(Vegetation, Data4to13yo, by="Plot")

### Soil variables ####
NSmod <- lme(SoilNstock ~ CCage*Fertilized, random = ~ 1|Pair/Clearcut, data = Data4to13yo, na.action = na.exclude)
summary(NSmod)
anova(NSmod)

CNmod <- lme(CN ~ CCage*Fertilized, random = ~ 1|Pair/Clearcut, data = Data4to13yo, na.action = na.exclude)
summary(CNmod)
anova(CNmod)

pHmod <- lme(pH ~ CCage*Fertilized, random = ~ 1|Pair/Clearcut, data = Data4to13yo, na.action = na.exclude)
summary(pHmod)
anova(pHmod)

OLmod <- lme(OrgLayer ~ CCage*Fertilized, random = ~ 1|Pair/Clearcut, data = Data4to13yo, na.action = na.exclude)
summary(OLmod)
anova(OLmod)

CSmod <- lme(SoilCstock ~ CCage*Fertilized, random = ~ 1|Pair/Clearcut, data = Data4to13yo, na.action = na.exclude)
summary(CSmod)
anova(CSmod)


## Soil respiration 2022 sampling ####
SRmod <- lm(flux.quad ~ poly(SoilTemp,2, raw=TRUE), data=SoilRespBFS)
summary(SRmod)
SoilRespBFS$SResp_cor <- residuals(SRmod) - 628.528 + 142.6836*mean(SoilRespBFS$SoilTemp) - 4.5611*(mean(SoilRespBFS$SoilTemp)^2)

SoilRespBFS2 <- aggregate(SoilRespBFS$SResp_cor, by=list(SoilRespBFS$Pair, SoilRespBFS$Clearcut, SoilRespBFS$Plot, SoilRespBFS$Round), FUN=mean, na.rm=TRUE)
colnames(SoilRespBFS2) <- c("Pair", "Clearcut", "Plot", "Round", "MeanSR")
SoilRespBFS2 <- left_join(SoilRespBFS2,DataBFS[,c(2,8,26)], by="Plot")
SoilRespBFS2$Pair <- as.factor(SoilRespBFS2$Pair)
SoilRespBFS2$Clearcut <- as.factor(SoilRespBFS2$Clearcut)
SoilResp1yo <- SoilRespBFS2 %>% filter(CCage < 2)
SoilResp4to13yo <- SoilRespBFS2 %>% filter(CCage > 2)

SRmod2 <- lme(MeanSR ~ Round*CCage + Fertilized, random = ~ 1|Pair/Clearcut, data = SoilResp4to13yo, na.action = na.exclude)
summary(SRmod2)
anova(SRmod2)

SRmod3 <- lme(MeanSR ~ Fertilized, random = ~ 1|Pair/Clearcut, data = SoilResp1yo, na.action = na.exclude)
summary(SRmod3)
anova(SRmod3)

# overview plot soil variables & soil respiration ####
CSplot <- ggplot(data = Data4to13yo, aes(x=Fertilized, y=SoilCstock)) + geom_boxplot() + geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y=bquote("Organic layer C stock (t "~ha^-1~")")) + theme(legend.position = "none") +
  coord_cartesian(ylim=c(0,97))
Orgplot <- ggplot(data = Data4to13yo, aes(x=Fertilized, y=OrgLayer)) + geom_boxplot() + geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
 labs(x=NULL, y="Depth organic layer (cm)") + theme(legend.position = "none") +
  coord_cartesian(ylim=c(0,14))
CNplot <- ggplot(data=Data4to13yo, aes(x=Fertilized,y=CN)) + geom_boxplot() + geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
 labs(x=NULL, y = "C:N organic layer") + theme(legend.position = "none") +
  coord_cartesian(ylim=c(0,40))
NSplot <- ggplot(data=Data4to13yo, aes(x=Fertilized,y=SoilNstock)) + geom_boxplot() + geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + 
  labs(x=NULL, y = bquote("Organic layer N stock (t "~ha^-1~")")) + theme(legend.position = "none") +
  coord_cartesian(ylim=c(0,3.1))
pHplot <- ggplot(data=Data4to13yo, aes(x=Fertilized,y=pH)) + geom_boxplot() + geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y = "pH organic layer") + theme(legend.position = "none") +
  coord_cartesian(ylim=c(0,5))
CRplot <- ggplot(data = SoilRespBFS2, aes(x=Fertilized, y=MeanSR)) + geom_boxplot() +   geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y=bquote("Soil "~CO[2]~" efflux (mg C"~m^-2~ha^-1~")")) + facet_wrap(~Round) + theme(legend.position = "none")
ggarrange(CSplot,NSplot, CNplot,pHplot,Orgplot,CRplot, nrow = 2,ncol = 3, labels = "AUTO")

### Ground vegetation ####
# vegetation cover
aggregate(Vegetation$Mosses/1.5, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Vegetation$Mosses/1.5, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)
aggregate(Vegetation$Lichens/1.5, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Vegetation$Lichens/1.5, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)
aggregate(Vegetation$VascularPlants/1.5, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Vegetation$VascularPlants/1.5, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)

Vegetation2 <- Vegetation %>% filter(VascularPlants>60)
VPmod <- lme(sqrt(VascularPlants/1.5) ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Vegetation)
anova(VPmod)

Mmod <- lme(sqrt(Mosses/1.5) ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Vegetation)
anova(Mmod)

Lmod <- lme(sqrt(Lichens/1.5) ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Vegetation)
anova(Lmod)

# diversity
VascPl <- Vegetation[,17:95]
Data4to13yo$VP_S <- specnumber(VascPl)
Data4to13yo$VP_H <- diversity(VascPl, index="shannon")
Data4to13yo$VP_PE <- Data4to13yo$VP_H/log(Data4to13yo$VP_S) # Pielou's evenness

aggregate(Data4to13yo$VP_S, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Data4to13yo$VP_S, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)
aggregate(Data4to13yo$VP_H, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Data4to13yo$VP_H, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)
aggregate(Data4to13yo$VP_PE, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Data4to13yo$VP_PE, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)

Smod <- lme(VP_S ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Data4to13yo)
anova(Smod)

Hmod <- lme(VP_H ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Data4to13yo)
anova(Hmod)

PEmod <- lme(VP_PE ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Data4to13yo)
anova(PEmod)

VPplot <- ggplot(data = Vegetation, aes(x=Data4to13yo$Fertilized, y=VascularPlants/1.5)) + geom_boxplot() + 
  geom_jitter(aes(col=Data4to13yo$Fertilized), width=0.1, alpha=0.5) +   scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="Cover (%)", title = "Vascular plant cover") + theme(legend.position = "none", plot.title=element_text(size=10)) +
  coord_cartesian(ylim = c(0,100))
Mplot <- ggplot(data = Vegetation, aes(x=Data4to13yo$Fertilized, y=Mosses/1.5)) + geom_boxplot() + 
  geom_jitter(aes(col=Data4to13yo$Fertilized), width=0.1, alpha=0.5) +   scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y=NULL, title = "Moss cover") + theme(legend.position = "none", plot.title=element_text(size=10)) + 
  coord_cartesian(ylim = c(0,100))
Lplot <- ggplot(data = Vegetation, aes(x=Data4to13yo$Fertilized, y=Lichens/1.5)) + geom_boxplot() + 
  geom_jitter(aes(col=Data4to13yo$Fertilized), width=0.1, alpha=0.5) +   scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y=NULL, title = "Lichen cover") + theme(legend.position = "none", plot.title=element_text(size=10)) + 
  coord_cartesian(ylim = c(0,100))
Splot <- ggplot(data = Data4to13yo, aes(x=Fertilized, y=VP_S)) + geom_boxplot() +
  geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) +   scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y=" ", title = "Vascular plant species richness") + theme(legend.position = "none", plot.title=element_text(size=10)) + 
  coord_cartesian(ylim = c(0,25))
Hplot <- ggplot(data = Data4to13yo, aes(x=Fertilized, y=VP_H)) + geom_boxplot() + 
  geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) +   scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y=NULL, title = "Vascular plant Shannon diversity") + theme(legend.position = "none", plot.title=element_text(size=10)) + 
  coord_cartesian(ylim = c(0,2.5))
PEplot <- ggplot(data = Data4to13yo, aes(x=Fertilized, y=VP_PE)) + geom_boxplot() + 
  geom_jitter(aes(col = Fertilized), width=0.1, alpha=0.5) +   scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y=NULL, title = "Vascular plant Pielou's evenness") + theme(legend.position = "none", plot.title=element_text(size=10)) + 
  coord_cartesian(ylim = c(0,0.9))
ggarrange(VPplot, Mplot,Lplot,Splot,Hplot,PEplot, nrow=2, ncol=3, labels = "AUTO")

# Community composition
# NMDS
VP_nmds <- metaMDS(VascPl, distance="bray")
en = envfit(VP_nmds, Data4to13yo[,8:15], permutations = 999, na.rm = TRUE)

data.scores <- as.data.frame(VP_nmds$points)  #Using the scores function from vegan to extract the site scores and convert to a data.frame
data.scores$Plot <- Data4to13yo$Plot  # create a column of site names
data.scores$Fertilized <- Data4to13yo$Fertilized  
species.scores <- as.data.frame(scores(VP_nmds, "species"))  #Using the scores function from vegan to extract the species scores and convert to a data.frame
species.scores$species <- rownames(species.scores)  # create a column of species, from the rownames of species.scores
en_coord_cont = as.data.frame(scores(en, "vectors")) * ordiArrowMul(en)
en_coord_cat = as.data.frame(scores(en, "factors")) * ordiArrowMul(en)
plot(VP_nmds)
# prepare data to add centroids per clearcut to the NMDS
OrdEl <- ordiellipse(VP_nmds, groups = Data4to13yo$Clearcut)
list_names <- names(OrdEl)
# Create a matrix to store the extracted numbers
extracted_numbers <- matrix(0, nrow = length(OrdEl), ncol = 2)

# Loop through the list of lists and extract the numbers
for (i in 1:length(OrdEl)) {
  sublist <- OrdEl[[i]]
  extracted_numbers[i, ] <- sublist$center
}
# Convert the matrix to a DataFrame
centroids <- data.frame(Clearcut = list_names, NMDS1 = extracted_numbers[, 1], NMDS2 = extracted_numbers[, 2])
centroids <- left_join(centroids, Data4to13yo, by="Clearcut")

ggplot() + 
  geom_point(data=data.scores,aes(x=MDS1,y=MDS2,colour=Fertilized),shape=8,size=1) + # add the point markers
  geom_point(data=centroids,aes(x=NMDS1,y=NMDS2,colour=Fertilized),size=5) + # add the point markers
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +  
  geom_text(aes(x=0.95, y=1.1, label=paste("stress = ", format(VP_nmds$stress, digits=3), sep="")), size=5) +
  coord_equal() + labs(x = "NMDS 1", y = "NMDS2") + theme(legend.title = element_blank(), text = element_text(size=16), legend.position = "top")

VPdist <- vegdist(x=VascPl,method = "bray")
adonis2(VPdist ~ Fertilized*CCage, permutations = 10000, method = "bray", data = Data4to13yo, strata=Data4to13yo$Pair)
adonis2(VascPl ~ Fertilized*CCage, permutations = 999, method = "bray", data = Data4to13yo)

### cover of ericoids and grasses ####
aggregate(Vegetation$Vaccinium_myrtilis/1.5, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Vegetation$Vaccinium_myrtilis/1.5, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)
aggregate(Vegetation$Vaccinium_vitisidaea/1.5, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Vegetation$Vaccinium_vitisidaea/1.5, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)
aggregate(Vegetation$Calluna_vulgaris/1.5, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Vegetation$Calluna_vulgaris/1.5, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)
aggregate(Vegetation$Deschampsia_flexuosa/1.5, by=list(Data4to13yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(Vegetation$Deschampsia_flexuosa/1.5, by=list(Data4to13yo$Fertilized), FUN=sd, na.rm=TRUE)

Vaccinium_vitisidaea <- Vegetation$Vaccinium_vitisidaea
VacVImod <- lme(Vaccinium_vitisidaea ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Data4to13yo)
anova(VacVImod)

Vaccinium_myrtilis <- Vegetation$Vaccinium_myrtilis
VacMmod <- lme(Vaccinium_myrtilis ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Data4to13yo)
anova(VacMmod)

Calluna_vulgaris <- Vegetation$Calluna_vulgaris
CalVmod <- lme(sqrt(Calluna_vulgaris) ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Data4to13yo)
anova(CalVmod)

Deschampsia_flexuosa <- Vegetation$Deschampsia_flexuosa
DesFmod <- lme(Deschampsia_flexuosa ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = Data4to13yo)
anova(DesFmod)

### Tree growth & height ####
TreeGrowth2 <- left_join(TreeGrowth,Data4to13yo, by="Plot")
TreeGrowth3 <- aggregate(TreeGrowth2$GrowthRate, by=list(TreeGrowth2$Species, TreeGrowth2$Plot), FUN=mean, na.rm=TRUE)
colnames(TreeGrowth3) <- c("Species", "Plot", "MeanGR")
TreeGrowth3 <- left_join(TreeGrowth3,Data4to13yo[,c(2,3,4,8)], by="Plot")
TreeGrowth3 <- unique(TreeGrowth3)
TGmod2 <-lme(MeanGR ~ Species + Fertilized, random = ~  1|Pair/Clearcut, data=TreeGrowth3)
summary(TGmod2)
anova(TGmod2)
ggplot(data= TreeGrowth3, aes(x=Fertilized, y=MeanGR)) + geom_boxplot() + geom_jitter(aes(col=Fertilized), width=0.2, alpha=0.5) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(legend.position = "none") +
  facet_wrap(~Species) + labs(x=NULL, y="Average yearly diameter increase (mm)")


PaGrowth <- TreeGrowth2 %>% filter(Species == "Picea abies")
PaGmod <-lme(GrowthRate ~ Fertilized*CCage, random = ~  1|Pair/Clearcut, data=PaGrowth)
summary(PaGmod)
anova(PaGmod)

TreeHeight2 <- left_join(TreeHeight,Data4to13yo, by="Plot") %>% 
  filter(Species=="Betula pendula"|Species=="Betula pubescens"|Species=="Picea abies"|
           Species=="Pinus sylvestris")
TreeHeight3 <- aggregate(TreeHeight2$height, by=list(TreeHeight2$Species, TreeHeight2$Plot), FUN=mean, na.rm=TRUE)
colnames(TreeHeight3) <- c("Species", "Plot", "MeanHeight")
TreeHeight3 <- left_join(TreeHeight3,Data4to13yo[,c(2,3,4,8,26)], by="Plot")

THmod <-lme(log(MeanHeight) ~ Species*CCage + Fertilized, random = ~ 1|Pair/Clearcut, data=TreeHeight3)
summary(THmod)
anova(THmod)

ggplot(data = TreeHeight3, aes(x=CCage, y=MeanHeight, col=Fertilized, fill=Fertilized)) + geom_jitter(width = 0.1, alpha=0.5) + facet_wrap(~Species) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + scale_fill_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x="Time since clearcutting (years)", y="Average tree height (cm)") + geom_smooth(method = "lm", formula = 'y ~ x') + theme(legend.title = element_blank(),strip.text = element_text(face = "italic"))
aggregate(TreeHeight3$MeanHeight, by=list(TreeHeight3$Species, TreeHeight3$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(TreeHeight3$MeanHeight, by=list(TreeHeight3$Species, TreeHeight3$Fertilized), FUN=sd, na.rm=TRUE)

# tree density, diversity and composition
TreeDensity <- TreeHeight %>% 
  #filter(Species=="Betula pendula"|Species=="Betula pubescens"|Species=="Picea abies"|Species=="Pinus sylvestris"|Species=="Populus tremula"|Species=="Salix caprea"|Species=="Sorbus aucuparia") %>% 
  count(Plot,Species)
TreeDensity$n_m2 <- TreeDensity$n/28.2743338823
TreeDensity$n_ha <- TreeDensity$n_m2*10000

TreeComp <- pivot_wider(TreeDensity[,c(1,2,4)], names_from = Species, values_from = n_m2, values_fill = 0)
TreeComp2 <- left_join(TreeComp, Data4to13yo, by="Plot")
TreeComp2$Tree_S <- specnumber(TreeComp[,2:17])
TreeComp2$Tree_dens <- rowSums(TreeComp[,2:17])

trmod1 <- lme(Tree_S ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data=TreeComp2)
anova(trmod1)

trmod2 <- lme(Tree_dens ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data=TreeComp2)
anova(trmod2)

### 2023 sampling - C and N stocks ####
SoilResp1yo <- filter(SoilRespBFS, SamplingYear=="2023")
SoilResp1yo <- filter(SoilResp1yo, SoilResp1yo$OrganicLayer != "NA")

SRmod <- lm(flux.quad ~ poly(SoilTemp,2, raw=TRUE), data=SoilResp1yo)
summary(SRmod)
SoilResp1yo$SResp_cor <- residuals(SRmod) + 4039.686 - 578.792*mean(SoilResp1yo$SoilTemp) + 23.327*(mean(SoilResp1yo$SoilTemp)^2)

SRmod2 <- lme(SResp_cor ~ Fertilized, random =  ~ 1|Pair/Clearcut/Plot, data = SoilResp1yo)
summary(SRmod2)
anova(SRmod2)

aggregate(SoilResp1yo$SResp_cor, by=list(SoilResp1yo$Fertilized), FUN=mean, na.rm=TRUE)
aggregate(SoilResp1yo$SResp_cor, by=list(SoilResp1yo$Fertilized), FUN=sd, na.rm=TRUE)

SCSmod <- lme(log(SoilCstock) ~ Fertilized, random = ~ 1|Pair/Clearcut, data = Data1yo)
summary(SCSmod)
anova(SCSmod)
emmeans(SCSmod,specs = pairwise ~ Fertilized, type="response")

SNSmod <- lme(SoilNstock ~ Fertilized, random = ~ 1|Pair/Clearcut, data = Data1yo)
summary(SNSmod)
anova(SNSmod)
emmeans(SNSmod,specs = pairwise ~ Fertilized, type="response")

CNmod <- lme(CN ~ Fertilized, random = ~ 1|Pair/Clearcut, data = Data1yo)
summary(CNmod)
anova(CNmod)

OLmod <- lme(OrgLayer ~ Fertilized, random = ~ 1|Pair/Clearcut, data = Data1yo)
summary(OLmod)
anova(OLmod)

### Fungal communities - both sampling campaigns ####
data <- read.delim("CarryOverNeffects_OTUtable.txt")
fungi <- data %>% filter(Kingdom == "k__Fungi")
EcMF <- data %>% filter(Lifestyle == "ectomycorrhizal")
SapAsc <- data %>% filter(Lifestyle == "saprotrophic ascomycetes")
RAasc <- data %>% filter(Lifestyle == "root-associated ascomycete")
Sap_ag <- data %>% filter(Lifestyle == "saprotrophic agaricomycetes")
MyPa <- data %>% filter(Lifestyle == "mycoparasite")

Sample <- colnames(data[2:145])
fungireads <- data.frame(colSums(fungi[2:145]))
TotReads <- data.frame(colSums(data[2:145]))
EcMFreads <- data.frame(colSums(EcMF[2:145]))
SapAscreads <- data.frame(colSums(SapAsc[2:145]))
RAascreads <- data.frame(colSums(RAasc[2:145]))
Sap_agreads <- data.frame(colSums(Sap_ag[2:145]))
MyPareads <- data.frame(colSums(MyPa[2:145]))

# ITS copy nr quantifications ####
# correct the qPCR ITS quantifications to exclude plant reads & change it to copynr per g OM (multiply by SOM ratio)
PercFungReads <- fungireads$colSums.fungi.2.145../TotReads$colSums.data.2.145..
DataBFS$ITScopyNR_OM <- DataBFS$ITScopyNR*DataBFS$SOM
ITScopyNRcor <- data.frame(DataBFS$ITScopyNR_OM*PercFungReads)
DataBFS$ITScopyNRcor <- ITScopyNRcor$DataBFS.ITScopyNR_OM...PercFungReads

PercEcMFReads <- EcMFreads$colSums.EcMF.2.145../TotReads$colSums.data.2.145..
PercSapAscReads <- SapAscreads$colSums.SapAsc.2.145../TotReads$colSums.data.2.145..
PercRAascReads <- RAascreads$colSums.RAasc.2.145../TotReads$colSums.data.2.145..
PercSap_agReads <- Sap_agreads$colSums.Sap_ag.2.145../TotReads$colSums.data.2.145..
PercMyPaReads <- MyPareads$colSums.MyPa.2.145../TotReads$colSums.data.2.145..

ITScopyNR_EcMF <- data.frame(DataBFS$ITScopyNR_OM*PercEcMFReads)
DataBFS$ITScopyNR_EcMF <- ITScopyNR_EcMF$DataBFS.ITScopyNR_OM...PercEcMFReads
ITScopyNR_SapAsc <- data.frame(DataBFS$ITScopyNR_OM*PercSapAscReads)
DataBFS$ITScopyNR_SapAsc <- ITScopyNR_SapAsc$DataBFS.ITScopyNR_OM...PercSapAscReads
ITScopyNR_RAasc <- data.frame(DataBFS$ITScopyNR_OM*PercRAascReads)
DataBFS$ITScopyNR_RAasc <- ITScopyNR_RAasc$DataBFS.ITScopyNR_OM...PercRAascReads
ITScopyNR_Sap_ag <- data.frame(DataBFS$ITScopyNR_OM*PercSap_agReads)
DataBFS$ITScopyNR_Sap_ag <- ITScopyNR_Sap_ag$DataBFS.ITScopyNR_OM...PercSap_agReads
ITScopyNR_MyPa <- data.frame(DataBFS$ITScopyNR_OM*PercMyPaReads)
DataBFS$ITScopyNR_MyPa <- ITScopyNR_MyPa$DataBFS.ITScopyNR_OM...PercMyPaReads

Data1yo <- DataBFS %>% filter(SamplingYear == 2023)
Data4to13yo <- DataBFS %>% filter(SamplingYear == 2022)

# test for differences with fertilization
ITSmod1 <- lme(sqrt(ITScopyNRcor) ~ Fertilized*CCage, random= ~ 1|Pair/Clearcut, data = DataBFS, na.action = na.exclude)
summary(ITSmod1)
ITSmod1_output <- anova(ITSmod1)
EMM <- emmeans(ITSmod1, "Fertilized")

ITSmod3 <- lme(sqrt(ITScopyNR_EcMF) ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = DataBFS, na.action = na.exclude)
summary(ITSmod3)
ITSmod3_output <- anova(ITSmod3)

ITSmod4 <- lme(sqrt(ITScopyNR_SapAsc) ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = DataBFS, na.action = na.exclude)
summary(ITSmod4)
ITSmod4_output <- anova(ITSmod4)

ITSmod5 <- lme(sqrt(ITScopyNR_RAasc) ~ Fertilized*CCage,  random = ~ 1|Pair/Clearcut, data = DataBFS, na.action = na.exclude)
summary(ITSmod5)
ITSmod5_output <- anova(ITSmod5)

ITSmod6 <- lme(sqrt(ITScopyNR_Sap_ag) ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = DataBFS, na.action = na.exclude)
summary(ITSmod6)
ITSmod6_output <- anova(ITSmod6)

ITSmod7 <- lme(sqrt(ITScopyNR_MyPa) ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data = DataBFS, na.action = na.exclude)
summary(ITSmod7)
ITSmod7_output <- anova(ITSmod7)


# plots
p0 <- ggplot(DataBFS, aes(x=Fertilized, y=ITScopyNRcor, col=Fertilized)) + geom_boxplot() +  geom_jitter(width = 0.1, alpha = 0.2) +
  geom_text(x=1.5, y=1550,label=paste("F = ",format(ITSmod1_output[3,3], digits=3), ", p = ",format(ITSmod1_output[3,4], digits=2),sep=""), col = "black") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=8), legend.position = "none") +
  labs(x=NULL, y=bquote("Fungal ITS2 copy numbers m"~g^-1~"OM"), title = "All fungi")
p1 <- ggplot(DataBFS, aes(x=Fertilized, y=ITScopyNR_SapAsc, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1, alpha = 0.2) +
  geom_text(x=1.5, y=1550,label=paste("F = ",format(ITSmod4_output[2,3], digits=3), ", p = ",format(ITSmod4_output[2,4], digits=2),sep=""), col = "black") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=8), legend.position = "none") +
  labs(x=NULL, y=NULL, title = "Saprotrophic ascomycetes")
p2 <- ggplot(DataBFS, aes(x=Fertilized, y=ITScopyNR_RAasc, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1, alpha = 0.2) +
  geom_text(x=1.5, y=1550,label=paste("F = ",format(ITSmod5_output[2,3], digits=3), ", p = ",format(ITSmod5_output[2,4], digits=2),sep=""), col = "black") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=8), legend.position = "none") +
  labs(x=NULL, y=NULL, title = "Root-associated ascomycetes")
p3 <- ggplot(DataBFS, aes(x=Fertilized, y=ITScopyNR_Sap_ag, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1, alpha = 0.2) +
  geom_text(x=1.5, y=1550,label=paste("F = ",format(ITSmod6_output[2,3], digits=3), ", p = ",format(ITSmod6_output[2,4], digits=2),sep=""), col = "black") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=8), legend.position = "none") +
  labs(x=NULL, y=NULL, title = "Saprotrophic agaricomycetes")
p4 <- ggplot(DataBFS, aes(x=Fertilized, y=ITScopyNR_EcMF, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1, alpha = 0.2) +
  geom_text(x=1.5, y=1550,label=paste("F = ",format(ITSmod3_output[2,3], digits=3), ", p = ",format(ITSmod3_output[2,4], digits=2),sep=""), col = "black") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=8), legend.position = "none") +
  labs(x=NULL, y=NULL, title = "Ectomycorrhizal fungi")
ggarrange(p0,p1,p2,p3,p4,ncol=5, widths = c(1.1,1,1,1,1))


# Diversity metrics ####
OTUtab <- data.frame(t(fungi[2:145]))
colnames(OTUtab) <- fungi$SampleID
SeqDepth <- data.frame(colSums(fungi[2:145]))
DataBFS$SeqDepth <- SeqDepth$colSums.fungi.2.145..

OTUtab <- filter(OTUtab, DataBFS$SeqDepth > 500)
Data <- filter(DataBFS, DataBFS$SeqDepth > 500)
Fungi_ratios <- OTUtab/DataBFS$SeqDepth
Fungi_CNr <- Fungi_ratios*Data$ITScopyNRcor

Fungi_CNr <- Fungi_CNr %>% filter(Fungi_CNr$scata6378_0 != "NA")
tOTUtab <- data.frame(t(Fungi_CNr))
taxa <- data.frame(fungi[c(1,146:156)])
colnames(taxa)[1] <- c("OTU_ID")
tOTUtab <- filter(tOTUtab, taxa$SumOfSeqs > 10) 
taxasub <- filter(taxa, taxa$SumOfSeqs > 10) 
OTUtab <- data.frame(t(tOTUtab))

S_OTU <- data.frame(specnumber(OTUtab))
Data <- Data %>% filter(ITScopyNRcor != "NA")
Data$S_OTU <- S_OTU$specnumber.OTUtab.
Shannon <- data.frame(diversity(OTUtab, index="shannon"))
Data$Shannon <- Shannon$diversity.OTUtab..index....shannon..
Data$OTU_PE <- Data$Shannon/log(Data$S_OTU)

# OTU richness
Smod2 <- lme(S_OTU ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data=Data)
summary(Smod2)
Smod_output <- anova(Smod2)

# Shannon diversity
Shmod <- lme(Shannon ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data=Data)
summary(Shmod)
Shmod_output <- anova(Shmod)

# Pielou's evenness
PEmod <- lme(OTU_PE ~ Fertilized*CCage, random = ~ 1|Pair/Clearcut, data=Data)
summary(PEmod)
PEmod_output <- anova(PEmod)


# plots fertilization effect
Splot <- ggplot(Data, aes(x=Fertilized, y=S_OTU, col=Fertilized)) + geom_boxplot() + geom_jitter(width=0.2, alpha = 0.5) + labs(x=NULL, y =NULL, title="OTU richness") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=11), legend.title = element_blank())
Shplot <- ggplot(Data, aes(x=Fertilized, y=Shannon, col=Fertilized)) + geom_boxplot() + geom_jitter(width=0.2, alpha = 0.5) + labs(x=NULL, y =NULL, title="Shannon diversity") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=11), legend.title = element_blank())
PEplot <- ggplot(Data, aes(x=Fertilized, y=OTU_PE, col=Fertilized)) + geom_boxplot() + geom_jitter(width=0.2, alpha = 0.5) + labs(x=NULL, y =NULL, title="Pielou's evenness") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=11), legend.title = element_blank())
ggarrange(Splot,Shplot,PEplot, ncol=3, labels="AUTO", common.legend = TRUE)

# Community composition Fungi ####
Fguilds <- data.frame(aggregate(tOTUtab, by=list(Guild=taxasub$Lifestyle), FUN=sum, na.action=NULL), row.names = 1)
Fguilds$Total <- rowSums(Fguilds)
Fguilds <- Fguilds[order(Fguilds$Total),]
tFguilds <- data.frame(t(Fguilds[,1:138]))

Fguilds_fert <- aggregate(tFguilds, by=list(Fertilized = Data$Fertilized), FUN=mean)
Fguilds_bp3 <- pivot_longer(Fguilds_fert, cols = 2:11, names_to = "Lifestyle",
                            values_to = "Sequences", values_drop_na = FALSE)
Fguilds_bp3 <- Fguilds_bp3[order(Fguilds_bp3$Sequences),]
Fguilds_bp3$Lifestyle <- factor(Fguilds_bp3$Lifestyle, 
                                        levels = c("Other","Unknown","lichenized","mycoparasite","plant.pathogen","other.saprotrophs", "saprotrophic.agaricomycetes", "ectomycorrhizal","root.associated.ascomycete", "saprotrophic.ascomycetes"), 
                                        labels = c("Other", "Unknown", "Lichenized", "Mycoparasites", "Plant pathogens","Other saprotrophs", "Saprotrophic agaricomycetes", "Ectomycorrhizal fungi", "Root-associated ascomycetes",  "Saprotrophic ascomycetes"))
Guilds_fert_plot <- ggplot(Fguilds_bp3, aes(x=Fertilized, y=Sequences, fill=Lifestyle)) + geom_bar(stat = "identity") + 
  labs(x=NULL, y=bquote("Fungal ITS2 copy numbers m"~g^-1~"dry soil")) + theme(legend.title= element_blank()) + 
  scale_fill_manual(values=c("#333333ff","#b3b3b3ff", "#5f8dd3ff", "#e9afafff", "#aa4400ff", "#ff6600ff", "#ff9955ff", "#ffd42aff", "#aa8800ff", "#ffeeaaff", "#554400ff"))

Fclass <- data.frame(aggregate(t(OTUtab), by=list(Class=taxasub$Class), FUN=sum), row.names = 1)
Fclass$Total <- rowSums(Fclass)
Fclass <- Fclass[order(Fclass$Total, decreasing = FALSE),]
tFclass <- data.frame(t(Fclass[,1:138]))
tFclass$Other <- rowSums(tFclass[,1:57])
tFclass2 <- tFclass[,58:68]

Fclass_fert <- aggregate(tFclass2, by=list(Fertilized = Data$Fertilized), FUN=mean)
Fclass_bp3 <- pivot_longer(Fclass_fert, cols = 2:12, names_to = "Class",
                           values_to = "Sequences", values_drop_na = FALSE)
Fclass_bp3$Class <- factor(Fclass_bp3$Class, 
                           levels = c("Other","c__Lecanoromycetes", "c__Umbelopsidomycetes","c__Tremellomycetes","c__Mortierellomycetes","c__Sordariomycetes","c__Dothideomycetes", "c__Archaeorhizomycetes","c__Eurotiomycetes","c__Agaricomycetes", "c__Leotiomycetes"), 
                           labels = c("Other", "Lecanoromycetes", "Umbelopsidomycetes", "Tremellomycetes", "Mortierellomycetes", "Sordariomycetes", "Dothideomycetes",  "Archaeorhizomycetes", "Eurotiomycetes","Agaricomycetes", "Leotiomycetes"))
Class_fert_plot <- ggplot(Fclass_bp3, aes(x=Fertilized, y=Sequences, fill=Class)) + geom_bar(stat = "identity") + 
  labs(x=NULL, y=bquote("Fungal ITS2 copy numbers m"~g^-1~"dry soil")) + theme(legend.title= element_blank()) + 
  scale_fill_manual(values=c("#b3b3b3ff","#351431","#8E5572",  "#4F4789", "#775253", "#BDC696", "#636940", "#D1D3C4", "#201335", "#628395", "#474A2C"))
ggarrange(Guilds_fert_plot,Class_fert_plot, ncol=2, labels="AUTO", widths = c(1.05,0.95))


## NMDS Fungal communities ####
Fungi_nmds <- metaMDS(OTUtab)
stressplot(Fungi_nmds)
plot(Fungi_nmds)

data.scores <- as.data.frame(Fungi_nmds$points)  #Using the scores function from vegan to extract the site scores and convert to a data.frame
data.species <- as.data.frame(Fungi_nmds$species)
data.species$OTU_ID <- row.names(data.species)
data.species <- left_join(data.species, taxasub, by = "OTU_ID")
EcM.species <- data.species %>% filter(primary_lifestyle == "ectomycorrhizal")
Sapr.species <- data.species %>% filter(primary_lifestyle == "soil_saprotroph"|primary_lifestyle == "litter_saprotroph"|primary_lifestyle == "wood_saprotroph"|primary_lifestyle == "unspecified_saprotroph")
# prepare data to add centroids per clearcut to the NMDS
OrdEl <- ordiellipse(Fungi_nmds, groups = Data$Clearcut)
list_names <- names(OrdEl)
# Create a matrix to store the extracted numbers
extracted_numbers <- matrix(0, nrow = length(OrdEl), ncol = 2)

# Loop through the list of lists and extract the numbers
for (i in 1:length(OrdEl)) {
  sublist <- OrdEl[[i]]
  extracted_numbers[i, ] <- sublist$center
}
# Convert the matrix to a DataFrame
centroids <- data.frame(Clearcut = list_names, NMDS1 = extracted_numbers[, 1], NMDS2 = extracted_numbers[, 2])
centroids <- left_join(centroids, Data, by="Clearcut")

ggplot() +   #geom_point(data=data.species, aes(x=MDS1,y=MDS2), shape=4, color="grey",alpha =0.3) + 
  geom_point(data=data.scores,aes(x=MDS1,y=MDS2,col=Data$Fertilized), shape=8,size=1) +
  geom_point(data=centroids,aes(x=NMDS1,y=NMDS2,col=Fertilized, shape=SamplingYear),size=5) +
  geom_text(aes(x=0.7, y=0.85, label=paste("stress = ", format(Fungi_nmds$stress, digits=3), sep="")), size = 5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + theme(text = element_text(size=16), legend.position = "top") +
  labs(x = "NMDS 1", y = "NMDS2", color=NULL, shape =NULL)
ggplot() +   #geom_point(data=data.species, aes(x=MDS1,y=MDS2), shape=4, color="grey",alpha =0.3) + 
  geom_point(data=data.scores,aes(x=MDS1,y=MDS2,col=Data$CCage, shape=Data$Fertilized), size=2, alpha=0.5) +
  geom_point(data=centroids,aes(x=NMDS1,y=NMDS2,col=CCage, shape=Fertilized),size=5) +
  geom_text(aes(x=0.7, y=0.95, label=paste("stress = ", format(Fungi_nmds$stress, digits=3), sep="")), size = 4) +
  theme(text = element_text(size=12),) +
  labs(x = "NMDS 1", y = "NMDS2", color="Time since", shape =NULL)

# Permanova Fungal communities ####
Fungidist <- vegdist(x=OTUtab,method = "bray")

mod1 <- betadisper(Fungidist,Data$Fertilized)
mod2 <- betadisper(Fungidist,Data$Pair)
anova(mod1)

adonis2(Fungidist ~ Fertilized, permutations = 10000, by="margin",
        method = "bray", data = Data, strata=Data$Pair)
adonis2(Fungidist ~ Fertilized + CCage, permutations = 10000, by="margin",
        method = "bray", data = Data)

OTUtab1yo <- OTUtab %>% filter(Data$SamplingYear == 2023)
Fungidist1yo <- vegdist(x=OTUtab1yo,method = "bray")
adonis2(Fungidist1yo ~ Fertilized, permutations = 10000, by="margin",
        method = "bray", data = Data1yo, strata=Data1yo$Pair)

FungiDist2 <- vegdist(x = centroids[,2:3], method = "euclidian")
adonis2(FungiDist2 ~ Fertilized + CCage, permutations = 10000, by="margin",
        method = "bray", data = Data, strata=Data$Pair)

## Ectomycorrhizal fungi ####
EcMF <- data.frame(t(OTUtab)) %>% filter(taxasub$Lifestyle=="ectomycorrhizal")
EcMFtax <- taxasub %>% filter(primary_lifestyle=="ectomycorrhizal")

# calculate diversity indices & ITS copy nrs
S_EcMF <- data.frame(specnumber(t(EcMF)))
Data$S_EcMF <- S_EcMF$specnumber.t.EcMF..

H_EcMF <- data.frame(diversity(t(EcMF), "shannon"))
Data$H_EcMF <- H_EcMF$diversity.t.EcMF....shannon..
PE_EcMF <- H_EcMF/log(Data$S_EcMF)
Data$PE_EcMF <- PE_EcMF$diversity.t.EcMF....shannon..

EcMmod1 <- lme(sqrt(S_EcMF) ~ Fertilized + poly(CCage,2), random= ~ 1|Pair/Clearcut, data = Data)
summary(EcMmod1)
anova(EcMmod1)

EcM_plot1 <- ggplot(Data, aes(x=Fertilized, y=S_EcMF, col=Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + 
  labs(x=NULL, y="OTU richness", title="Ectomycorrhizal fungi")

EcMFmod2 <- lme(H_EcMF ~ Fertilized*CCage, random= ~ 1|Pair/Clearcut, Data)
summary(EcMFmod2)
EcMFmod_output2 <- anova(EcMFmod2)

EcMF_plot2 <- ggplot(Data, aes(x=Fertilized, y=H_EcMF, col=Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + 
  labs(x=NULL, y="Shannon diversity", title="Ectomycorrhizal fungi")

EcMFmod3 <- lme(PE_EcMF ~ Fertilized*CCage, random= ~ 1|Pair/Clearcut, Data)
summary(EcMFmod3)
EcMFmod_output3 <- anova(EcMFmod3)

EcMF_plot3 <- ggplot(Data, aes(x=Fertilized, y=PE_EcMF, col=Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) + 
  labs(x=NULL, y="Pielou's evenness", title="Ectomycorrhizal fungi")
ggarrange(EcM_plot1, EcMF_plot2, EcMF_plot3, ncol=3, labels = "AUTO", common.legend = TRUE)

EcMFdist <- vegdist(t(EcMF), method="bray")
adonis2(EcMFdist~ Fertilized, permutations = 999, method = "bray", data = Data, strata=Data$Pair)


## Saprotrophic agaricomycetes ####
Sapr_ag <- data.frame(t(OTUtab)) %>% filter(taxasub$Lifestyle=="saprotrophic agaricomycetes")
Saprtax_ag <- taxasub %>% filter(taxasub$Lifestyle=="saprotrophic agaricomycetes")

S_Sapr_ag <- data.frame(specnumber(t(Sapr_ag)))
H_Sapr_ag <- data.frame(diversity(t(Sapr_ag), "shannon"))
Data$S_Sapr_ag <- S_Sapr_ag$specnumber.t.Sapr_ag..
Data$H_Sapr_ag <- H_Sapr_ag$diversity.t.Sapr_ag....shannon..
PE_Sapr_ag <- data.frame(H_Sapr_ag/log(Data$S_Sapr_ag))
Data$PE_Sapr_ag <- PE_Sapr_ag$diversity.t.Sapr_ag....shannon..

aggregate(Data$S_Sapr_ag, by = list(Data$Fertilized), FUN = mean)
aggregate(Data$S_Sapr_ag, by = list(Data$Fertilized), FUN = sd)
aggregate(Data$H_Sapr_ag, by = list(Data$Fertilized), FUN = mean)
aggregate(Data$H_Sapr_ag, by = list(Data$Fertilized), FUN = sd)
aggregate(Data$PE_Sapr_ag, by = list(Data$Fertilized), FUN = mean)
aggregate(Data$PE_Sapr_ag, by = list(Data$Fertilized), FUN = sd)

Sapr_agmod1 <- lme(sqrt(S_Sapr_ag) ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(Sapr_agmod1)
anova(Sapr_agmod1)

Sapr_ag_plot <- ggplot(Data, aes(x=Fertilized, y=S_Sapr_ag, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.5) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="OTU richness", title="Saprotrophic Agaricomycetes")


Sapr_agmod2 <- lme(H_Sapr_ag ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(Sapr_agmod2)
Sapr_agmod_output2 <- anova(Sapr_agmod2)

Sapr_ag_plot2 <- ggplot(Data, aes(x=Fertilized, y=H_Sapr_ag, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + facet_wrap(~SamplingYear) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="Shannon diversity", title="Saprotrophic Agaricomycetes")

Sapr_agmod3 <- lme(PE_Sapr_ag ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(Sapr_agmod3)
Sapr_agmod_output3 <- anova(Sapr_agmod3)

Sapr_ag_plot3 <- ggplot(Data, aes(x=Fertilized, y=PE_Sapr_ag, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + facet_wrap(~SamplingYear) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="Pielou's evenness", title="Saprotrophic Agaricomycetes")
ggarrange(Sapr_ag_plot, Sapr_ag_plot2, Sapr_ag_plot3, ncol = 3, labels="AUTO", common.legend = TRUE)

Sapr_ag_dist <- vegdist(t(Sapr_ag),"bray")
adonis2(Sapr_ag_dist~ Fertilized + CCage, permutations = 999, method = "bray", data = Data, strata=Data$Pair)

## Root-associated ascomycetes ####
RAasc <- data.frame(t(OTUtab)) %>% filter(taxasub$Lifestyle=="root-associated ascomycete")
RAasc_tax <- taxasub %>% filter(taxasub$Lifestyle=="root-associated ascomycete")

S_RAasc <- data.frame(specnumber(t(RAasc)))
H_RAasc <- data.frame(diversity(t(RAasc), "shannon"))
Data$S_RAasc <- S_RAasc$specnumber.t.RAasc..
Data$H_RAasc <- H_RAasc$diversity.t.RAasc....shannon..
Data$PE_RAasc <- Data$H_RAasc/log(Data$S_RAasc)

RAascmod1 <- lme(sqrt(S_RAasc) ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(RAascmod1)
anova(RAascmod1)

RAasc_plot <- ggplot(Data, aes(x=Fertilized, y=sqrt(S_RAasc), col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="OTU richness", title="Root-associated ascomycetes")

RAascmod2 <- lme(H_RAasc ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(RAascmod2)
anova(RAascmod2)

RAasc_plot2 <- ggplot(Data, aes(x=Fertilized, y=H_RAasc, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="Shannon diversity", title="Root-associated ascomycetes")

RAascmod3 <- lme(PE_RAasc ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(RAascmod3)
anova(RAascmod3)

RAasc_plot3 <- ggplot(Data, aes(x=Fertilized, y=PE_RAasc, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="Pielou's evenness", title="Root-associated ascomycetes")
ggarrange(RAasc_plot, RAasc_plot2, RAasc_plot3, ncol = 3, labels="AUTO", common.legend = TRUE)

RAasc_dist <- vegdist(t(RAasc),"bray")
adonis2(RAasc_dist~ Fertilized + CCage, permutations = 999, method = "bray", data = Data, strata=Data$Pair)
adonis2(RAasc_dist~ Fertilized + CCage, permutations = 999, method = "bray", data = Data)

## Saprotrophic ascomycetes ####
SapAsc <- data.frame(t(OTUtab)) %>% filter(taxasub$Lifestyle=="saprotrophic ascomycetes")
SapAsc_tax <- taxasub %>% filter(taxasub$Lifestyle=="saprotrophic ascomycetes")

S_SapAsc <- data.frame(specnumber(t(SapAsc)))
H_SapAsc <- data.frame(diversity(t(SapAsc), "shannon"))
Data$S_SapAsc <- S_SapAsc$specnumber.t.SapAsc..
Data$H_SapAsc <- H_SapAsc$diversity.t.SapAsc....shannon..
Data$PE_SapAsc <- Data$H_SapAsc/log(Data$S_SapAsc)

SapAscmod1 <- lme(S_SapAsc ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(SapAscmod1)
anova(SapAscmod1)

SapAsc_plot <- ggplot(Data, aes(x=Fertilized, y=sqrt(S_SapAsc), col=Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="OTU richness", title="Saprotrophic ascomycetes")

SapAscmod2 <- lme(H_SapAsc ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(SapAscmod2)
anova(SapAscmod2)

SapAsc_plot2 <- ggplot(Data, aes(x=Fertilized, y=H_SapAsc, col=Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="Shannon diversity", title="Root-associated ascomycetes")

SapAscmod3 <- lme(PE_SapAsc ~ Fertilized*CCage, random = ~ 1|Pair, Data)
summary(SapAscmod3)
anova(SapAscmod3)

SapAsc_plot3 <- ggplot(Data, aes(x=Fertilized, y=PE_SapAsc, col=Fertilized)) + geom_boxplot() + geom_jitter(width = 0.2, alpha=0.2) + 
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y="Pielou's evenness", title="Root-associated ascomycetes")
ggarrange(SapAsc_plot, SapAsc_plot2, SapAsc_plot3, ncol = 3, labels="AUTO", common.legend = TRUE)

SapAsc_dist <- vegdist(t(SapAsc),"bray")
adonis2(SapAsc_dist~ Fertilized + CCage, permutations = 999, method = "bray", data = Data, strata=Data$Pair)
adonis2(SapAsc_dist~ Fertilized + CCage, permutations = 999, method = "bray", data = Data)

## Enzymatic activities ####
HydEnz <- read.delim("CarryOverNeffects_EnzymeAssayData.txt")
HydEnz$fert <- factor(HydEnz$fert, labels = c("Unfertilized", "Fertilized"))

HydEnz2 <- aggregate(HydEnz[,c(5:9)], by=list(HydEnz$Plot), FUN=mean, na.rm=TRUE)
colnames(HydEnz2) <- c("Plot", "CBH", "BG", "NAG", "aP", "BXD")
HydEnz2 <- left_join(HydEnz2, Data1yo, by = "Plot")

HydEnz$CBH_tr <- log(HydEnz$CBH_SOM)
HydEnz$BG_tr <- log(HydEnz$BG_SOM)
HydEnz$NAG_tr <- log(HydEnz$NAG_SOM)
HydEnz$aP_tr <- log(HydEnz$aP_SOM)
HydEnz$BXD_tr <- log(HydEnz$BXD_SOM)

CBHmod <- lme(CBH ~ Fertilized, random = ~ 1|Pair/Clearcut, HydEnz2, na.action = na.exclude)
summary(CBHmod)
anova(CBHmod)

BGmod <- lme(BG ~ Fertilized, random = ~ 1|Pair/Clearcut, HydEnz2, na.action = na.exclude)
summary(BGmod)
anova(BGmod)

NAGmod <- lme(NAG ~ Fertilized, random = ~ 1|Pair/Clearcut, HydEnz2, na.action = na.exclude)
summary(NAGmod)
anova(NAGmod)

aPmod <- lme(aP ~ Fertilized, random = ~ 1|Pair/Clearcut, HydEnz2, na.action = na.exclude)
summary(aPmod)
anova(aPmod)

BXDmod <- lme(BXD ~ Fertilized, random = ~ 1|Pair/Clearcut, HydEnz2, na.action = na.exclude)
summary(BXDmod)
anova(BXDmod)


p1 <- ggplot(HydEnz2, aes(x=Fertilized, y=CBH, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1,alpha=0.5, size=2) + 
  theme(plot.title = element_text(size=11), legend.position = "none") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y = bquote(" Enzyme activity (?mol "~g^-1~"SOM"~min^-1~")"), title = "Cellobiohydrolase")
p2 <- ggplot(HydEnz2, aes(x=Fertilized, y=BG, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1,alpha=0.5, size=2) + 
  theme(plot.title = element_text(size=11), legend.position = "none") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y = NULL, title = expression(~beta*"-glucosidase"))
p3 <- ggplot(HydEnz2, aes(x=Fertilized, y=NAG, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1,alpha=0.5, size=2) + 
  theme(plot.title = element_text(size=11), legend.position = "none") + coord_cartesian(ylim = c(0,410)) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y = bquote(" Enzyme activity (?mol "~g^-1~"SOM"~min^-1~")"), title = expression("N-acetyl-"~beta*"-D-glucosaminidase"))
p4 <- ggplot(HydEnz2, aes(x=Fertilized, y=aP, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1,alpha=0.5, size=2) + 
  theme(plot.title = element_text(size=11), legend.position = "none") + coord_cartesian(ylim = c(0,2000)) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y = NULL, title = "Acid phosphatase")
p5 <- ggplot(HydEnz2, aes(x=Fertilized, y=BXD, col = Fertilized)) + geom_boxplot() + geom_jitter(width = 0.1,alpha=0.5, size=2) + 
  theme(plot.title = element_text(size=11), legend.position = "none") +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  labs(x=NULL, y = NULL, title = expression(~beta*"-xylosidase"))
ggarrange(p1,p5,p2,p3,p4, widths = c(1.05, 0.95, 0.95, 1.05, 0.95,1.05), nrow = 2,ncol = 3, labels = "AUTO")


# Figures for in the manuscript ####
Data_mean <- aggregate(Data4to13yo[,c(11:17,27:28)], by=list(Clearcut = Data4to13yo$Clearcut), FUN=mean)
Data_mean <- left_join(Data_mean, Data4to13yo[,c(3,4,8)], by = "Clearcut")
Data_mean <- unique(Data_mean)
Data_wide <- pivot_wider(Data_mean[,2:12], names_from = Fertilized, values_from = c(TotC_., TotN_., pH, Moisture, SOM, SoilCstock, SoilNstock, OrgLayer,CN))

Data_mean2 <- aggregate(Data1yo[,c(11:17,27:28)], by=list(Clearcut = Data1yo$Clearcut), FUN=mean)
Data_mean2 <- left_join(Data_mean2, Data1yo[,c(3,4,8)], by = "Clearcut")
Data_mean2 <- unique(Data_mean2)
Data_wide2 <- pivot_wider(Data_mean2[,2:12], names_from = Fertilized, values_from = c(TotC_., TotN_., pH, Moisture, SOM, SoilCstock, SoilNstock, OrgLayer,CN))

DataNNA <- Data4to13yo %>% filter(!is.na(SoilCstock))
Data_mean3 <- aggregate(DataNNA[,c(11:17,27:28)], by=list(Fertilized = DataNNA$Fertilized), FUN=mean)
Data_mean4 <- aggregate(Data1yo[,c(11:17,27:28)], by=list(Fertilized = Data1yo$Fertilized), FUN=mean)
Data_sd1 <- aggregate(DataNNA[,c(11:17,27:28)], by=list(Fertilized = DataNNA$Fertilized), FUN=sd)
Data_sd2 <- aggregate(Data1yo[,c(11:17,27:28)], by=list(Fertilized = Data1yo$Fertilized), FUN=sd)

ggplot(Data_wide, aes(x=1, xend=2, y=SoilCstock_Unfertilized, yend=SoilCstock_Fertilized)) + geom_segment()
ggplot(Data_wide2, aes(x=1, xend=2, y=SoilCstock_Unfertilized, yend=SoilCstock_Fertilized)) + geom_segment()


PR_upper1 <- Data_mean4$SoilCstock + Data_sd2$SoilCstock
PR_lower1 <- Data_mean4$SoilCstock - Data_sd2$SoilCstock
PR_upper2 <- Data_mean3$SoilCstock + Data_sd1$SoilCstock
PR_lower2 <- Data_mean3$SoilCstock - Data_sd1$SoilCstock
p1 <- ggplot() + 
  geom_point(data = Data1yo, mapping = aes(x=Fertilized, y=SoilCstock, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide2, mapping = aes(x=1, xend=2, y=SoilCstock_Unfertilized, yend=SoilCstock_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  theme(text = element_text(size=9), legend.position = "none") + 
  geom_linerange(aes(x=c(0.9,2.1), ymin = PR_lower1, ymax=PR_upper1), col = c("#214478ff","#ff6600ff"), linewidth =1.25) +
  geom_point(aes(x=c(0.9,2.1), y=Data_mean4$SoilCstock), col = c("#214478ff","#ff6600ff"), size=3) +
  coord_cartesian(ylim = c(0,105)) + labs(x=NULL, y = bquote("Organic layer C stock (t "~ha^-1~")"), title = "1 year after clearcutting")
p2 <- ggplot() + 
  geom_point(data = Data4to13yo, mapping = aes(x=Fertilized, y=SoilCstock, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide, mapping = aes(x=1, xend=2, y=SoilCstock_Unfertilized, yend=SoilCstock_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  geom_linerange(aes(x=c(0.9,2.1), ymin = PR_lower2, ymax=PR_upper2), col = c("#214478ff","#ff6600ff"), linewidth =1.25) +
  geom_point(aes(x=c(0.9,2.1), y=Data_mean3$SoilCstock), col = c("#214478ff","#ff6600ff"), size=3) +
  theme(text = element_text(size=9), legend.position = "none") + #facet_wrap(~Pair) +
  coord_cartesian(ylim = c(0,105)) + labs(x=NULL, y = NULL, title = "4 - 13 years after clearcutting")
ggarrange(p1, p2, ncol=2)

PR_upper3 <- Data_mean4$SoilNstock + Data_sd2$SoilNstock
PR_lower3 <- Data_mean4$SoilNstock - Data_sd2$SoilNstock
PR_upper4 <- Data_mean3$SoilNstock + Data_sd1$SoilNstock
PR_lower4 <- Data_mean3$SoilNstock - Data_sd1$SoilNstock
p3 <- ggplot() + 
  geom_point(data = Data1yo, mapping = aes(x=Fertilized, y=SoilNstock, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide2, mapping = aes(x=1, xend=2, y=SoilNstock_Unfertilized, yend=SoilNstock_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  geom_linerange(aes(x=c(0.9,2.1), ymin = PR_lower3, ymax=PR_upper3), col = c("#214478ff","#ff6600ff"), linewidth =1.25) +
  geom_point(aes(x=c(0.9,2.1), y=Data_mean4$SoilNstock), col = c("#214478ff","#ff6600ff"), size=3) +
  theme(text = element_text(size=9), legend.position = "none") +
  coord_cartesian(ylim = c(0,3)) + labs(x=NULL, y= bquote("Organic layer N stock (t "~ha^-1~")"), title = "1 year after clearcutting")
p4 <- ggplot() + 
  geom_point(data = Data4to13yo, mapping = aes(x=Fertilized, y=SoilNstock, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide, mapping = aes(x=1, xend=2, y=SoilNstock_Unfertilized, yend=SoilNstock_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  geom_linerange(aes(x=c(0.9,2.1), ymin = PR_lower4, ymax=PR_upper4), col = c("#214478ff","#ff6600ff"), linewidth =1.25) +
  geom_point(aes(x=c(0.9,2.1), y=Data_mean3$SoilNstock), col = c("#214478ff","#ff6600ff"), size=3) +
  theme(text = element_text(size=9), legend.position = "none") +
  coord_cartesian(ylim = c(0,3)) + labs(x=NULL, y= NULL, title = "4 - 13 years after clearcutting")
ggarrange(p1, p2, p3, p4, ncol=4, widths = c(1.1,1,1.1,1))


p5 <- ggplot() + 
  geom_point(data = Data1yo, mapping = aes(x=Fertilized, y=pH, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide2, mapping = aes(x=1, xend=2, y=pH_Unfertilized, yend=pH_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  theme(text = element_text(size=9), legend.position = "none") +
  coord_cartesian(ylim = c(0,6)) + labs(x=NULL, y= "Soil pH", title = "1 year after clearcutting")
p6 <- ggplot() + 
  geom_point(data = Data4to13yo, mapping = aes(x=Fertilized, y=pH, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide, mapping = aes(x=1, xend=2, y=pH_Unfertilized, yend=pH_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  theme(text = element_text(size=9), legend.position = "none") +
  coord_cartesian(ylim = c(0,6)) + labs(x=NULL, y= "Soil pH", title = "4 - 13 years after clearcutting")
ggarrange(p5,p6,ncol=2)

p7 <- ggplot() + 
  geom_point(data = Data1yo, mapping = aes(x=Fertilized, y=OrgLayer, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide2, mapping = aes(x=1, xend=2, y=OrgLayer_Unfertilized, yend=OrgLayer_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  theme(text = element_text(size=9), legend.position = "none") +
  coord_cartesian(ylim = c(0,17)) + labs(x=NULL, y= "Soil organic layer depth", title = "1 year after clearcutting")
p8 <- ggplot() + 
  geom_point(data = Data4to13yo, mapping = aes(x=Fertilized, y=OrgLayer, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide, mapping = aes(x=1, xend=2, y=OrgLayer_Unfertilized, yend=OrgLayer_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  theme(text = element_text(size=9), legend.position = "none") +
  coord_cartesian(ylim = c(0,17)) + labs(x=NULL, y= "Soil organic layer depth", title = "4 - 13 years after clearcutting")
ggarrange(p7,p8,ncol=2)

p9 <- ggplot() + 
  geom_point(data = Data1yo, mapping = aes(x=Fertilized, y=CN, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide2, mapping = aes(x=1, xend=2, y=CN_Unfertilized, yend=CN_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  theme(text = element_text(size=9), legend.position = "none") +
  coord_cartesian(ylim = c(0,45)) + labs(x=NULL, y= "C:N Soil organic layer", title = "1 year after clearcutting")
p10 <- ggplot() + 
  geom_point(data = Data4to13yo, mapping = aes(x=Fertilized, y=CN, col=Fertilized), alpha=0.2, size=2) +
  geom_segment(data = Data_wide, mapping = aes(x=1, xend=2, y=CN_Unfertilized, yend=CN_Fertilized), alpha=0.5) +
  scale_color_manual(values=c("Fertilized" = "#ff6600ff", "Unfertilized" = "#214478ff")) +
  theme(text = element_text(size=9), legend.position = "none") +
  coord_cartesian(ylim = c(0,45)) + labs(x=NULL, y= "C:N Soil organic layer", title = "4 - 13 years after clearcutting")
ggarrange(p9,p10,ncol=2)