---
title: "Supplementary Material 3"
author: "Moritz Hallama"
output:
  html_document: 
   toc: true
   toc_float: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,warning = FALSE, message = FALSE)

```
# Online Resource S2: Fitted models and ANOVAs

This document contains details about the statistical analysis, fitted models, ANOVA tables and figures of the publication  
_"The role of microbes in the increase of organic phosphorus availability in the rhizosheath of cover crops"_  
(Plant and Soil 2022)
DOI: 10.1007/s11104-022-05340-5
For any questions, ask Moritz Hallama (https://orcid.org/0000-0003-4209-6760).

**No gods, no masters, no borders (A) :P**

```{r packages&dataset, include=F}
##### initializing R and dataset ####
library(readxl)
library(writexl)

library(lme4)
library(lmerTest)
library(plyr)
#library(reshape2)
library(emmeans)
library(multcomp)
library(tidyverse)
library(stringi)
library(viridis)
library(kableExtra)
library(performance)

library(highr)

#library(cowplot)

#citation("here") #how to cite packages
#install.packages("here") #to install packages not in your library

#detach(package:plyr)

library(ggpmisc)
# equation necessary to add regression coefficients to the scatter plots
lm_eqn <- function(df){
  m <- lm(y ~ x, df);
  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
                   list(a = format(unname(coef(m)[1]), digits = 2),
                        b = format(unname(coef(m)[2]), digits = 2),
                        r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq));
}
my.formula <- y ~ x

setwd("~/Nextcloud/PhD/!WP3_Article/Submission_PlantSoil/Submission_September") #this is one of the most delicate points. 
# You have to fill in the path of your saved files
# alternative: you call the package "here"
#library(here)

dataAll <- read_excel("Supplementary_Material_1.xlsx",sheet = "Data")
dataEAA <- read_excel("Supplementary_Material_1.xlsx",sheet = "EAA")
#View(dataAll)
#whole dataset
dataAll$Site<- as.factor(dataAll$Site)
dataAll$Type<- as.factor(dataAll$Type)
dataAll$Plot<- as.factor(dataAll$Plot)
dataAll$Date<- as.factor(dataAll$Date)
dataAll$Block<- as.factor(dataAll$Block)
dataAll$CoverCrop<- as.factor(dataAll$CoverCrop)
str(dataAll)

dataAll$Month<-revalue(dataAll$Date, c(T0="August",T1="November", T2="March", T3="June"))
dataAll$Month<-factor(dataAll$Month, c("August","November", "March", "June"))
dataAll$Date2<-revalue(dataAll$Date, c(T0="0",T1="10", T2="24", T3="40"))

dataAll$CoverCrop <- factor(dataAll$CoverCrop, levels=c("Control", "Buckwheat", "Mustard", "Phacelia"))

dataAll$TypeID<-as.factor(dataAll$Type:dataAll$Block)
dataAll$MonthID<-as.factor(dataAll$Month:dataAll$Block)

# Phospho- and neutral lipid fatty acid biomarkers (acording to Ruess & Chamberlain 2010)
dataAll$Gpos<-dataAll$`i15.0`+dataAll$`a15.0`+dataAll$`i16.0`+dataAll$`i17.0`
dataAll$Gneg<-dataAll$`cy17.0`+dataAll$`cy19.0`
dataAll$Bacteria<-dataAll$Gpos+dataAll$Gneg+dataAll$`P15.0`+dataAll$`P16.1n7`
#dataAll$Zygomycota <- dataAll$`P18.3n6` # not used
dataAll$Fungi<-dataAll$`P18.2n6`#+dataAll$Zygomycota
dataAll$AMF<-dataAll$`N16.1n5`
dataAll$Microbial<-dataAll$Bacteria+dataAll$Fungi+dataAll$`P16.1n5`

#ratios
dataAll$Fungi_bacteria<-dataAll$Fungi/dataAll$Bacteria
dataAll$Gneg_Gpos <- dataAll$Gneg/dataAll$Gpos
dataAll$CtoPmic <- dataAll$Microbial/dataAll$Pmic

dataAll_bulk=dataAll%>%subset(Type=="bulk"&Date!="T0")
dataAll_rhizoT1=dataAll%>%subset(Type2=="rhizo"&Date=="T1")
dataAll_rhizoT3=dataAll%>%subset(Type2=="rhizo"&Date=="T3")
dataAll_plant=dataAll%>%subset(CoverCrop!="Control")
dataAll_T1=dataAll%>%subset(Date=="T1")
dataAll_T3=dataAll%>%subset(Date=="T3")

# EAA (the enzyme addition dataset has analytical replicates)
dataEAA$Site<- as.factor(dataEAA$Site)
dataEAA$Date<- as.factor(dataEAA$Date)
dataEAA$Type<- as.factor(dataEAA$Type)
dataEAA$Plot<- as.factor(dataEAA$Plot)
dataEAA$Date<- as.factor(dataEAA$Date)
dataEAA$CoverCrop<- as.factor(dataEAA$CoverCrop)
dataEAA$Analysis<- as.factor(dataEAA$Analysis)

str(dataEAA)
dataEAA$Month<-revalue(dataEAA$Date, c(T1="November"))

dataEAA$CoverCrop <- factor(dataEAA$CoverCrop, levels=c("Control", "Buckwheat", "Mustard", "Phacelia"))

#It is necessary to resume the pseudoreplicates of the enzyme addition assay into a model to be able to use these values as together with the enzyme activity data for correlations

mod_pred_GP <-lmer(GP ~(1|Analysis)+Plot*Type, data=dataEAA, REML=T)

pred_GP<-as.data.frame(emmeans(mod_pred_GP, ~Plot|Type, adjust="bon", nesting=NULL))

pred_GP$Date<-as.factor("T1")

pred_GP$mrg<-pred_GP$Date:pred_GP$Type:pred_GP$Plot
dataAll_EAA<-dataAll
dataAll_EAA$mrg<-dataAll_EAA$Date:dataAll_EAA$Type:dataAll_EAA$Plot
#View(pred_GP)
pred_GP<- rename(pred_GP, "GP"="emmean")

dataAll_EAA<-merge(dataAll_EAA,pred_GP[,c(3,9)],by=c("mrg"))#, all.x=T)
str(dataAll_EAA)
mod_pred_DNA <-lmer(DNA ~(1|Analysis)+Plot*Type, data=dataEAA, REML=T)

pred_DNA<-as.data.frame(emmeans(mod_pred_DNA, ~Plot|Type, adjust="bon", nesting=NULL))

pred_DNA$Date<-as.factor("T1")

pred_DNA$mrg<-pred_DNA$Date:pred_DNA$Type:pred_DNA$Plot

pred_DNA<- rename(pred_DNA, "DNA"="emmean")
#View(pred_DNA)
dataAll_EAA<-merge(dataAll_EAA,pred_DNA[,c(3,9)],by=c("mrg"))#, all.x=T)

mod_pred_PHYF <-lmer(PHYF ~(1|Analysis)+Plot*Type, data=dataEAA, REML=T)

pred_PHYF<-as.data.frame(emmeans(mod_pred_PHYF, ~Plot|Type, adjust="bon", nesting=NULL))

pred_PHYF$Date<-as.factor("T1")

pred_PHYF$mrg<-pred_PHYF$Date:pred_PHYF$Type:pred_PHYF$Plot

pred_PHYF<- rename(pred_PHYF, "PHYF"="emmean")
dataAll_EAA<-merge(dataAll_EAA,pred_PHYF[,c(3,9)],by=c("mrg"))#, all.x=T)

mod_pred_PHYB <-lmer(PHYB ~(1|Analysis)+Plot*Type, data=dataEAA, REML=T)

pred_PHYB<-as.data.frame(emmeans(mod_pred_PHYB, ~Plot|Type, adjust="bon", nesting=NULL))
pred_PHYB$Date<-as.factor("T1")

pred_PHYB$mrg<-pred_PHYB$Date:pred_PHYB$Type:pred_PHYB$Plot
#View(pred_PHYB)
pred_PHYB<- rename(pred_PHYB, "PHYB"="emmean")

dataAll_EAA<-merge(dataAll_EAA,pred_PHYB[,c(3,9)],by=c("mrg"))#, all.x=T)

mod_pred_total_enzyme_labile_P <-lmer(total_enzyme_labile_P ~(1|Analysis)+Plot*Type, data=dataEAA, REML=T)

pred_total_enzyme_labile_P<-as.data.frame(emmeans(mod_pred_total_enzyme_labile_P, ~Plot|Type, adjust="bon", nesting=NULL))
pred_total_enzyme_labile_P$Date<-as.factor("T1")

pred_total_enzyme_labile_P$mrg<-pred_total_enzyme_labile_P$Date:pred_total_enzyme_labile_P$Type:pred_total_enzyme_labile_P$Plot

pred_total_enzyme_labile_P<- rename(pred_total_enzyme_labile_P, "total_enzyme_labile_P"="emmean")

dataAll_EAA<-merge(dataAll_EAA,pred_total_enzyme_labile_P[,c(3,9)],by=c("mrg"))#, all.x=T)

#str(dataAll_EAA)
#subset(dataAll_EAA, Date=="T1")
#str(dataAll)
#View(dataAll_EAA)

#visual comparison of raw data (with analytical pseudoreplicates) and modelled averages
boxPHYF_pred <-ggplot(dataAll_EAA , aes(x=CoverCrop, y=PHYF, fill=CoverCrop))+
  geom_boxplot()+
  geom_text(data = dataAll_EAA, aes(x = CoverCrop, group=CoverCrop, y = PHYF, label = Plot), position=position_dodge(1), hjust=.5)+
  #Motheme_all+
  xlab(" Phytase-labile Porg (averaged) ")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
boxPHYF_pred

boxPHYF <-ggplot(dataEAA , aes(x=CoverCrop, y=PHYF, fill=CoverCrop))+
  geom_boxplot()+
  geom_text(data = dataEAA, aes(x = CoverCrop, group=CoverCrop, y = PHYF, label = Plot:Analysis), position=position_dodge(1), hjust=.5)+
  # Motheme_all+
  xlab(" Phytase-labile Porg (Analytical replicates) ")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
boxPHYF

#acceptable

###### Plot appearance ########
library(viridis)
n_distinct(dataAll$CoverCrop) #richness for figures (to choose colors)
viridis(n_distinct(dataAll$CoverCrop))# viridis selected for aesthetic reasons (reordered and purple lighted up slightly)
#scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
viridis(n_distinct(dataAll$CoverCrop:dataAll$Type))# viridis selected for aesthetic reasons (reordered and purple lighted up slightly)

# colors for shoot vs root
library(RColorBrewer)
display.brewer.all(n=10, exact.n=FALSE)
brewer.pal(11,"BrBG")
display.brewer.pal(11,"BrBG")
#scale_fill_manual("Soil compartment",values=c("#35978F","#BF812D"))+

# plot appearance for most figures
Motheme_all <- theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1,size=12, color="black"), 
                     axis.text.y= element_text(color="black", face="bold", size=rel(1.3)), 
                     plot.title= element_text(face="bold",size=rel(1.2)), 
                     axis.title.y= element_text(face="bold",size=rel(1.1)), 
                     strip.text = element_text(face="bold", size=rel(1.1)),
                     legend.text=element_text(size=12),
                     panel.background=element_blank(),
                     axis.line=element_line(color="black"),
                     panel.grid.major.y=element_line(color="light grey", size=0.2),
                     panel.grid.major.x = element_blank())

# plot appearance for scatter plots
Motheme_all2 <- theme(axis.text.x = element_text(size=20, color="black"), 
                      axis.text.y= element_text(color="black", size=20), 
                      axis.title= element_text(size=21), 
                      plot.title= element_text(size=22),
                      strip.text = element_text( size=rel(1.1)),
                      legend.text=element_text(size=20),
                      legend.title=element_text(size=20),
                      panel.background=element_blank(),
                      axis.line=element_line(color="black"),
                      panel.grid.major.y=element_line(color="light grey", size=0.2),
                      panel.grid.major.x = element_blank())

#plot appearance for Fig3
Motheme_EAA <- theme(plot.title= element_text(face="bold",size=14), 
                     axis.text.x = element_text(angle=30, hjust=1, vjust=1,size=18, color="black"), 
                     axis.text.y= element_text(color="black", size=18),
                     axis.title.y= element_text(face="bold",size=18), 
                     strip.text = element_text(face="bold", size=18, 
                                               margin = margin(0.15,0,0.15,0, "cm")),
                     legend.text=element_text(size=18),
                     panel.grid.major.y = element_line(color="grey"))


```

# Fig 2 Cover crop biomass properties ####

###  Fig 2a1 Cover crop root biomass

``` {r Fig 2a1,fig.width = 4, fig.height=4, echo=F} 
######## 2a.1 CC_root_BM_kg_ha ########
#Boxplot for data inspection
ggplot(dataAll_plant , aes(x=CoverCrop, y=CC_root_BM_kg_ha, fill=CoverCrop))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(kg h",a^-1,")")))+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
# model
CC_root_BM_kg_ha_model.full.lm<- lm((CC_root_BM_kg_ha)~CoverCrop+Block, na.action=na.omit,
                                    data=dataAll_plant)
anova(CC_root_BM_kg_ha_model.full.lm )
# CoverCrop is  sig!
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modCC_root_BM_kg_ha<-CC_root_BM_kg_ha_model.full.lm
plot(fitted(modCC_root_BM_kg_ha),resid(modCC_root_BM_kg_ha),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modCC_root_BM_kg_ha)) 
qqline(y=resid(modCC_root_BM_kg_ha))
# histogram of residuals
hist(resid(modCC_root_BM_kg_ha))

shapiro.test(resid(modCC_root_BM_kg_ha))

# estimating marginal means
modCC_root_BM_kg_ha_graph <-as.data.frame(multcomp::cld(emmeans(modCC_root_BM_kg_ha, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```

```{r, fig.width = 6, echo=F}
CC_root_BM_kg_ha_graph <- ggplot(subset(modCC_root_BM_kg_ha_graph), aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  ggtitle("") +
  ylab(expression(paste("Cover crop root biomass (kg h",a^-1, ")")))+
  theme(legend.position="none")+xlab("")+
  ylim( min(modCC_root_BM_kg_ha_graph$lower.CL),max(modCC_root_BM_kg_ha_graph$upper.CL)+max(modCC_root_BM_kg_ha_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
CC_root_BM_kg_ha_graph
```

##### Model fit

```{r, echo=F}
formula(modCC_root_BM_kg_ha)

anova(modCC_root_BM_kg_ha) %>%
  kable(caption= "Cover crop root biomass: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,  include=F}
svg("CC_root_BM_kg_ha_graph_leg.svg",width=8,height=8)
CC_root_BM_kg_ha_graph
dev.off()

png("CC_root_BM_kg_ha_graph.png",width=12,height=12, units = "cm", res=800)
CC_root_BM_kg_ha_graph
dev.off()
```
\newpage

###  Fig 2a2 Cover crop shoot biomass

``` {r Fig 2a2,fig.width = 4, fig.height=4,echo=F }  
######## 2a.2 CC_shoot_BM_kg_ha ########
ggplot(dataAll_plant , aes(x=CoverCrop, y=CC_shoot_BM_kg_ha, fill=CoverCrop))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(kg h",a^-1,")")))+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
CC_shoot_BM_kg_ha_model.full.lm<- lm(log(CC_shoot_BM_kg_ha)~CoverCrop+Block, na.action=na.omit,
                                     data=dataAll_plant)
anova(CC_shoot_BM_kg_ha_model.full.lm )
# CoverCrop is  sig!
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modCC_shoot_BM_kg_ha<-CC_shoot_BM_kg_ha_model.full.lm
plot(fitted(modCC_shoot_BM_kg_ha),resid(modCC_shoot_BM_kg_ha),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modCC_shoot_BM_kg_ha)) 
qqline(y=resid(modCC_shoot_BM_kg_ha))
# histogram of residuals
hist(resid(modCC_shoot_BM_kg_ha))

shapiro.test(resid(modCC_shoot_BM_kg_ha))

modCC_shoot_BM_kg_ha_graph <-as.data.frame(multcomp::cld(emmeans(modCC_shoot_BM_kg_ha, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```

```{r, fig.width = 6, echo=F }

CC_shoot_BM_kg_ha_graph <- ggplot(subset(modCC_shoot_BM_kg_ha_graph), aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste("Cover crop shoot biomass (kg h",a^-1, ")")))+
  ylim( min(modCC_shoot_BM_kg_ha_graph$lower.CL),max(modCC_shoot_BM_kg_ha_graph$upper.CL)+max(modCC_shoot_BM_kg_ha_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
CC_shoot_BM_kg_ha_graph
```

##### Model fit

```{r, echo=F }
formula(modCC_shoot_BM_kg_ha)

anova(modCC_shoot_BM_kg_ha) %>%
  kable(caption= "Cover crop shoot biomass: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,  include=F}
svg("CC_shoot_BM_kg_ha_graph.svg",width=8,height=8)
CC_shoot_BM_kg_ha_graph
dev.off()

png("CC_shoot_BM_kg_ha_graph.png",width=12,height=12, units = "cm", res=800)
CC_shoot_BM_kg_ha_graph+  theme(legend.position="none")+xlab("")
dev.off()
```
\newpage
###  Fig 2a3 Total cover crop biomass

``` {r Fig 2a3,fig.width = 4, fig.height=4,echo=F}  
######## 2a.3 CC_total_dry_BM_kg_ha ########
ggplot(dataAll_plant , aes(x=CoverCrop, y=CC_total_dry_BM_kg_ha, fill=CoverCrop))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(kg h",a^-1,")")))+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
CC_total_dry_BM_kg_ha_model.full.lm<- lm(log(CC_total_dry_BM_kg_ha)~CoverCrop+Block, na.action=na.omit,
                                         data=dataAll_plant)
anova(CC_total_dry_BM_kg_ha_model.full.lm )
# CoverCrop is  sig!

modCC_total_dry_BM_kg_ha<-CC_total_dry_BM_kg_ha_model.full.lm
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modCC_total_dry_BM_kg_ha),resid(modCC_total_dry_BM_kg_ha),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modCC_total_dry_BM_kg_ha)) 
qqline(y=resid(modCC_total_dry_BM_kg_ha))
# histogram of residuals
hist(resid(modCC_total_dry_BM_kg_ha))

shapiro.test(resid(modCC_total_dry_BM_kg_ha))

modCC_total_dry_BM_kg_ha_graph <-as.data.frame(multcomp::cld(emmeans(modCC_total_dry_BM_kg_ha, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```

```{r, fig.width = 6, echo=F }

CC_total_dry_BM_kg_ha_graph <- ggplot(subset(modCC_total_dry_BM_kg_ha_graph), aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste( "Total cover crop biomass (kg h",a^-1,")")))+
  theme(legend.position="none")+xlab("")+
  ylim( min(modCC_total_dry_BM_kg_ha_graph$lower.CL),max(modCC_total_dry_BM_kg_ha_graph$upper.CL)+max(modCC_total_dry_BM_kg_ha_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
CC_total_dry_BM_kg_ha_graph
```

##### Model fit

```{r, echo=F}
formula(modCC_total_dry_BM_kg_ha)

anova(modCC_total_dry_BM_kg_ha) %>%
  kable(caption= "Total cover crop biomass: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,  include=F}
svg("CC_total_dry_BM_kg_ha_graph.svg",width=8,height=8)
CC_total_dry_BM_kg_ha_graph
dev.off()

png("CC_total_dry_BM_kg_ha_graph.png",width=12,height=12, units = "cm", res=800)
CC_total_dry_BM_kg_ha_graph
dev.off()
```
\newpage
## Fig 2a stacked bars: Cover crop biomass #####

``` {r Fig 2a, echo=F}
modCC_root_BM_kg_ha_graph <-as.data.frame(multcomp::cld(emmeans(modCC_root_BM_kg_ha, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
modCC_shoot_BM_kg_ha_graph <-as.data.frame(multcomp::cld(emmeans(modCC_shoot_BM_kg_ha, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

modCC_root_BM_kg_ha_graph$PlantPart<-"root"
modCC_shoot_BM_kg_ha_graph$PlantPart<-"shoot"

modCC_shoot_BM_kg_ha_graph$lower.CL <-modCC_shoot_BM_kg_ha_graph$lower.CL+modCC_root_BM_kg_ha_graph$response
modCC_shoot_BM_kg_ha_graph$upper.CL <-modCC_shoot_BM_kg_ha_graph$upper.CL+modCC_root_BM_kg_ha_graph$response

stacked_CC_BM<-rbind(modCC_root_BM_kg_ha_graph, modCC_shoot_BM_kg_ha_graph)
stacked_CC_BM$PlantPart<-as.factor(stacked_CC_BM$PlantPart)
stacked_CC_BM$PlantPart<-factor(stacked_CC_BM$PlantPart, levels=c("shoot", "root"))

stacked_CC_BM_graph <- ggplot(stacked_CC_BM, aes(x=CoverCrop, y=response, fill=PlantPart, order=desc(PlantPart)))+
  geom_bar(stat="identity")+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste( "Cover crop biomass (kg DM h",a^-1,")")))+
  xlab("")+
  scale_fill_manual("Plant part",values=c("#35978F","#BF812D"))+
  #ylim( min(modCC_root_graph$lower.CL),max(modCC_root_graph$upper.CL)+max(modCC_root_graph$upper.CL)/10)+
  Motheme_all
stacked_CC_BM_graph
```

```{r, include=F}
svg("stacked_CC_BM.svg",width=8,height=8)
stacked_CC_BM_graph
dev.off()

png("stacked_CC_BM.png",width=16,height=16, units = "cm", res=800)
stacked_CC_BM_graph#+theme(legend.position="none")
dev.off()
```

\newpage
####  Fig 2b1 P concentration in cover crop shoot biomass

``` {r Fig 2b1, fig.width = 4, fig.height=4,echo=F}  
######## 2b.1 P_CC_shoot concentration ########
ggplot(dataAll_plant , aes(x=CoverCrop, y=P_CC_shoot, fill=CoverCrop))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(mg P k",g^-1,")")))+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
P_CC_shoot_model.full.lm<- lm((P_CC_shoot)~CoverCrop+Block, na.action=na.omit,
                              data=dataAll_plant)
anova(P_CC_shoot_model.full.lm )
# CoverCrop is  sig!
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modP_CC_shoot<-P_CC_shoot_model.full.lm
plot(fitted(modP_CC_shoot),resid(modP_CC_shoot),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modP_CC_shoot)) 
qqline(y=resid(modP_CC_shoot))
# histogram of residuals
hist(resid(modP_CC_shoot))

shapiro.test(resid(modP_CC_shoot))

modP_CC_shoot_graph <-as.data.frame(multcomp::cld(emmeans(modP_CC_shoot, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```

```{r, fig.width = 6, echo=F }

P_CC_shoot_graph <- ggplot(subset(modP_CC_shoot_graph), aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual("Cover Crop",values=c(22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste( "Shoot P concentration (mg P k",g^-1,")")))+
  theme(legend.position="none")+
  xlab("")+
  ylim( min(modP_CC_shoot_graph$lower.CL),max(modP_CC_shoot_graph$upper.CL)+max(modP_CC_shoot_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
P_CC_shoot_graph

```

##### Model fit

```{r, echo=F}
formula(modP_CC_shoot)

anova(modP_CC_shoot) %>%
  kable(caption= "P concentration in cover crop shoots: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,   include=F}
svg("P_CC_shoot_graph.svg",width=8,height=8)
P_CC_shoot_graph
dev.off()

png("P_CC_shoot_graph.png",width=12,height=12, units = "cm", res=800)
P_CC_shoot_graph
dev.off()
```
\newpage
###  Fig 2b2 P concentration in cover crop root biomass

``` {r Fig 2b2,fig.width = 4, fig.height=4,echo=F}  
######## 2b.2 P_CC_root concentration ########
ggplot(dataAll_plant , aes(x=CoverCrop, y=P_CC_root, fill=CoverCrop))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(mg P k",g^-1,")")))+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
P_CC_root_model.full.lm<- lm((P_CC_root)~CoverCrop+Block, na.action=na.omit,
                             data=dataAll_plant)
anova(P_CC_root_model.full.lm )
# CoverCrop is  sig!
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modP_CC_root<-P_CC_root_model.full.lm
plot(fitted(modP_CC_root),resid(modP_CC_root),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modP_CC_root)) 
qqline(y=resid(modP_CC_root))
# histogram of residuals
hist(resid(modP_CC_root))

shapiro.test(resid(modP_CC_root))

modP_CC_root_graph <-as.data.frame(multcomp::cld(emmeans(modP_CC_root, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```

```{r, fig.width = 6, echo=F }

P_CC_root_graph <- ggplot(subset(modP_CC_root_graph), aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual("Cover Crop",values=c(22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste( "Root P concentration (mg P k",g^-1,")")))+
  #theme(legend.position="none")+
  xlab("")+
  ylim( min(modP_CC_root_graph$lower.CL),max(modP_CC_root_graph$upper.CL)+max(modP_CC_root_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
P_CC_root_graph
```

##### Model fit

```{r, echo=F}
formula(modP_CC_root)

anova(modP_CC_root) %>%
  kable(caption= "P concentration in cover crop roots: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,   include=F}
svg("P_CC_root_graph.svg",width=8,height=8)
P_CC_root_graph
dev.off()

png("P_CC_root_graph.png",width=12,height=12, units = "cm", res=800)
P_CC_root_graph
dev.off()
```
\newpage
## Fig 2b stacked bars P concentration #####

``` {r Fig 2b, echo=F}
modP_CC_shoot_graph <-as.data.frame(multcomp::cld(emmeans(modP_CC_shoot, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

modP_CC_root_graph <-as.data.frame(multcomp::cld(emmeans(modP_CC_root, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

modP_CC_root_graph$PlantPart<-"root"
modP_CC_shoot_graph$PlantPart<-"shoot"
P_in_CC_biomass<-rbind(modP_CC_root_graph, modP_CC_shoot_graph)
P_in_CC_biomass$PlantPart<-factor(P_in_CC_biomass$PlantPart, levels=c( "shoot","root"))

P_in_CC_biomass_graph <- ggplot(P_in_CC_biomass, aes(x=CoverCrop, y=response/1000, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual("Cover Crop",values=c(22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL/1000, ymax=upper.CL/1000), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste( "Cover crop P concentration (g P k",g^-1,")")))+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( min(modP_CC_root_graph$lower.CL),max(modP_CC_root_graph$upper.CL)+max(modP_CC_root_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  facet_grid(PlantPart~.)+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
P_in_CC_biomass_graph
```

``` {r,  include=F}
svg("P_in_CC_biomass.svg",width=8,height=8)
P_in_CC_biomass_graph
dev.off()

png("P_in_CC_biomass.png",width=16,height=16, units = "cm", res=800)
P_in_CC_biomass_graph
dev.off()
```

\newpage
## Fig2c1 P content cover crop shoot biomass

``` {r Fig 2c1,fig.width = 4, fig.height=4,echo=F}  
######## 2c.1 P_in_CC_shoot ########
ggplot(dataAll_plant , aes(x=CoverCrop, y=P_in_CC_shoot, fill=CoverCrop))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(kg  h",a^-1,")")))+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
P_in_CC_shoot_model.full.lm<- lm(log(P_in_CC_shoot)~CoverCrop+Block, na.action=na.omit,
                                 data=dataAll_plant)
anova(P_in_CC_shoot_model.full.lm )
# CoverCrop is  sig!
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modP_in_CC_shoot<-P_in_CC_shoot_model.full.lm
plot(fitted(modP_in_CC_shoot),resid(modP_in_CC_shoot),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modP_in_CC_shoot)) 
qqline(y=resid(modP_in_CC_shoot))
# histogram of residuals
hist(resid(modP_in_CC_shoot))

shapiro.test(resid(modP_in_CC_shoot))

modP_in_CC_shoot_graph <-as.data.frame(multcomp::cld(emmeans(modP_in_CC_shoot, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
```

```{r, fig.width = 6, echo=F }

P_in_CC_shoot_graph <- ggplot(subset(modP_in_CC_shoot_graph), aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste( "Shoot P (kg h",a^-1,")")))+
  theme(legend.position="none")+xlab("")+
  ylim( min(modP_in_CC_shoot_graph$lower.CL),max(modP_in_CC_shoot_graph$upper.CL)+max(modP_in_CC_shoot_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
P_in_CC_shoot_graph
```

##### Model fit

```{r, echo=F}

formula(modP_in_CC_shoot)

anova(modP_in_CC_shoot) %>%
  kable(caption= "P content in cover crop shoots: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,  include=F}
svg("P_in_CC_shoot_graph.svg",width=8,height=8)
P_in_CC_shoot_graph
dev.off()

png("P_in_CC_shoot_graph.png",width=12,height=12, units = "cm", res=800)
P_in_CC_shoot_graph
dev.off()
```
\newpage
## Fig 2c2 P content in cover crop root biomass

``` {r Fig 2c2,fig.width = 4, fig.height=4,echo=F}  
######## 2c.2 P_in_CC_root ########
ggplot(dataAll_plant , aes(x=CoverCrop, y=P_in_CC_root, fill=CoverCrop))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(kg P h",a^-1,")")))+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
P_in_CC_root_model.full.lm<- lm(log(P_in_CC_root)~CoverCrop+Block, na.action=na.omit,
                                data=dataAll_plant)
anova(P_in_CC_root_model.full.lm )
# CoverCrop is  sig!
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modP_in_CC_root<-P_in_CC_root_model.full.lm
plot(fitted(modP_in_CC_root),resid(modP_in_CC_root),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modP_in_CC_root)) 
qqline(y=resid(modP_in_CC_root))
# histogram of residuals
hist(resid(modP_in_CC_root))

shapiro.test(resid(modP_in_CC_root))

modP_in_CC_root_graph <-as.data.frame(multcomp::cld(emmeans(modP_in_CC_root, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
```

```{r, fig.width = 6, echo=F }

P_in_CC_root_graph <- ggplot(subset(modP_in_CC_root_graph), aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste( "Root P (kg h",a^-1,")")))+
  theme(legend.position="none")+xlab("")+
  ylim( min(modP_in_CC_root_graph$lower.CL),max(modP_in_CC_root_graph$upper.CL)+max(modP_in_CC_root_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#693376FF", "#FDE725FF", "#35B779FF"))
P_in_CC_root_graph
```

##### Model fit

```{r, echo=F}
formula(modP_in_CC_root)

anova(modP_in_CC_root) %>%
  kable(caption= "P content in cover crop roots: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,  include=F}
svg("P_in_CC_root_graph.svg",width=8,height=8)
P_in_CC_root_graph
dev.off()

png("P_in_CC_root_graph.png",width=12,height=12, units = "cm", res=800)
P_in_CC_root_graph
dev.off()
```
\newpage
## Fig 2c stacked bars: CC P content #####

``` {r Fig 2c,echo=F}
modP_in_CC_root_graph <-as.data.frame(multcomp::cld(emmeans(modP_in_CC_root, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
modP_in_CC_shoot_graph <-as.data.frame(multcomp::cld(emmeans(modP_in_CC_shoot, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

modP_in_CC_root_graph$PlantPart<-"root"
modP_in_CC_shoot_graph$PlantPart<-"shoot"

modP_in_CC_shoot_graph$lower.CL <-modP_in_CC_shoot_graph$lower.CL+modP_in_CC_root_graph$response
modP_in_CC_shoot_graph$upper.CL <-modP_in_CC_shoot_graph$upper.CL+modP_in_CC_root_graph$response

stacked_P_in_CC_BM<-rbind(modP_in_CC_root_graph, modP_in_CC_shoot_graph)
stacked_P_in_CC_BM$PlantPart<-as.factor(stacked_P_in_CC_BM$PlantPart)
stacked_P_in_CC_BM$PlantPart<-factor(stacked_P_in_CC_BM$PlantPart, levels=c("shoot", "root"))

stacked_P_in_CC_BM_graph <- ggplot(stacked_P_in_CC_BM, aes(x=CoverCrop, y=response, fill=PlantPart, order=desc(PlantPart)))+
  geom_bar(stat="identity")+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste( "Cover crop P content (kg h",a^-1,")")))+
  #theme(legend.position="none")+
  xlab("")+
  scale_fill_manual("Plant part",values=c("#35978F","#BF812D"))+
  #ylim( min(modP_in_CC_root_graph$lower.CL),max(modP_in_CC_root_graph$upper.CL)+max(modP_in_CC_root_graph$upper.CL)/10)+
  Motheme_all
stacked_P_in_CC_BM_graph

```

```{r, include=F}

svg("stacked_P_in_CC_BM.svg",width=8,height=8)
stacked_P_in_CC_BM_graph
dev.off()

png("stacked_P_in_CC_BM.png",width=16,height=16, units = "cm", res=800)
stacked_P_in_CC_BM_graph
dev.off()
```
\newpage
# Fig 3. Soil P pools in the rhizosphere of cover crops

### Total P

``` {r Total P,fig.width = 4, fig.height=4,echo=F}  
######## 3.1 EAA_Pt ########
ggplot(dataEAA , aes(x=CoverCrop, y=Pt, fill=CoverCrop))+
  geom_boxplot()+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```
```{r, include=F}
range(dataEAA$Pt)
EAA_Pt_model.full<- lm(Pt~CoverCrop+Block, na.action=na.omit,
                       data=dataAll)
anova(EAA_Pt_model.full)
EAA_Pt_model_type.full<- lm(Pt~Type+Block, na.action=na.omit,
                            data=dataAll)
anova(EAA_Pt_model_type.full)

anova(EAA_Pt_model.full,EAA_Pt_model_type.full )

modPt<-EAA_Pt_model.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modPt),resid(modPt),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modPt)) 
qqline(y=resid(modPt))
# histogram of residuals
hist(resid(modPt))

shapiro.test(resid(modPt))

modPt_graph <-as.data.frame(multcomp::cld(emmeans(modPt, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```
##### Model fit 

```{r, echo=F}
formula(modPt)

anova(modPt) %>%
  kable(caption= "Total P: ANOVA output of fitted model") %>%
  kable_styling()
```
\newpage
### Organic P
``` {r Organic P,fig.width = 4, fig.height=4,echo=F}  
######## 3.2 EAA_Porg ########
ggplot(dataEAA , aes(x=CoverCrop, y=Porg, fill=CoverCrop))+
  geom_text(data = dataEAA, aes(x = CoverCrop, group=CoverCrop, y = Porg, label = Block), position=position_dodge(1), hjust=.5)+
  geom_boxplot()+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
EAA_Porg_model.full<- lmer((Porg)~CoverCrop+Block+(1|Plot), na.action=na.omit,
                           data=dataEAA)
anova(EAA_Porg_model.full)
VarCorr(EAA_Porg_model.full)
Porg_redstep<-get_model(step(EAA_Porg_model.full, direction="both",  reduce.random=F,keep=c("Block")))

anova(EAA_Porg_model.full,Porg_redstep)

modPorg<-EAA_Porg_model.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modPorg),resid(modPorg),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modPorg)) 
qqline(y=resid(modPorg))
# histogram of residuals
hist(resid(modPorg))

shapiro.test(resid(modPorg))
```
##### Model fit 

```{r, echo=F}
formula(modPorg)

anova(modPorg) %>%
  kable(caption= "Organic P: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
### inorganic P

``` {r Inorganic P,fig.width = 4, fig.height=4, echo=F}  
######## 3 EAA_Pi ########
ggplot(dataEAA , aes(x=CoverCrop, y=Pi, fill=CoverCrop))+
  geom_text(data = dataEAA, aes(x = CoverCrop, group=CoverCrop, y = Pi, label = Plot:Analysis), position=position_dodge(1), hjust=.5)+
  geom_boxplot()+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
EAA_Pi_model.full<- lmer(Pi~CoverCrop+Block+(1|Plot), na.action=na.omit,
                         data=dataEAA)
anova(EAA_Pi_model.full)
VarCorr(EAA_Pi_model.full)

# CoverCrop not sig

modPi<-EAA_Pi_model.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modPi),resid(modPi),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modPi)) 
qqline(y=resid(modPi))
# histogram of residuals
hist(resid(modPi))

shapiro.test(resid(modPi))

modPi_graph <-as.data.frame(multcomp::cld(emmeans(modPi, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```
\newpage
## Fig 3a1 Inorganic P
```{r Fig 3a1, fig.width = 6, echo=F }

ggplot(subset(modPi_graph), aes(x=CoverCrop, y=emmean, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(21,22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)+10), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("Inorganic P") +ylab(expression(paste("(",mu,"g ", P[i]," ",g^-1,")",sep="")))+ labs(fill="")+
  theme(legend.position="none")+xlab("")+
  ylim( min(modPi_graph$lower.CL),max(modPi_graph$upper.CL)+max(modPi_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

##### Model fit 

```{r, echo=F}
formula(modPi)

anova(modPi) %>%
  kable(caption= "Inorganic P: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig 3a2 Enzyme-labile P

``` {r Fig 3a2,fig.width = 4, fig.height=4,echo=F}  
######## 3.4 EAA_total_enzyme_labile_P ########
ggplot(dataEAA , aes(x=CoverCrop, y=total_enzyme_labile_P, fill=CoverCrop))+
  geom_boxplot()+
  geom_text(data = dataEAA, aes(x = CoverCrop, group=CoverCrop, y = total_enzyme_labile_P, label = Plot:Analysis), position=position_dodge(1), hjust=.5)+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
range(na.exclude(dataEAA$total_enzyme_labile_P))
summary(dataEAA$total_enzyme_labile_P)

EAA_total_enzyme_labile_P_model.full<- lmer((total_enzyme_labile_P)~CoverCrop+Block+(1|Plot)+(1|Analysis), na.action=na.omit, data=dataEAA)
anova(EAA_total_enzyme_labile_P_model.full)
VarCorr(EAA_total_enzyme_labile_P_model.full)

EAA_total_enzyme_labile_P_model.full2<- lmer(total_enzyme_labile_P~CoverCrop+Block+(1|Plot), na.action=na.omit, data=dataEAA)

anova(EAA_total_enzyme_labile_P_model.full2)
VarCorr(EAA_total_enzyme_labile_P_model.full2)

EAA_total_enzyme_labile_P_model.full3<- lmer(total_enzyme_labile_P~CoverCrop+Block+(1|Analysis), na.action=na.omit, data=dataEAA)

anova(EAA_total_enzyme_labile_P_model.full3)
VarCorr(EAA_total_enzyme_labile_P_model.full3)

multcomp::cld(emmeans(lmer(total_enzyme_labile_P~Type+Block+(1|Analysis), na.action=na.omit, data=dataEAA),~Type))
269/216 #average change between rhizosphere and bulk soil

# using soil compartment ("Type") instead of cover crops
EAA_total_enzyme_labile_P_model_type<- lmer(total_enzyme_labile_P~Type+Block+(1|Analysis), na.action=na.omit,
                                            data=dataEAA)
anova(EAA_total_enzyme_labile_P_model_type)
cld(emmeans(EAA_total_enzyme_labile_P_model_type, ~Type, adjust="bon"),Letters=letters, type = "response")

anova(EAA_total_enzyme_labile_P_model.full, EAA_total_enzyme_labile_P_model.full2, EAA_total_enzyme_labile_P_model.full3 )

# full3 has lowest AIC: CoverCrops is sig! 

modtotal_enzyme_labile_P<-EAA_total_enzyme_labile_P_model.full3
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modtotal_enzyme_labile_P),resid(modtotal_enzyme_labile_P),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modtotal_enzyme_labile_P)) 
qqline(y=resid(modtotal_enzyme_labile_P))
# histogram of residuals
hist(resid(modtotal_enzyme_labile_P))

shapiro.test(resid(modtotal_enzyme_labile_P))

modtotal_enzyme_labile_P_graph <-as.data.frame(multcomp::cld(emmeans(modtotal_enzyme_labile_P, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```

```{r, fig.width = 6, echo=F }
ggplot(subset(modtotal_enzyme_labile_P_graph), aes(x=CoverCrop, y=emmean, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(21,22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)+10), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("Total enzyme-labile P") +ylab(expression(paste("(",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  theme(legend.position="none")+xlab("")+
  ylim( min(modtotal_enzyme_labile_P_graph$lower.CL),max(modtotal_enzyme_labile_P_graph$upper.CL)+max(modtotal_enzyme_labile_P_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))

```

##### Model fit

```{r, echo=F}
formula(modtotal_enzyme_labile_P)

anova(modtotal_enzyme_labile_P) %>%
  kable(caption= "Total enzyme-labile organic P: ANOVA output of fitted model") %>%
  kable_styling()
```
\newpage
## Fig 3a3 Enzyme stable P

``` {r Fig 3a3,fig.width = 4, fig.height=4,echo=F}  
######## 3.5 EAA_enzyme_stable_MUP ########
ggplot(dataEAA , aes(x=CoverCrop, y=enzyme_stable_MUP, fill=CoverCrop))+
  geom_boxplot()+
  geom_text(data = dataEAA, aes(x = CoverCrop, group=CoverCrop, y = enzyme_stable_MUP, label = Plot:Analysis), position=position_dodge(1), hjust=.5)+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
EAA_enzyme_stable_MUP_model.full<- lmer(enzyme_stable_MUP~CoverCrop+Block+(1|Plot)+(1|Analysis), na.action=na.omit, data=dataEAA)
anova(EAA_enzyme_stable_MUP_model.full)
VarCorr(EAA_enzyme_stable_MUP_model.full)

EAA_enzyme_stable_MUP_model.full<- lmer(enzyme_stable_MUP~CoverCrop+Block+(1|Plot), na.action=na.omit, data=dataEAA)

anova(EAA_enzyme_stable_MUP_model.full)
VarCorr(EAA_enzyme_stable_MUP_model.full)

EAA_enzyme_stable_MUP_model.full3<- lmer(enzyme_stable_MUP~CoverCrop+Block+(1|Analysis), na.action=na.omit, data=dataEAA)

anova(EAA_enzyme_stable_MUP_model.full3)
VarCorr(EAA_enzyme_stable_MUP_model.full3)

anova(EAA_enzyme_stable_MUP_model.full, EAA_enzyme_stable_MUP_model.full, EAA_enzyme_stable_MUP_model.full3 )

# Cover Crops not sig!

#modenzyme_stable_MUP<-enzyme_stable_MUP_redstep
modenzyme_stable_MUP<-EAA_enzyme_stable_MUP_model.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modenzyme_stable_MUP),resid(modenzyme_stable_MUP),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modenzyme_stable_MUP)) 
qqline(y=resid(modenzyme_stable_MUP))
# histogram of residuals
hist(resid(modenzyme_stable_MUP))

shapiro.test(resid(modenzyme_stable_MUP))

modenzyme_stable_MUP_graph <-as.data.frame(multcomp::cld(emmeans(modenzyme_stable_MUP, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
```

```{r, fig.width = 6, echo=F }

ggplot(subset(modenzyme_stable_MUP_graph), aes(x=CoverCrop, y=emmean, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(21,22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)+10), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("Enzyme-labile P") +ylab(expression(paste("(",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  theme(legend.position="none")+xlab("")+
  ylim( min(modenzyme_stable_MUP_graph$lower.CL),max(modenzyme_stable_MUP_graph$upper.CL)+max(modenzyme_stable_MUP_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

##### Model fit

```{r, echo=F}
formula(modenzyme_stable_MUP)

anova(modenzyme_stable_MUP) %>%
  kable(caption= "Total enzyme-stable organic P: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r, echo=F}  
########  + estimating means P figure ####################

#cld from full_model (for Pi)
Pi_cld <-as.data.frame(multcomp::cld(emmeans(EAA_Pi_model.full, ~CoverCrop, adjust="bon", nesting=NULL),Letters=letters, type = "response"))
modEAA_Pi_cld<-Pi_cld

#cld from full_model (for stable)
Enzyme_stable_P_cld <-as.data.frame(multcomp::cld(emmeans(EAA_enzyme_stable_MUP_model.full, ~CoverCrop, adjust="bon", nesting=NULL),Letters=letters, type = "response"))
modEAA_stable_cld<-Enzyme_stable_P_cld

#cld from red_model (for labile)
modEAA_labile_red_cld <-as.data.frame(multcomp::cld(emmeans(modtotal_enzyme_labile_P, ~CoverCrop, adjust="bon", nesting=NULL),Letters=letters, type = "response"))
modEAA_labile_red_cld2<- modEAA_labile_red_cld[,c(1,7)]

modEAA_labile_full <-as.data.frame(multcomp::cld(emmeans(EAA_total_enzyme_labile_P_model.full, ~CoverCrop, adjust="bon", nesting=NULL),type="response"))
modEAA_labile_full$.group<-NULL
colnames(modEAA_labile_full)[which(names(modEAA_labile_full) == "response")] <- "emmean"

modEAA_labile_graph<-merge(modEAA_labile_red_cld2, modEAA_labile_full)

modEAA_labile_graph$.group<-stri_replace_all_fixed(modEAA_labile_graph$.group, " ", "")

modEAA_labile_cld<-modEAA_labile_graph

# generating dataframe with Pi, enzyme labile and enzyme-stable organic P
modEAA_Pi_cld$type <-"inorganic P"
modEAA_labile_cld$type <-"enzyme labile Porg"
modEAA_stable_cld$type <-"enzyme stable Porg"

modEAA_Pi_cld$.group<-stri_replace_all_fixed(modEAA_Pi_cld$.group, " a ", "")
modEAA_stable_cld$.group<-stri_replace_all_fixed(modEAA_stable_cld$.group, " a ", "")

modEAA_Pi_cld<-arrange(modEAA_Pi_cld, CoverCrop)
modEAA_labile_cld<-arrange(modEAA_labile_cld, CoverCrop)
modEAA_stable_cld<-arrange(modEAA_stable_cld, CoverCrop)

modEAA_stable_cld$upper.CL <-modEAA_stable_cld$upper.CL+modEAA_labile_cld$emmean
modEAA_stable_cld$lower.CL <-modEAA_stable_cld$lower.CL+modEAA_labile_cld$emmean

modEAA_Pi_cld$upper.CL <-modEAA_Pi_cld$upper.CL+modEAA_stable_cld$emmean+modEAA_labile_cld$emmean
modEAA_Pi_cld$lower.CL <-modEAA_Pi_cld$lower.CL+modEAA_stable_cld$emmean+modEAA_labile_cld$emmean

EAA_merged_2<-rbind(modEAA_labile_cld,modEAA_stable_cld, modEAA_Pi_cld)
EAA_merged_2$type<-as.factor(EAA_merged_2$type)
EAA_merged_2$type<-factor(EAA_merged_2$type, levels=c("inorganic P", "enzyme stable Porg","enzyme labile Porg" ))

###### + stacked bar chart ########## 
library(cowplot)

EAA_merged <- ggplot(EAA_merged_2, aes(x=CoverCrop, y=emmean, fill=type))+
  geom_bar(stat="identity", colour="black",width=0.7)+
  geom_text(aes(label=.group, x=CoverCrop, y=lower.CL-lower.CL/10), size=6,color="black")+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), position=position_dodge(0.3),width=.2)+
  #theme(legend.position = "none")+
  ggtitle("Soil P pools \n(inorganic, stable and enzyme-labile P)") +
  ylab(expression(paste("(",mu, "g  P ",g^-1, " soil)",sep="")))+ xlab("")+
  Motheme_EAA+scale_fill_manual(values=c( "#FEE5D9", "#FCAE91", "#FB6A4A"))
EAA_merged
```

``` {r,  include=F}
png("EAA_P.png",width=16,height=16, units = "cm", res=800)
EAA_merged+theme(legend.position = "none")
dev.off()

svg("EAA_P.svg",width=8,height=8)
EAA_merged+theme(legend.position = "none")
dev.off()

svg("EAA_P_Legend.svg",width=8,height=8)
plot(get_legend(EAA_merged))
dev.off()
```

\newpage
## Fig 3b1 Phytase-labile P

``` {r Fig 3b1,fig.width = 4, fig.height=4,echo=F}  
######## 3.6 EAA_PHYF ########
ggplot(dataAll_EAA , aes(x=CoverCrop, y=PHYF, fill=CoverCrop))+
  geom_boxplot()+
  geom_text(data = dataAll_EAA, aes(x = CoverCrop, group=CoverCrop, y = PHYF, label = Plot), position=position_dodge(1), hjust=.5)+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
EAA_PHYF_model.full<- lm(PHYF~CoverCrop+Block, na.action=na.omit, data=dataAll_EAA)
anova(EAA_PHYF_model.full)

EAA_PHYF_model.Type<- lm(PHYF~Type+Block, na.action=na.omit, data=dataAll_EAA)
anova(EAA_PHYF_model.Type)

# Cover Crop not sig! But Type tended (p=0.6)

multcomp::cld(emmeans(lm(PHYF~Type+Block, na.action=na.omit, data=dataAll_EAA),~Type))
# change between rhizosphere and bulk soil

modPHYF<-EAA_PHYF_model.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modPHYF),resid(modPHYF),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modPHYF)) 
qqline(y=resid(modPHYF))
# histogram of residuals
hist(resid(modPHYF))

shapiro.test(resid(modPHYF))

modPHYF_graph <-as.data.frame(multcomp::cld(emmeans(EAA_PHYF_model.full, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
```

```{r, fig.width = 6, echo=F }

PHYF_graph <- ggplot(subset(modPHYF_graph), aes(x=CoverCrop, y=emmean, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(21,22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)+10), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("Fungal phytase-labile P") +
  ylab(expression(paste("(",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  theme(legend.position="none")+xlab("")+
  ylim( min(modPHYF_graph$lower.CL),max(modPHYF_graph$upper.CL)+max(modPHYF_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
PHYF_graph
```

##### Model fit 

```{r, echo=F}
formula(modPHYF)

anova(modPHYF) %>%
  kable(caption= "Fungal phytase-labile organic P: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
### Bacterial phytase-labile P

``` {r Bacterial Phytase,fig.width = 4, fig.height=4,echo=F}  
######## 3.7 EAA_PHYB ########
ggplot(dataAll_EAA , aes(x=CoverCrop, y=PHYB, fill=CoverCrop))+
  geom_boxplot()+
  geom_text(data = dataAll_EAA, aes(x = CoverCrop, group=CoverCrop, y = PHYB, label = Plot), position=position_dodge(1), hjust=.5)+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
EAA_PHYB_model.full<- lm((PHYB)~CoverCrop+Block, na.action=na.omit, data=dataAll_EAA)
anova(EAA_PHYB_model.full)

modPHYB<-EAA_PHYB_model.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modPHYB),resid(modPHYB),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modPHYB)) 
qqline(y=resid(modPHYB))
# histogram of residuals
hist(resid(modPHYB))

shapiro.test(resid(modPHYB))

modPHYB_graph <-as.data.frame(multcomp::cld(emmeans(modPHYB, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

```

```{r, fig.width = 6, echo=F }

PHYB_graph <- ggplot(modPHYB_graph, aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(21,22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)+10), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("Enzyme-labile P") +ylab(expression(paste("(",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  theme(legend.position="none")+xlab("")+
  ylim( min(modPHYB_graph$lower.CL),max(modPHYB_graph$upper.CL)+max(modPHYB_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
PHYB_graph
```

##### Model fit 

```{r, echo=F}
formula(modPHYB)

anova(modPHYB) %>%
  kable(caption= "Bacterial phytase-labile organic P: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig 3b2 Phosphomonoesterase-labile P

``` {r Fig 3b2,fig.width = 4, fig.height=4,echo=F}  
######## 3.8 EAA_GP ########
ggplot(dataAll_EAA , aes(x=CoverCrop, y=GP, fill=CoverCrop))+
  geom_boxplot()+
  geom_text(data = dataAll_EAA, aes(x = CoverCrop, group=CoverCrop, y = GP, label = Plot), position=position_dodge(1), hjust=.5)+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
EAA_GP_model.full<- lm(GP~CoverCrop+Block, na.action=na.omit, data=dataAll_EAA)
anova(EAA_GP_model.full)

EAA_GP_model.Type<- lm(GP~Type+Block, na.action=na.omit, data=dataAll_EAA)
anova(EAA_GP_model.Type)

anova(EAA_GP_model.full,EAA_GP_model.Type)
# EAA_GP_model.full has lower RSS

modGP<-EAA_GP_model.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modGP),resid(modGP),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modGP)) 
qqline(y=resid(modGP))
# histogram of residuals
hist(resid(modGP))

shapiro.test(resid(modGP))

modGP_graph <-as.data.frame(multcomp::cld(emmeans(modGP, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))

#multcomp::cld(emmeans(lm(GP~Type+Block, na.action=na.omit, data=dataAll_EAA), ~Type, adjust="bon"),Letters=letters, type = "response")
# 63.1/44.5 #increae of GP in rhizosphere
```

```{r, fig.width = 6, echo=F }

GP_graph <- ggplot(subset(modGP_graph), aes(x=CoverCrop, y=emmean, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(21,22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("Phosphomonoesterase-labile P") +ylab(expression(paste("(",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  theme(legend.position="none")+xlab("")+
  ylim( min(modGP_graph$lower.CL),max(modGP_graph$upper.CL)+max(modGP_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
GP_graph

```

##### Model fit

```{r, echo=F}
formula(modGP)

anova(modGP) %>%
  kable(caption= "Phosphomonoesterase-labile P: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig 3b3 Phosphodiesterase-labile P

``` {r Fig 3b3,fig.width = 4, fig.height=4,echo=F}  
######## 3.9 EAA_DNA (diesterase-labile Porg) ########
ggplot(dataAll_EAA , aes(x=CoverCrop, y=DNA, fill=CoverCrop))+
  geom_boxplot()+
  geom_text(data = dataAll_EAA, aes(x = CoverCrop, group=CoverCrop, y = DNA, label = Plot), position=position_dodge(1), hjust=.5)+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
EAA_DNA_model.full<- lm(DNA~CoverCrop+Block, na.action=na.omit, data=dataAll_EAA)
anova(EAA_DNA_model.full)

EAA_DNA_model.Type<- lm(DNA~Type+Block, na.action=na.omit, data=dataAll_EAA)

anova(EAA_DNA_model.Type)

anova(EAA_DNA_model.full,EAA_DNA_model.Type )
#RSS not sig different

modDNA<-EAA_DNA_model.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F}
plot(fitted(modDNA),resid(modDNA),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modDNA)) 
qqline(y=resid(modDNA))
# histogram of residuals
hist(resid(modDNA))

shapiro.test(resid(modDNA))

modDNA_graph <-as.data.frame(multcomp::cld(emmeans(modDNA, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
```

```{r, fig.width = 6, echo=F }

DNA_graph <- ggplot(subset(modDNA_graph), aes(x=CoverCrop, y=emmean, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(21,22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)+10), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("Phosphodiesterase-labile P") +ylab(expression(paste("(",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  theme(legend.position="none")+xlab("")+
  ylim( min(modDNA_graph$lower.CL),max(modDNA_graph$upper.CL)+max(modDNA_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
DNA_graph
```

##### Model fit

```{r, echo=F}
formula(modDNA)

anova(modDNA) %>%
  kable(caption= "Phosphodiesterase-labile P: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r, include=F}  
########## + estimating means Enzyme-labile P####### 
# (with full model, calculating clds with reduced models and merging both)

modPHYF_full_cld <-as.data.frame(multcomp::cld(emmeans(EAA_PHYF_model.full, ~CoverCrop, adjust="bon", nesting=NULL),Letters=letters, type = "response"))

# in case of transformation, the estimated mean is stored in a column named "response" (whyever...)
colnames(modPHYF_full_cld)[which(names(modPHYF_full_cld) == "response")] <- "emmean"

#
modGP_red_cld <-as.data.frame(multcomp::cld(emmeans(modGP, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
modGP_red_cld2<- modGP_red_cld[,c(1,7)]

EAA_GP_model_full <-as.data.frame(multcomp::cld(emmeans(EAA_GP_model.full, ~CoverCrop, adjust="bon"), type="response"))
EAA_GP_model_full$.group<-NULL
colnames(EAA_GP_model_full)[which(names(EAA_GP_model_full) == "response")] <- "emmean"

modGP_graph<-merge(modGP_red_cld2, EAA_GP_model_full)

modGP_graph$.group<-stri_replace_all_fixed(modGP_graph$.group, " ", "")

#
modEAA_DNA_red_cld <-as.data.frame(multcomp::cld(emmeans(modDNA, ~CoverCrop, adjust="bon", nesting=NULL),Letters=letters, type = "response"))
modEAA_DNA_red_cld2<- modEAA_DNA_red_cld[,c(1,7)]

EAA_DNA_model_full <-as.data.frame(multcomp::cld(emmeans(EAA_DNA_model.full, ~CoverCrop, adjust="bon"), type="response"))
EAA_DNA_model_full$.group<-NULL
colnames(EAA_DNA_model_full)[which(names(EAA_DNA_model_full) == "response")] <- "emmean"

modEAA_DNA_graph<-merge(modEAA_DNA_red_cld2, EAA_DNA_model_full)

# emmeans of full models with clds from reduced models
modPHYF_cld <-modPHYF_full_cld
modGP_cld <-modGP_graph
modEAA_DNA_cld <-modEAA_DNA_graph

modPHYF_cld$type <-"phytase-labile P"
modGP_cld$type <-"monoesterase-labile P"
modEAA_DNA_cld$type <-"diesterase-labile P"

modPHYF_cld<-arrange(modPHYF_cld, CoverCrop)
modGP_cld<-arrange(modGP_cld, CoverCrop)
modEAA_DNA_cld<-arrange(modEAA_DNA_cld, CoverCrop)

# stacked bar chart requires to "lift" CLs
modEAA_DNA_cld$upper.CL <-modEAA_DNA_cld$upper.CL+modPHYF_cld$emmean+modGP_cld$emmean
modEAA_DNA_cld$lower.CL <-modEAA_DNA_cld$lower.CL+modPHYF_cld$emmean+modGP_cld$emmean

modGP_cld$upper.CL <-modGP_cld$upper.CL+modPHYF_cld$emmean
modGP_cld$lower.CL <-modGP_cld$lower.CL+modPHYF_cld$emmean

EAA_enzymes<-rbind(modPHYF_cld,modEAA_DNA_cld, modGP_cld)
EAA_enzymes$type<-as.factor(EAA_enzymes$type)
EAA_enzymes$type<-factor(EAA_enzymes$type, levels=c("diesterase-labile P","monoesterase-labile P","phytase-labile P" ))
```
+ stacked bar chart: Enzymatic availability of Porg
```{r, echo=F}
EAA_enzymes_graph <- ggplot(EAA_enzymes, aes(x=CoverCrop, y=emmean, fill=type))+
  geom_bar(stat="identity", colour="black",width=0.7)+
  geom_text(aes(label=.group, x=CoverCrop, y=lower.CL-10), size=6,color="black")+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), position=position_dodge(0.4),width=.2)+
  # theme(legend.position="none")+
  ggtitle("Enzymatically available soil organic P") +ylab(expression(paste("(",mu,"g ",P[org]," ", g^-1," soil)", sep="")))+ xlab("")+
  Motheme_EAA+
  scale_fill_manual(values=c( "#E7D4E8","#C2A5CF","#9970AB"))
EAA_enzymes_graph
```

``` {r,  include=F}
# saving graph with and without legend
png("EAA_enzymelabile.png",width=16,height=16, units = "cm", res=800)
EAA_enzymes_graph+theme(legend.position = "none")
dev.off()

svg("EAA_enzymelabile.svg",width=6,height=6)
EAA_enzymes_graph+theme(legend.position = "none")
dev.off()

svg("EAA_enzymelabile_Legend.svg",width=8,height=8)
plot(get_legend(EAA_enzymes_graph))
dev.off()


EAA_enzymes_leg<- get_legend(EAA_enzymes_graph+theme(legend.text=element_text(size=10)))
EAA_merged_leg<- get_legend(EAA_merged+theme(legend.text=element_text(size=10)))

# creating panel with graphs and legend
plot2by2 <- plot_grid( EAA_merged+theme(legend.position = "none"),
                       EAA_enzymes_graph+theme(legend.position = "none"), 
                       EAA_merged_leg, 
                       EAA_enzymes_leg,
                       ncol = 2, nrow=2, 
                       rel_widths = c(16, 16,8,8), rel_heights = c(16, 16,8,8) )
plot2by2

save_plot("EAA_all.svg", plot2by2,
          ncol = 2, # we're saving a grid plot of 2 columns
          nrow = 2, # and 2 rows
          # each individual subplot should have an aspect ratio of 1.3
          base_aspect_ratio = 1, base_height=6, base_width=6
)

```

\newpage

# Fig 4 Microbial P

``` {r Fig 4,fig.width = 8, fig.height=4, echo=F}  

box_Pmic <-ggplot(dataAll, aes(x=CoverCrop:Type, y=Pmic, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(",mu, "g", P[mic], " ",g^-1,")")))+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = Pmic, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ggtitle("Microbial P")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
box_Pmic
```

``` {r, include=F}
modPmic.full<-
  lmer(sqrt(Pmic)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modPmic.full)
VarCorr(modPmic.full)
summary(modPmic.full)
vcov(modPmic.full)

modPmic.red<-
  lmer(sqrt(Pmic)~Type*CoverCrop*Month+
         Block+(1|TypeID), na.action=na.omit, data=dataAll)
VarCorr(modPmic.red)
anova(modPmic.red)
multcomp::cld(emmeans(modPmic.red, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response")

# CoverCrop not sig! Type sig!

anova(modPmic.full,modPmic.red)

#full has lowest AIC
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modPmic<-modPmic.red
plot(fitted(modPmic),resid(modPmic),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modPmic)) 
qqline(y=resid(modPmic))
# histogram of residuals
hist(resid(modPmic))

shapiro.test(resid(modPmic))
#cld(emmeans(modPmic, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modPmic_cld <-as.data.frame(multcomp::cld(emmeans(modPmic, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))

modPmic_cld$.group<-stri_replace_all_fixed(modPmic_cld$.group, " ", "")

modPmic_cld$rhizo<-revalue(modPmic_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
Pmic_graph <- ggplot(na.omit(modPmic_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)+1), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("")+
  ylab(expression(paste( "Microbial P (",mu, "g ", P[mic], " ",g^-1,")")))+
  #theme(legend.position="none")+
  xlab("")+
  ylim( NA,max(modPmic_cld$upper.CL)+max(modPmic_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
Pmic_graph

```


``` {r,  include=F}
svg("Fig4_Pmic.svg",width=12, height=8)
Pmic_graph
dev.off()

png("Fig4_Pmic.png",width=32,height=16, units = "cm", res=800)
Pmic_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modPmic.full)
anova(modPmic) %>%
  kable(caption= "Microbial P: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
# Fig 5 Microbial abundance and community structure

## Fig 5a Gram-positive bacteria

``` {r Fig 5a,fig.width = 8, fig.height=4, echo=F}  
######## 5a  Gram+ (Gpos) ########
ggplot(dataAll, aes(x=CoverCrop:Type, y=Gpos, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = Gpos, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ggtitle("Gram+ abundance in bulk soil")+
  ylab(expression(paste( "(nmol PLFA ",g^-1,")")))+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modGpos.full<-
  lmer((Gpos)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modGpos.full)
VarCorr(modGpos.full)
summary(modGpos.full)
vcov(modGpos.full)

modGpos.red<-
  lmer((Gpos)~Type*CoverCrop*Month+
         Block+(1|TypeID), na.action=na.omit, data=dataAll)

# CoverCrop not sig! Type:Month sig!

anova(modGpos.full,modGpos.red)

#red has lowest AIC

modGpos<-modGpos.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modGpos),resid(modGpos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modGpos)) 
qqline(y=resid(modGpos))
# histogram of residuals
hist(resid(modGpos))

shapiro.test(resid(modGpos))
#cld(emmeans(modGpos, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modGpos_cld <-as.data.frame(multcomp::cld(emmeans(modGpos, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))
#multcomp::cld(emmeans(modGpos, ~CoverCrop|Type|Month, adjust="bon"),Letters=letters, type = "response")

modGpos_cld$rhizo<-revalue(modGpos_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
Gpos_graph <- ggplot(na.omit(modGpos_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  ##geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("")+
  ylab(expression(paste( "Gram-positive bacteria (nmol PLFA ",g^-1,")")))+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(modGpos_cld$upper.CL)+max(modGpos_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
Gpos_graph

```


``` {r,  include=F}
svg("Fig5a_Gpos.svg",width=12, height=8)
Gpos_graph
dev.off()

png("Fig5a_Gpos.png",width=32,height=16, units = "cm", res=800)
Gpos_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modGpos.full)
anova(modGpos) %>%
  kable(caption= "Gram-positive abundance: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig 5b Gram-negative bacteria

``` {r Fig 5b,fig.width = 8, fig.height=4,echo=F}  
######## 5b Gram- (Gneg) ########

ggplot(dataAll, aes(x=CoverCrop:Type, y=Gneg, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = Gneg, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ggtitle("Gram- abundance in bulk soil")+
  ylab(expression(paste( "(nmol PLFA ",g^-1,")")))+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modGneg.full<-
  lmer((Gneg)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modGneg.full)
VarCorr(modGneg.full)
summary(modGneg.full)
vcov(modGneg.full)

modGneg.red<-
  lmer((Gneg)~Type*CoverCrop*Month+
         Block+(1|TypeID), na.action=na.omit, data=dataAll)
anova(modGneg.red)
multcomp::cld(emmeans(lmer(Paq~CoverCrop*Type+Month+
                             Block+(1|TypeID), na.action=na.omit, data=dataAll), 
                      ~CoverCrop|Type:Month, adjust="bon"),Letters=letters, type = "response")

# CoverCrop not sig! Type sig!


anova(modGneg.full,modGneg.red)

#red has lowest AIC
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modGneg<-modGneg.red
#modGneg<-modGneg.full_redstep
plot(fitted(modGneg),resid(modGneg),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modGneg)) 
qqline(y=resid(modGneg))
# histogram of residuals
hist(resid(modGneg))

shapiro.test(resid(modGneg))
#cld(emmeans(modGneg, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modGneg_cld <-as.data.frame(multcomp::cld(emmeans(modGneg, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))
#
#multcomp::cld(emmeans(modGneg, ~CoverCrop, adjust="bon"),Letters=letters, type = "response")

modGneg_cld$rhizo<-revalue(modGneg_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
Gneg_graph <- ggplot(na.omit(modGneg_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  ##geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("")+
  ylab(expression(paste( "Gram-negative bacteria (nmol PLFA ",g^-1,")")))+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(modGneg_cld$upper.CL)+max(modGneg_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
Gneg_graph

```


``` {r,  include=F}
svg("Fig5b_Gneg.svg",width=12, height=8)
Gneg_graph
dev.off()

png("Fig5b_Gneg.png",width=32,height=16, units = "cm", res=800)
Gneg_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modGneg.full)
anova(modGneg) %>%
  kable(caption= "Gram-negative abundance: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig 5c Fungal abundance

``` {r Fig 5c,fig.width = 8, fig.height=4,echo=F}  
######## 5c Fungi ########
ggplot(dataAll, aes(x=CoverCrop:Type, y=Fungi, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = Fungi, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ggtitle("Fungal abundance in bulk soil")+
  ylab(expression(paste( "(nmol PLFA ",g^-1,")")))+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modFungi.full<-
  lmer(log(Fungi)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modFungi.full)
VarCorr(modFungi.full)
summary(modFungi.full)
vcov(modFungi.full)
```

``` {r ,fig.width = 8, fig.height=4, echo=F} 
modFungi<-modFungi.full
#modFungi<-modFungi.full_redstep
plot(fitted(modFungi),resid(modFungi),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modFungi)) 
qqline(y=resid(modFungi))
# histogram of residuals
hist(resid(modFungi))

shapiro.test(resid(modFungi))

#cld(emmeans(modFungi, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modFungi_cld <-as.data.frame(multcomp::cld(emmeans(modFungi, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))
#multcomp::cld(emmeans(modFungi, ~CoverCrop|Type|Month, adjust="bon"),Letters=letters, type = "response")

modFungi_cld$rhizo<-revalue(modFungi_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
Fungi_graph <- ggplot(na.omit(modFungi_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=0.1, size=rel(6), color="black")+
  ggtitle("")+
  ylab(expression(paste( "Fungi (nmol PLFA ",g^-1,")")))+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(modFungi_cld$upper.CL)+max(modFungi_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
Fungi_graph
```

``` {r,  include=F}
svg("Fig5c_Fungi.svg",width=12, height=8)
Fungi_graph
dev.off()

png("Fig5c_Fungi.png",width=32,height=16, units = "cm", res=800)
Fungi_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modFungi.full)
anova(modFungi) %>%
  kable(caption= "Fungal abundance: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
# Fig 6 Potential enzymatic activity

## Fig 6a Acid phosphomonoesterase activity

``` {r Fig 6a,fig.width = 8, fig.height=4, echo=F}  
####### Fig 6a ac.phos activity (acid phosphomonoesterase) ########
ggplot(dataAll, aes(x=CoverCrop:Type, y=ac.phos, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = ac.phos, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ylab(expression(paste("(nmol MUB ",g^-1,h^-1,")")))+
  ggtitle("Acid phosphomonoesterase activity under\n cover crops")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modac.phos.full<-
  lmer(log(ac.phos)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modac.phos.full)
VarCorr(modac.phos.full)
summary(modac.phos.full)
vcov(modac.phos.full)

# Type p<0.0001 Type sig!

anova(modac.phos.full)

modac.phos.red<-
  lmer((ac.phos)~Type+CoverCrop*Month+
         Block+(1|TypeID), na.action=na.omit, data=dataAll)

anova(modac.phos.red)
anova(modac.phos.full,modac.phos.red )
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modac.phos<-modac.phos.full
plot(fitted(modac.phos),resid(modac.phos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modac.phos)) 
qqline(y=resid(modac.phos))
# histogram of residuals
hist(resid(modac.phos))

shapiro.test(resid(modac.phos))
#cld(emmeans(modac.phos, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modac.phos_cld <-as.data.frame(multcomp::cld(emmeans(modac.phos, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))
#multcomp::cld(emmeans(modac.phos, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response")


modac.phos_cld$rhizo<-revalue(modac.phos_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
ac.phos_graph <- ggplot(na.omit(modac.phos_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ylab(expression(paste("Acid phosphomonoesterase activity (nmol MUB ",g^-1,h^-1,")")))+
  ggtitle(" ")+
  #theme(legend.position="none")+
  xlab("")+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
ac.phos_graph
```


``` {r,  include=F}
svg("modac.phos_graph.svg",width=12, height=8)
ac.phos_graph
dev.off()

png("modac.phos_graph.png",width=32,height=16, units = "cm", res=800)
ac.phos_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modac.phos.full)
anova(modac.phos) %>%
  kable(caption= "Acid phosphomonoesterase activity: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig 6b Alkaline phosphomonoesterase activity

``` {r Fig 6b,fig.width = 8, fig.height=4,echo=F}  
####### Fig 6b alk.phos (alkaline phosphomonoesterase) ########
ggplot(dataAll, aes(x=CoverCrop:Type, y=alk.phos, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = alk.phos, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ylab(expression(paste("(nmol MUB ",g^-1,h^-1,")")))+
  ggtitle("Alkaline phosphomonoesterase activity under\n cover crops")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modalk.phos.full<-
  lmer((alk.phos)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modalk.phos.full)
VarCorr(modalk.phos.full)
summary(modalk.phos.full)
vcov(modalk.phos.full)

modalk.phos.red<-
  lmer((alk.phos)~Type*CoverCrop*Month+
         Block+(1|TypeID), na.action=na.omit, data=dataAll)

# Type sig!

anova(modalk.phos.full, modalk.phos.red)

#red has lowest AIC
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modalk.phos<-modalk.phos.red
#modalk.phos<-modalk.phos.full_redstep
plot(fitted(modalk.phos),resid(modalk.phos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modalk.phos)) 
qqline(y=resid(modalk.phos))
# histogram of residuals
hist(resid(modalk.phos))

shapiro.test(resid(modalk.phos))
#cld(emmeans(modalk.phos, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modalk.phos_cld <-as.data.frame(multcomp::cld(emmeans(modalk.phos, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))

modalk.phos_cld$rhizo<-revalue(modalk.phos_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
alk.phos_graph <- ggplot(na.omit(modalk.phos_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ylab(expression(paste("Alk. phosphomonoesterase activity (nmol MUB ",g^-1,h^-1,")")))+
  ggtitle(" ")+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(modalk.phos_cld$upper.CL)+max(modalk.phos_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
alk.phos_graph
```


``` {r,  include=F}
svg("modalk.phos_graph.svg",width=12, height=8)
alk.phos_graph
dev.off()

png("modalk.phos_graph.png",width=32,height=16, units = "cm", res=800)
alk.phos_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modalk.phos.full)
anova(modalk.phos) %>%
  kable(caption= "Alkaline phosphomonoesterase activity: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig 6c Phosphodiesterase acticity

``` {r Fig 6c,fig.width = 8, fig.height=4,echo=F}  
####### Fig 6c diphos activity (phosphodiesterase) ########
ggplot(dataAll, aes(x=CoverCrop:Type, y=diphos, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = diphos, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ylab(expression(paste("(nmol MUB ",g^-1,h^-1,")")))+
  ggtitle("Phosphodiesterase activity under\n cover crops")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
moddiphos.full<-lmer(log(diphos)~Type*CoverCrop*Month+
                       Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(moddiphos.full)
VarCorr(moddiphos.full)
# Type:CoverCrop:Month p=0.076 Type:Month sig!

moddiphos.red<-lm(log(diphos)~Type*CoverCrop*Month+
                    Block, na.action=na.omit, data=dataAll)
anova(moddiphos.red)

anova(moddiphos.full,moddiphos.red)

#red has lowest AIC
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
moddiphos<-moddiphos.red
plot(fitted(moddiphos),resid(moddiphos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(moddiphos)) 
qqline(y=resid(moddiphos))
# histogram of residuals
hist(resid(moddiphos))

shapiro.test(resid(moddiphos))
moddiphos_cld <-as.data.frame(multcomp::cld(emmeans(moddiphos, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))

moddiphos_cld$rhizo<-revalue(moddiphos_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
diphos_graph <- ggplot(na.omit(moddiphos_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(5), color="black")+
  ylab(expression(paste("Phosphodiesterase activity (nmol MUB ",g^-1,h^-1,")")))+
  ggtitle(" ")+
  #theme(legend.position="none")+
  xlab("")+
  ylim( NA,420)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
diphos_graph
```

``` {r,  include=F}
svg("moddiphos_graph.svg",width=12, height=8)
diphos_graph
dev.off()

png("moddiphos_graph.png",width=32,height=16, units = "cm", res=800)
diphos_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(moddiphos.full)
anova(moddiphos) %>%
  kable(caption= "Phosphodiesterase activity: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig 6d N-acetyl-hexosaminidase  activity

``` {r Fig 6d,fig.width = 8, fig.height=4,echo=F}  
####### Fig 6d NAG (N-acetyl-hexosaminidase ) ########

ggplot(dataAll, aes(x=CoverCrop:Type, y=NAG, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = NAG, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ylab(expression(paste("(nmol MUB ",g^-1,h^-1,")")))+
  ggtitle("N-acetyl-hexosaminidase  activity ")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modNAG.full<-
  lmer(log(NAG)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modNAG.full)
VarCorr(modNAG.full)
summary(modNAG.full)
vcov(modNAG.full)

modNAG.red<-
  lm(log(NAG)~Type*CoverCrop*Month+
       Block, na.action=na.omit, data=dataAll)

# Type:CoverCrop:Month sig!

anova(modNAG.full, modNAG.red)

#full has lowest AIC

modNAG<-modNAG.full
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
plot(fitted(modNAG),resid(modNAG),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modNAG)) 
qqline(y=resid(modNAG))
# histogram of residuals
hist(resid(modNAG))

shapiro.test(resid(modNAG))
#cld(emmeans(modNAG, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modNAG_cld <-as.data.frame(multcomp::cld(emmeans(modNAG, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))
#multcomp::cld(emmeans(modNAG, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response")

modNAG_cld$rhizo<-revalue(modNAG_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
NAG_graph <- ggplot(na.omit(modNAG_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ylab(expression(paste("N-acetyl-hexosaminidase activity (nmol MUB ",g^-1,h^-1,")")))+
  ggtitle("")+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(modNAG_cld$upper.CL)+max(modNAG_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
NAG_graph
```


``` {r,  include=F}
svg("modNAG_graph.svg",width=12, height=8)
NAG_graph
dev.off()

png("modNAG_graph.png",width=32,height=16, units = "cm", res=800)
NAG_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modNAG)
anova(modNAG) %>%
  kable(caption= "N-acetyl-hexosaminidase  activity: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
# Fig 7 Correlations

## Fig 7a phoD vs alk. phosphomonoesterase

``` {r Fig 7a,fig.width = 4, fig.height=4, echo=F}  
####### Fig 7a: phoD_dw_x_alk.phos ########

#str(dataAll)

ggplot(dataAll,aes(y=log(phoD_dw),x=alk.phos, color=Type))+
  geom_point(aes(color=Type))+stat_smooth(method="lm")
```

```{r, include=F}
phoD_dw_x_alk.phos_model.full<- lmer(alk.phos~CoverCrop*Type*log(phoD_dw)+Block+(1|Plot), na.action=na.exclude, data=dataAll)
anova(phoD_dw_x_alk.phos_model.full)
VarCorr(phoD_dw_x_alk.phos_model.full)

phoD_dw_x_alk.phos_model.red.lm<- lm(alk.phos~Type*log(phoD_dw)+Block, na.action=na.exclude, data=dataAll)
anova(phoD_dw_x_alk.phos_model.red.lm)

phoD_dw_x_alk.phos_model.red.lm_step<-step(phoD_dw_x_alk.phos_model.red.lm, direction="both")
anova(phoD_dw_x_alk.phos_model.red.lm_step)
summary(phoD_dw_x_alk.phos_model.red.lm_step)

phoD_dw_x_alk.phos_model.figure<- lm(alk.phos~Type*log(phoD_dw), na.action=na.exclude, data=dataAll)
anova(phoD_dw_x_alk.phos_model.figure)

anova(phoD_dw_x_alk.phos_model.full, phoD_dw_x_alk.phos_model.red.lm, phoD_dw_x_alk.phos_model.figure,phoD_dw_x_alk.phos_model.red.lm_step)

# phoD_dw_x_alk.phos_model.red.lm_step has lowest AIC, but F-test results are similar. sticking to .red to conserve random effects
anova(phoD_dw_x_alk.phos_model.red.lm)
summary(phoD_dw_x_alk.phos_model.red.lm)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modphoD_dw_x_alk.phos<-phoD_dw_x_alk.phos_model.red.lm
plot(fitted(modphoD_dw_x_alk.phos),resid(modphoD_dw_x_alk.phos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues
qqnorm(resid(modphoD_dw_x_alk.phos)) 
qqline(y=resid(modphoD_dw_x_alk.phos))
# histogram of residuals
hist(resid(modphoD_dw_x_alk.phos))

shapiro.test(resid(modphoD_dw_x_alk.phos))
```

##### Model fit 

```{r, echo=F}
formula(modphoD_dw_x_alk.phos)

anova(modphoD_dw_x_alk.phos) %>%
  kable(caption= "Relation between alk.phosphomonoesterase activity and phoD per g soil") %>%
  kable_styling()
```

``` {r,  echo=F}
library(cowplot)
detach(package:cowplot)

# representing data with regression lines
phoD_dw_x_alk.phos_graph <-ggplot(dataAll,aes(y=log(phoD_dw),x=alk.phos))+
  geom_point(size=3, aes(shape=Type,fill=Type))+
  geom_smooth(method = "lm", se=T, color="black",linetype="dashed",
              formula = my.formula, size=1.5) +
  stat_poly_eq(formula = my.formula, label.x=0.9,size=6,
               aes(label = paste(..rr.label.., ..p.value.label..,sep = "~~~")), 
               parse = TRUE,
               rr.digits = 2,
               p.digits = 2) + 
  #ggtitle("Alkaline phosphomonoesterase activity vs phoD abundance") +
  scale_linetype_manual("Soil compartment",values=c(1,8))+
  scale_fill_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_shape_manual("Soil compartment",values=c(21,23))+  
  ylab(expression(paste("phoD (ln[copies ",  g^-1,"soil])",sep="")))+ labs(fill="")+
  xlab(expression(paste("Alk. phomonoesterase activity (nmol MUB ",g^-1,h^-1,")",sep="")))+ labs(fill="")+
  Motheme_all2
phoD_dw_x_alk.phos_graph
```

``` {r,  include=F}
svg("Fig8a.svg",width=8,height=8)
phoD_dw_x_alk.phos_graph
dev.off()

png("Fig8a.png",width=12,height=12, units = "cm", res=800)
phoD_dw_x_alk.phos_graph
dev.off()
```

\newpage
## Fig 7b: phoD vs bacterial PLFA

``` {r Fig 7b,fig.width = 4, fig.height=4,echo=F}  

#str(dataAll)
ggplot(dataAll_T1,aes(y=log(phoD_dw),x=Bacteria, color=Type))+
  geom_point(aes(color=Type))+stat_smooth(method="lm")
```

```{r, include=F}
phoD_dw_x_Bacteria_model.full<- lmer(log(phoD_dw)~CoverCrop*Type*Bacteria+Block+(1|Plot), na.action=na.exclude, data=dataAll)
anova(phoD_dw_x_Bacteria_model.full)
VarCorr(phoD_dw_x_Bacteria_model.full)

phoD_dw_x_Bacteria_model.red.lm<- lm(log(phoD_dw)~Type*Bacteria+Block, na.action=na.exclude, data=dataAll)
anova(phoD_dw_x_Bacteria_model.red.lm)
#phoD_dw_x_Bacteria_model.red.lm<- lm(Bacteria~Type*log(phoD_dw)+Type+Block, na.action=na.exclude, data=dataAll)

phoD_dw_x_Bacteria_model.red.lm_step<-step(phoD_dw_x_Bacteria_model.red.lm, direction="both")
anova(phoD_dw_x_Bacteria_model.red.lm_step)
summary(phoD_dw_x_Bacteria_model.red.lm_step)

anova(phoD_dw_x_Bacteria_model.full, phoD_dw_x_Bacteria_model.red.lm, phoD_dw_x_Bacteria_model.red.lm_step)

# phoD_dw_x_Bacteria_model.red.lm has lowest AIC
anova(phoD_dw_x_Bacteria_model.red.lm)
summary(phoD_dw_x_Bacteria_model.red.lm)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modphoD_dw_x_Bacteria<-phoD_dw_x_Bacteria_model.red.lm
#modphoD_dw<-phoD_dw_model.full
plot(fitted(modphoD_dw_x_Bacteria),resid(modphoD_dw_x_Bacteria),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modphoD_dw_x_Bacteria)) 
qqline(y=resid(modphoD_dw_x_Bacteria))
# histogram of residuals
hist(resid(modphoD_dw_x_Bacteria))

shapiro.test(resid(modphoD_dw_x_Bacteria))
```

##### Model fit

```{r, echo=F}
formula(modphoD_dw_x_Bacteria)

anova(modphoD_dw_x_Bacteria) %>%
  kable(caption= "Relation between phoD per g soil and bacterial abundance") %>%
  kable_styling()
```

``` {r, echo=F}
library(cowplot)
detach(package:cowplot)

# representing data with regression lines
phoD_dw_x_Bacteria_graph <-ggplot(dataAll,aes(y=log(phoD_dw),x=Bacteria))+
  geom_point(size=3, aes(shape=Type,fill=Type))+
  geom_smooth(method = "lm", se=T, color="black",linetype="dashed",
              formula = my.formula, size=1.5) +
  stat_poly_eq(formula = my.formula, label.x=0.9,size=6,
               aes(label = paste(..rr.label.., ..p.value.label..,sep = "~~~")), 
               parse = TRUE,
               rr.digits = 2,
               p.digits = 2) + 
  #ggtitle("Bacterial abundance vs phoD abundance") +
  scale_fill_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_shape_manual("Soil compartment",values=c(21,23))+
  ylab(expression(paste("phoD (ln[copies ",  g^-1,"soil])",sep="")))+ labs(fill="")+
  xlab(expression(paste("Bacteria (ng PLFA",  g^-1,"soil)",sep="")))+ labs(fill="")+
  Motheme_all2
phoD_dw_x_Bacteria_graph
```

``` {r,  include=F}
svg("Fig8b.svg",width=8,height=8)
phoD_dw_x_Bacteria_graph
dev.off()

png("Fig8b.png",width=12,height=12, units = "cm", res=800)
phoD_dw_x_Bacteria_graph
dev.off()
```

\newpage
## Fig 7c: Bacteria vs alk. phosphatase 

``` {r Fig 7c,fig.width = 4, fig.height=4,echo=F}  
#str(dataAll)

ggplot(dataAll,aes(y=alk.phos,x=log(Bacteria), color=Type))+
  geom_point(aes(color=Type))+stat_smooth(method="lm")
```

```{r, include=F}
Bacteria_x_alk.phos_model.full<- lmer(alk.phos~Date*Type*Bacteria+Block+(1|Plot), na.action=na.exclude, data=dataAll)
anova(Bacteria_x_alk.phos_model.full)
VarCorr(Bacteria_x_alk.phos_model.full)

Bacteria_x_alk.phos_model.red.lm<- lm(alk.phos~Date*Type*Bacteria+Block, na.action=na.exclude, data=dataAll)
anova(Bacteria_x_alk.phos_model.red.lm)

Bacteria_x_alk.phos_model.red.lm2<- lm(alk.phos~Bacteria+Block, na.action=na.exclude, data=dataAll)
anova(Bacteria_x_alk.phos_model.red.lm2)

Bacteria_x_alk.phos_model.step<-get_model(step(Bacteria_x_alk.phos_model.full, direction="both", keep="Block"))
anova(Bacteria_x_alk.phos_model.step)
summary(Bacteria_x_alk.phos_model.step)

anova(Bacteria_x_alk.phos_model.full, Bacteria_x_alk.phos_model.red.lm, Bacteria_x_alk.phos_model.red.lm2, Bacteria_x_alk.phos_model.step)

# Bacteria_x_alk.phos_model.step has lowest AIC
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modBacteria_x_alk.phos<-Bacteria_x_alk.phos_model.step
#modBacteria<-Bacteria_model.full
plot(fitted(modBacteria_x_alk.phos),resid(modBacteria_x_alk.phos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modBacteria_x_alk.phos)) 
qqline(y=resid(modBacteria_x_alk.phos))
# histogram of residuals
hist(resid(modBacteria_x_alk.phos))

shapiro.test(resid(modBacteria_x_alk.phos))
```

##### Model fit 

```{r, echo=F}
formula(modBacteria_x_alk.phos)

anova(modBacteria_x_alk.phos) %>%
  kable(caption= "Relation between Bacterial abundance and alk.phosphomonoesterase activity") %>%
  kable_styling()
```

``` {r, echo=F}
library(cowplot)
detach(package:cowplot)

# representing data with regression lines
Bacteria_x_alk.phos_graph<-ggplot(dataAll,aes(y=alk.phos,x=Bacteria))+
  geom_point(size=3, aes(shape=Type, fill=Type))+
  geom_smooth(method = "lm", se=T, color="black",linetype="dashed",
              formula = my.formula, size=1.5) +
  stat_poly_eq(formula = my.formula, label.x=0.9,size=6,
               aes(label = paste(..rr.label.., ..p.value.label..,sep = "~~~")), 
               parse = TRUE,
               rr.digits = 2,
               p.digits = 2) + 
  #ggtitle("Bacterial abundance vs alkaline phosphomonoesterase activity") +
  scale_fill_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_shape_manual("Soil compartment",values=c(21,23))+
  xlab(expression(paste("Bacteria (ng PLFA",  g^-1,"soil)",sep="")))+ labs(fill="")+
  ylab(expression(paste("Alk. phosphomonoesterase activity (nmol MUB ",g^-1,h^-1,")",sep="")))+ labs(fill="")+
  Motheme_all2
Bacteria_x_alk.phos_graph
```

``` {r,  include=F}
svg("Fig8c.svg",width=8,height=8)
Bacteria_x_alk.phos_graph
dev.off()

png("Fig8c.png",width=12,height=12, units = "cm", res=800)
Bacteria_x_alk.phos_graph+  theme(legend.position="none")
dev.off()
```


\newpage
## Fig 7d: phoD vs 16S

``` {r Fig 7d,fig.width = 4, fig.height=4,echo=F}  
#str(dataAll)
ggplot(dataAll_T1,aes(y=log(phoD_dw),x=log(bact_16S_dw), color=Type))+
  geom_point(aes(color=Type))+stat_smooth(method="lm")
```

```{r, include=F}
phoD_dw_x_bact_16S_dw_model.full<- lmer(log(phoD_dw)~CoverCrop*Type*log(bact_16S_dw)+Block+(1|Plot), na.action=na.exclude, data=dataAll)
anova(phoD_dw_x_bact_16S_dw_model.full)
VarCorr(phoD_dw_x_bact_16S_dw_model.full)

phoD_dw_x_bact_16S_dw_model.full.lm<- lm(log(phoD_dw)~Type*CoverCrop*log(bact_16S_dw)+Block, na.action=na.exclude, data=dataAll)
anova(phoD_dw_x_bact_16S_dw_model.full.lm)
#phoD_dw_x_bact_16S_dw_model.red.lm<- lm(bact_16S_dw~Type*log(phoD_dw)+Type+Block, na.action=na.exclude, data=dataAll)

phoD_dw_x_bact_16S_dw_model.red.lm_step<-step(phoD_dw_x_bact_16S_dw_model.full.lm, direction="both")
anova(phoD_dw_x_bact_16S_dw_model.red.lm_step)
summary(phoD_dw_x_bact_16S_dw_model.red.lm_step)

anova(phoD_dw_x_bact_16S_dw_model.full, phoD_dw_x_bact_16S_dw_model.red.lm_step)

# phoD_dw_x_bact_16S_dw_model.red.lm_step has lowest AIC
anova(phoD_dw_x_bact_16S_dw_model.red.lm_step)
summary(phoD_dw_x_bact_16S_dw_model.red.lm_step)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modphoD_dw_x_bact_16S_dw<-phoD_dw_x_bact_16S_dw_model.red.lm_step
#modphoD_dw<-phoD_dw_model.full
plot(fitted(modphoD_dw_x_bact_16S_dw),resid(modphoD_dw_x_bact_16S_dw),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modphoD_dw_x_bact_16S_dw)) 
qqline(y=resid(modphoD_dw_x_bact_16S_dw))
# histogram of residuals
hist(resid(modphoD_dw_x_bact_16S_dw))

shapiro.test(resid(modphoD_dw_x_bact_16S_dw))
```

##### Model fit 

```{r, echo=F}
formula(modphoD_dw_x_bact_16S_dw)

anova(modphoD_dw_x_bact_16S_dw) %>%
  kable(caption= "Relation between phoD per g soil and bacterial 16S abundance") %>%
  kable_styling()
```

``` {r,  echo=F}
library(cowplot)
detach(package:cowplot)

# representing data with regression lines
phoD_dw_x_bact_16S_dw_graph <-ggplot(dataAll_T1,aes(y=log(phoD_dw),x=log(bact_16S_dw), color=CoverCrop, shape=CoverCrop))+
  geom_point(size=3, aes(fill=CoverCrop, shape=CoverCrop), color="black")+
  geom_smooth(method = "lm", se=F, aes(color=CoverCrop),linetype="dashed",
              formula = my.formula, size=1.5) +
  stat_poly_eq(formula = my.formula, label.x=0.9,size=6,
               aes(label = paste(..rr.label.., ..p.value.label..,sep = "~~~")), 
               parse = TRUE,
               rr.digits = 2,
               p.digits = 2) + 
  #ggtitle("16S abundance vs phoD abundance") +
  #scale_fill_manual("Cover Crop",values=c("#a6611a","#018571"))+
  scale_color_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))+
  scale_shape_manual("Cover Crop",values=c(21,23, 22, 24))+
  ylab(expression(paste("phoD (ln[copies ",  g^-1,"soil])",sep="")))+ labs(fill="")+
  xlab(expression(paste("16S (ln[copies ",  g^-1,"soil])",sep="")))+ labs(fill="")+
  Motheme_all2
phoD_dw_x_bact_16S_dw_graph
```

``` {r,  include=F}
svg("Fig8d.svg",width=8,height=8)
phoD_dw_x_bact_16S_dw_graph
dev.off()

png("Fig8d.png",width=12,height=12, units = "cm", res=800)
phoD_dw_x_bact_16S_dw_graph
dev.off()
```

\newpage
# Fig 8 labile organic P vs enzyme activities
## Fig 8a: phosphomonoesterase-labile Porg vs ac. phosphomonoesterase activity

``` {r Fig 8a,fig.width = 4, fig.height=4, echo=F}  
#str(dataAll_EAA)

ggplot(dataAll_EAA,aes(y=GP,x=ac.phos, color=))+
  geom_point(aes(color=))+stat_smooth(method="lm")
```

```{r, include=F}
#GP_x_ac.phos_model.full<- lmer(ac.phos~CoverCrop*Type*GP+Block+(1|Plot), na.action=na.exclude, data=dataAll_EAA)
#anova(GP_x_ac.phos_model.full)
#VarCorr(GP_x_ac.phos_model.full)

GP_x_ac.phos_model.red.lm<- lm(ac.phos~Type*CoverCrop*GP+Block, na.action=na.exclude, data=dataAll_EAA)
anova(GP_x_ac.phos_model.red.lm)

GP_x_ac.phos_model.red.lm_step<-step(GP_x_ac.phos_model.red.lm, direction="both")
anova(GP_x_ac.phos_model.red.lm_step)
summary(GP_x_ac.phos_model.red.lm_step)

anova(GP_x_ac.phos_model.red.lm_step, GP_x_ac.phos_model.red.lm)

# No improvement of fit with step
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modGP_x_ac.phos<-GP_x_ac.phos_model.red.lm
plot(fitted(modGP_x_ac.phos),resid(modGP_x_ac.phos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues
qqnorm(resid(modGP_x_ac.phos)) 
qqline(y=resid(modGP_x_ac.phos))
# histogram of residuals
hist(resid(modGP_x_ac.phos))

shapiro.test(resid(modGP_x_ac.phos))

```

##### Model fit 

```{r, echo=F}
formula(modGP_x_ac.phos)

anova(modGP_x_ac.phos) %>%
  kable(caption= "Relation between ac.phosphomonoesterase activity and phoD per g soil") %>%
  kable_styling()
```

``` {r, echo=F}
library(cowplot)
detach(package:cowplot)

# representing data with regression lines
GP_x_ac.phos_graph <-ggplot(dataAll_EAA,aes(y=GP,x=ac.phos))+
  geom_point(size=3, aes(shape=Type,fill=Type))+
  geom_smooth(method = "lm", se=T, color="black",linetype="dashed",
              formula = my.formula, size=1.5) +
  stat_poly_eq(formula = my.formula, label.x=0.9,size=6,
               aes(label = paste(..rr.label.., ..p.value.label..,sep = "~~~")), 
               parse = TRUE,
               rr.digits = 2,
               p.digits = 2) + 
  #ggtitle("Acid phosphomonoesterases activity vs \nphosphomonoesterases-labile organic P") +
  scale_linetype_manual("Soil compartment",values=c(1,8))+
  scale_fill_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_shape_manual("Soil compartment",values=c(21,23))+  
  ylab(expression(paste("Phosphomonoesterase-labile P (",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  xlab(expression(paste("Acid phosphomonoesterase activity (nmol MUB ",g^-1,h^-1,")",sep="")))+ labs(fill="")+
  Motheme_all2
GP_x_ac.phos_graph
```

``` {r,  include=F}
svg("Fig9a.svg",width=8,height=8)
GP_x_ac.phos_graph
dev.off()

png("Fig9a.png",width=12,height=12, units = "cm", res=800)
GP_x_ac.phos_graph
dev.off()
```

\newpage
## Fig 8b: Phosphodiesterase-labile Porg vs phosphodiesterase activity
``` {r Fig 8b,fig.width = 4, fig.height=4,echo=F}  

ggplot(dataAll_EAA,aes(y=DNA,x=diphos, color=))+
  geom_point(aes(color=))+stat_smooth(method="lm")
```

```{r, include=F}
DNA_x_diphos_model.red.lm<- lm(diphos~Type*CoverCrop*DNA+Block, na.action=na.exclude, data=dataAll_EAA)
anova(DNA_x_diphos_model.red.lm)

#DNA_x_diphos_model.red.lm_step<-step(DNA_x_diphos_model.red.lm, direction="both")
#anova(DNA_x_diphos_model.red.lm_step)
#summary(DNA_x_diphos_model.red.lm_step)
#step seems not able to improve fit

#anova(DNA_x_diphos_model.red.lm_step, DNA_x_diphos_model.red.lm)

# DNA_x_diphos_model.red.lm has lowest AIC
anova(DNA_x_diphos_model.red.lm)
summary(DNA_x_diphos_model.red.lm)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modDNA_x_diphos<-DNA_x_diphos_model.red.lm
plot(fitted(modDNA_x_diphos),resid(modDNA_x_diphos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues
qqnorm(resid(modDNA_x_diphos)) 
qqline(y=resid(modDNA_x_diphos))
# histogram of residuals
hist(resid(modDNA_x_diphos))

shapiro.test(resid(modDNA_x_diphos))
```

##### Model fit 

```{r, echo=F}
formula(modDNA_x_diphos)

anova(modDNA_x_diphos) %>%
  kable(caption= "Relation between phosphodiesterasse activity and phoD per g soil") %>%
  kable_styling()
```

``` {r, echo=F}
# representing data with regression lines
DNA_x_diphos_graph <-ggplot(dataAll_EAA,aes(y=DNA,x=diphos))+
  geom_point(size=3, aes(shape=Type,fill=Type))+
  geom_smooth(method = "lm", se=T, color="black",linetype="dashed",
              formula = my.formula, size=1.5) +
  stat_poly_eq(formula = my.formula, label.x=0.9,size=6,
               aes(label = paste(..rr.label.., ..p.value.label..,sep = "~~~")), 
               parse = TRUE,
               rr.digits = 2,
               p.digits = 2) + 
  #ggtitle("Phosphodiesterase activity vs \nphosphodiesterase-labile organic P") +
  scale_linetype_manual("Soil compartment",values=c(1,8))+
  scale_fill_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_shape_manual("Soil compartment",values=c(21,23))+  
  ylab(expression(paste("Phosphodiesterase-labile P (",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  xlab(expression(paste("Phosphodiesterase activity (nmol MUB ",g^-1,h^-1,")",sep="")))+ labs(fill="")+
  Motheme_all2
DNA_x_diphos_graph
```

``` {r,  include=F}
svg("Fig9b.svg",width=8,height=8)
DNA_x_diphos_graph
dev.off()

png("Fig9b.png",width=12,height=12, units = "cm", res=800)
DNA_x_diphos_graph
dev.off()
```

\newpage
## Fig 8c: phosphomonoesterase-labile Porg vs alk. phosphomonoesterase

``` {r Fig 8c,fig.width = 4, fig.height=4,echo=F}  
#str(dataAll_EAA)

ggplot(dataAll_EAA,aes(y=GP,x=alk.phos, color=))+
  geom_point(aes(color=))+stat_smooth(method="lm")
```

```{r, include=F}
GP_x_alk.phos_model.red.lm<- lm(alk.phos~Type*CoverCrop*GP+Block, na.action=na.exclude, data=dataAll_EAA)
anova(GP_x_alk.phos_model.red.lm)

GP_x_alk.phos_model.red.lm_step<-step(GP_x_alk.phos_model.red.lm, direction="both")
anova(GP_x_alk.phos_model.red.lm_step)
summary(GP_x_alk.phos_model.red.lm_step)

anova(GP_x_alk.phos_model.red.lm_step, GP_x_alk.phos_model.red.lm)

# Step does not improve fit: GP_x_alk.phos_model.red.lm has lowest AIC
anova(GP_x_alk.phos_model.red.lm)
summary(GP_x_alk.phos_model.red.lm)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modGP_x_alk.phos<-GP_x_alk.phos_model.red.lm
plot(fitted(modGP_x_alk.phos),resid(modGP_x_alk.phos),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues
qqnorm(resid(modGP_x_alk.phos)) 
qqline(y=resid(modGP_x_alk.phos))
# histogram of residuals
hist(resid(modGP_x_alk.phos))

shapiro.test(resid(modGP_x_alk.phos))
```

##### Model fit 

```{r, echo=F}
formula(modGP_x_alk.phos)

anova(modGP_x_alk.phos) %>%
  kable(caption= "Relation between alk. phosphomonoesterase alktivity and phoD per g soil") %>%
  kable_styling()
```

representing data with regression lines
``` {r, echo=F}
GP_x_alk.phos_graph <-ggplot(dataAll_EAA,aes(y=GP,x=alk.phos))+
  geom_point(size=3, aes(shape=Type,fill=Type))+
  geom_smooth(method = "lm", se=T, color="black",linetype="dashed",
              formula = my.formula, size=1.5) +
  stat_poly_eq(formula = my.formula, label.x=0.9,size=6,
               aes(label = paste(..rr.label.., ..p.value.label..,sep = "~~~")), 
               parse = TRUE,
               rr.digits = 2,
               p.digits = 2) + 
  #ggtitle("Alkaline phosphomonoesterases alktivity vs \nphosphomonoesterases-labile organic P") +
  scale_linetype_manual("Soil compartment",values=c(1,8))+
  scale_fill_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_shape_manual("Soil compartment",values=c(21,23))+  
  ylab(expression(paste("Phosphomonoesterase-labile P (",mu,"g ", P[org]," ",g^-1,")",sep="")))+ labs(fill="")+
  xlab(expression(paste("Alk. phosphomonoesterase alktivity (nmol MUB ",g^-1,h^-1,")",sep="")))+ labs(fill="")+
  Motheme_all2
GP_x_alk.phos_graph
```

``` {r,  include=F}
svg("Fig9c.svg",width=8,height=8)
GP_x_alk.phos_graph
dev.off()

png("Fig9c.png",width=12,height=12, units = "cm", res=800)
GP_x_alk.phos_graph
dev.off()
```

\newpage
# Fig S4 P concentration in soybean grains

``` {r Fig S4,fig.width = 4, fig.height=4, echo=F}  
####### Fig S6 P_grain ########

# in contrast to cover crop biomass results, for soybean there are also the control results!
ggplot(dataAll , aes(x=CoverCrop, y=P_grain, fill=CoverCrop))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(g k",g^-1,")")))+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
```

``` {r, include=F}
summary(dataAll$P_grain)

P_grain_model.full.lm<- lm(log(P_grain)~CoverCrop+Block, na.action=na.omit,
                           data=dataAll)
anova(P_grain_model.full.lm )
# CoverCrop is not sig!
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modP_grain<-P_grain_model.full.lm
plot(fitted(modP_grain),resid(modP_grain),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modP_grain)) 
qqline(y=resid(modP_grain))
# histogram of residuals
hist(resid(modP_grain))

shapiro.test(resid(modP_grain))

modP_grain_graph <-as.data.frame(multcomp::cld(emmeans(modP_grain, ~CoverCrop, adjust="bon"),Letters=letters, type = "response"))
```

```{r, fig.width = 8, echo=F }

P_grain_graph <- ggplot(subset(modP_grain_graph), aes(x=CoverCrop, y=response, fill=CoverCrop, shape=CoverCrop))+
  scale_shape_manual(values=c(21,22, 23, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+
  #geom_text(aes(label=.group, x=CoverCrop, y=max(upper.CL)), position=position_dodge(0.8),vjust=-1, size=rel(6), color="black")+
  #ggtitle("P concentration in soybean grains") +
  ylab(expression(paste( "Soybean grain P concentration (g P k",g^-1,")")))+
  theme(legend.position="none")+xlab("")+
  ylim( min(modP_grain_graph$lower.CL),max(modP_grain_graph$upper.CL)+max(modP_grain_graph$upper.CL)/10)+
  geom_vline(xintercept=4.606, size=1.1, color="grey")+
  Motheme_all+
  scale_fill_manual("Cover Crop",values=c( "#31688EFF","#693376FF", "#FDE725FF", "#35B779FF"))
P_grain_graph
```

##### Model fit

```{r, echo=F}
formula(modP_grain)

anova(modP_grain) %>%
  kable(caption= "P concentration in soybean grains: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,  include=F}
svg("FigS4_Soybean_Pconc.svg",width=8,height=8)
P_grain_graph
dev.off()

png("FigS4_Soybean_Pconc.png",width=12,height=12, units = "cm", res=800)
P_grain_graph
dev.off()
```
\newpage
# Fig S5 Paq (Resin-P)

``` {r Fig S5,fig.width = 8, fig.height=4, echo=F}  
box_Paq <-ggplot(dataAll, aes(x=CoverCrop:Type, y=Paq, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(",mu, "g", P[mic], " ",g^-1,")")))+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = Paq, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ggtitle("Resin-P")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
box_Paq
```

``` {r, include=F}
modPaq.full<-
  lmer((Paq)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modPaq.full)
VarCorr(modPaq.full)
summary(modPaq.full)
vcov(modPaq.full)

modPaq.red<-
  lmer((Paq)~Type*CoverCrop*Month+
         Block+(1|TypeID), na.action=na.omit, data=dataAll)

# CoverCrop not sig! Type sig!

anova(modPaq.full, modPaq.red)

#red has lowest AIC
anova(modPaq.red)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modPaq<-modPaq.red
plot(fitted(modPaq),resid(modPaq),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modPaq)) 
qqline(y=resid(modPaq))
# histogram of residuals
hist(resid(modPaq))

shapiro.test(resid(modPaq))
#cld(emmeans(modPaq, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modPaq_cld <-as.data.frame(multcomp::cld(emmeans(modPaq, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))

modPaq_cld$.group<-stri_replace_all_fixed(modPaq_cld$.group, " ", "")

modPaq_cld$rhizo<-revalue(modPaq_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
Paq_graph <- ggplot(na.omit(modPaq_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)+1), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("")+
  ylab(expression(paste( "Available phosphate (",mu, "g ", P[resin], " ",g^-1,")")))+
  #theme(legend.position="none")+
  xlab("")+
  ylim( NA,max(modPaq_cld$upper.CL)+max(modPaq_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
Paq_graph

```


``` {r,  include=F}
svg("FigS5_Presin.svg",width=12, height=8)
Paq_graph
dev.off()

png("FigS5_Presin.png",width=32,height=16, units = "cm", res=800)
Paq_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modPaq.full)
anova(modPaq) %>%
  kable(caption= "Resin P: ANOVA output of fitted model") %>%
  kable_styling()
```

Resin-P only November to check if power increases by separating timepoints (does not)
``` {r,  include=F}

modPaq.November.T1<-
  lmer(Paq~Type*CoverCrop+
         Block+(1|Plot), na.action=na.omit, data=dataAll_T1)
anova(modPaq.November.T1)
multcomp::cld(emmeans(modPaq.November.T1, ~CoverCrop|Type, adjust="bon"),Letters=letters, type = "response")
# result is the same as with full model

modPaq.June.T3<-
  lm(Paq~Type*CoverCrop+
       Block, na.action=na.omit, data=dataAll_T3)
anova(modPaq.June.T3)
multcomp::cld(emmeans(modPaq.June.T3, ~CoverCrop|Type, adjust="bon"),Letters=letters, type = "response")
```

# (S6 is a Table with the p-values of the ANOVAs of a collection of the measured properties)

\newpage
# Fig S7 Prec (soil P sorption)

``` {r Fig S7,fig.width = 8, fig.height=4, echo=F}  
####### Fig S7 Prec (soil P sorption) ########

ggplot(dataAll, aes(x=CoverCrop:Type, y=Prec, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  ylab(expression(paste( "(% of added spike)")))+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = Prec, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ggtitle("P recovery")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modPrec.full<-
  lmer(sqrt(Prec)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modPrec.full)
VarCorr(modPrec.full)
summary(modPrec.full)
vcov(modPrec.full)

modPrec.red<-
  lmer(sqrt(Prec)~Type*CoverCrop*Month+
         Block+(1|MonthID), na.action=na.omit, data=dataAll)

# CoverCrop not sig! Type:Month sig!

anova(modPrec.full, modPrec.red)

#red has lowest AIC
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modPrec<-modPrec.red
plot(fitted(modPrec),resid(modPrec),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modPrec)) 
qqline(y=resid(modPrec))
# histogram of residuals
hist(resid(modPrec))

shapiro.test(resid(modPrec))

#cld(emmeans(modPrec, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modPrec_cld <-as.data.frame(multcomp::cld(emmeans(modPrec, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))

modPrec_cld$.group<-stri_replace_all_fixed(modPrec_cld$.group, " ", "")
modPrec_cld$rhizo<-revalue(modPrec_cld$Type, c("rhizo"="R", "bulk"="" ))
```
```{r, fig.width = 8, echo=F }
Prec_graph <- ggplot(na.omit(modPrec_cld), aes(x=CoverCrop:Type, y=response*100, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL*100, ymax=upper.CL*100), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response*100))+
  ##geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)*100), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("")+
  ylab(expression(paste( "recovery of added P-spike (%)")))+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(modPrec_cld$upper.CL)+max(modPrec_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
Prec_graph

```


``` {r,  include=F}
svg("FigS7_Prec.svg",width=12, height=8)
Prec_graph
dev.off()

png("FigS7_Prec.png",width=32,height=16, units = "cm", res=800)
Prec_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modPrec)
anova(modPrec) %>%
  kable(caption= "Prec: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
# Fig S8 Bacterial 16S rRNA

``` {r Fig S8,fig.width = 4, fig.height=4, echo=F}  
ggplot(dataAll , aes(x=CoverCrop:Type, y=bact_16S_dw, fill=CoverCrop:Type))+
  geom_boxplot()+
  geom_text(data = dataAll, aes(x = CoverCrop:Type, group=CoverCrop:Type, y = bact_16S_dw, label = Plot), position=position_dodge(1), hjust=.5)+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
bact_16S_dw_model.full<- lmer(log(bact_16S_dw)~CoverCrop*Type+Block+(1|Plot), na.action=na.exclude, data=dataAll)
anova(bact_16S_dw_model.full)
VarCorr(bact_16S_dw_model.full)

bact_16S_dw_model.full.lm<- lm(log(bact_16S_dw)~CoverCrop*Type+Block, na.action=na.exclude, data=dataAll)
multcomp::cld(emmeans(bact_16S_dw_model.full.lm, ~Type, adjust="bon"),Letters=letters, type = "response")
anova(bact_16S_dw_model.full.lm)

anova(bact_16S_dw_model.full, bact_16S_dw_model.full.lm)
# lm has lower AIC.
anova(bact_16S_dw_model.full.lm)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modbact_16S_dw<-bact_16S_dw_model.full.lm
plot(fitted(modbact_16S_dw),resid(modbact_16S_dw),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modbact_16S_dw)) 
qqline(y=resid(modbact_16S_dw))
# histogram of residuals
hist(resid(modbact_16S_dw))

shapiro.test(resid(modbact_16S_dw))

modbact_16S_dw_graph <-as.data.frame(multcomp::cld(emmeans(modbact_16S_dw, ~CoverCrop|Type, adjust="bon"),Letters=letters, type = "response"))

modbact_16S_dw_graph$rhizo<-revalue(modbact_16S_dw_graph$Type, c("rhizo"="R", "bulk"="" ))
#modbact_16S_dw_graph$rhizo[8]<-""
```

```{r, fig.width = 8, echo=F }

bact_16S_dw_graph <- ggplot(subset(modbact_16S_dw_graph), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22,22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  ##geom_text(aes(label=.group, x=CoverCrop:Type, y=max(modbact_16S_dw_graph$upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle(" ") +
  ylab(expression(paste("16S rRNA (copies ", g^-1, " soil)",sep="")))+ labs(fill="")+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( 60, 450)+#min(modbact_16S_dw_graph$lower.CL)*0.9,max(modbact_16S_dw_graph$upper.CL)+max(modbact_16S_dw_graph$upper.CL)/5)+
  Motheme_all+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
bact_16S_dw_graph
```

##### Model fit 

```{r, echo=F}
formula(modbact_16S_dw)

anova(modbact_16S_dw) %>%
  kable(caption= "Bacterial 16S copies per g soil dw: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,  include=F}
svg("FigS8_bact_16S_dw.svg",width=8,height=8)
bact_16S_dw_graph
dev.off()

png("FigS8_bact_16S_dw.png",width=32,height=16, units = "cm", res=800)
bact_16S_dw_graph
dev.off()
```

\newpage
# Fig S9 phoD (per g soil)

``` {r Fig S9,fig.width = 4, fig.height=4, echo=F}  
ggplot(dataAll , aes(x=CoverCrop:Type, y=phoD_dw, fill=CoverCrop:Type))+
  geom_boxplot()+
  geom_text(data = dataAll, aes(x = CoverCrop:Type, group=CoverCrop:Type, y = phoD_dw, label = Plot), position=position_dodge(1), hjust=.5)+
  Motheme_all+xlab("")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
phoD_dw_model.full<- lmer(log(phoD_dw)~CoverCrop*Type+Block+(1|Plot), na.action=na.exclude, data=dataAll)
anova(phoD_dw_model.full)
VarCorr(phoD_dw_model.full)

phoD_dw_model.full.lm<- lm(log(phoD_dw)~CoverCrop*Type+Block, na.action=na.exclude, data=dataAll)
anova(phoD_dw_model.full.lm)

phoD_dw_model.red.lm<- lm(log(phoD_dw)~Type+Block, na.action=na.exclude, data=dataAll)
anova(phoD_dw_model.red.lm)
anova(phoD_dw_model.full, phoD_dw_model.full.lm, phoD_dw_model.red.lm)

# .red.lm has lowest AIC, but 
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modphoD_dw<-phoD_dw_model.full.lm
plot(fitted(modphoD_dw),resid(modphoD_dw),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modphoD_dw)) 
qqline(y=resid(modphoD_dw))
# histogram of residuals
hist(resid(modphoD_dw))

shapiro.test(resid(modphoD_dw))

modphoD_dw_graph <-as.data.frame(multcomp::cld(emmeans(phoD_dw_model.full, ~CoverCrop:Type, adjust="bon"),Letters=letters, type = "response"))
modphoD_dw_graph$rhizo<-revalue(modphoD_dw_graph$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 6, echo=F }

phoD_dw_graph <- ggplot(subset(modphoD_dw_graph), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22,22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  ##geom_text(aes(label=.group, x=CoverCrop:Type, y=max(modphoD_dw_graph$upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ggtitle("") +
  ylab(expression(paste("phoD (copies ", g^-1,")",sep="")))+ labs(fill="")+
  #theme(legend.position="none")+
  xlab("")+
  Motheme_all+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
phoD_dw_graph
```

##### Model fit 

```{r, echo=F}
formula(modphoD_dw)

anova(modphoD_dw) %>%
  kable(caption= "Bacterial phoD copies per g soil dw: ANOVA output of fitted model") %>%
  kable_styling()
```

``` {r,  include=F}
svg("FigS9_phoD_dw.svg",width=8,height=8)
phoD_dw_graph
dev.off()

png("FigS9_phoD_dw.png",width=32,height=16, units = "cm", res=800)
phoD_dw_graph
dev.off()
```
# (S10 is a Table presenting the p-values of the ANOVAs of the preceeding figures)

\newpage
# Fig S11 specific potential enzyme activities
## Fig S11a ac.phos on Pmic basis 

``` {r Fig S11a,fig.width = 8, fig.height=4, echo=F}  
dataAll$ac.phos.per.Pmic<-dataAll$ac.phos/dataAll$Pmic

ggplot(dataAll, aes(x=CoverCrop:Type, y=ac.phos.per.Pmic, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = ac.phos.per.Pmic, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ylab(expression(paste("Specific acid phosphomonoesterase activity (",nmol, " MUB ","per ",mu, "g ",P[mic]," ",h^-1,")")))+  ggtitle("Alkaline phosphomonoesterase activity under\n cover crops")+
  ggtitle("")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modac.phos.per.Pmic.full<-
  lmer(log(ac.phos.per.Pmic)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modac.phos.per.Pmic.full)
VarCorr(modac.phos.per.Pmic.full)
summary(modac.phos.per.Pmic.full)
vcov(modac.phos.per.Pmic.full)

# no sig!

modac.phos.per.Pmic.red<-get_model(step(modac.phos.per.Pmic.full, direction="both",  reduce.random=F,keep=c("Block")))
anova(modac.phos.per.Pmic.full, modac.phos.per.Pmic.red)

# no sig! using full model to estimate means for figure

modac.phos.per.Pmic.red_cld <-as.data.frame(multcomp::cld(emmeans(modac.phos.per.Pmic.full, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))

modac.phos.per.Pmic.red_cld$rhizo<-revalue(modac.phos.per.Pmic.red_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
ac.phos.per.Pmic_graph <- ggplot(na.omit(modac.phos.per.Pmic.red_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ylab(expression(paste("Specific acid phosphomonoesterase activity (",nmol, " MUB ","per ",mu, "g ",P[mic]," ",h^-1,")")))+  ggtitle("Alkaline phosphomonoesterase activity under\n cover crops")+
  ggtitle(" ")+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(modac.phos.per.Pmic.red_cld$upper.CL)+max(modac.phos.per.Pmic.red_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
ac.phos.per.Pmic_graph
```

``` {r,  include=F}
svg("FigS11a.svg",width=12, height=8)
ac.phos.per.Pmic_graph
dev.off()

png("FigS11a.png",width=32,height=16, units = "cm", res=800)
ac.phos.per.Pmic_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modac.phos.per.Pmic.red)
anova(modac.phos.per.Pmic.red) %>%
  kable(caption= "Acid phosphomonoesterase activity by Pmic: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
## Fig S11b alk.phos on Pmic basis

``` {r Fig S11b,fig.width = 8, fig.height=4,echo=F}  
dataAll$alk.phos.per.Pmic<-dataAll$alk.phos/dataAll$Pmic
ggplot(dataAll, aes(x=CoverCrop:Type, y=alk.phos.per.Pmic, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = alk.phos.per.Pmic, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ylab(expression(paste("(",nmol, " MUB ","per ",mu, "g ",P[mic]," ",h^-1,")")))+  ggtitle("Alkaline phosphomonoesterase activity under\n cover crops")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
modalk.phos.per.Pmic.full<-
  lmer(log(alk.phos.per.Pmic)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(modalk.phos.per.Pmic.full)
VarCorr(modalk.phos.per.Pmic.full)
summary(modalk.phos.per.Pmic.full)
vcov(modalk.phos.per.Pmic.full)

# Type sig!
modalk.phos.per.Pmic.red<-get_model(step(modalk.phos.per.Pmic.full, direction="both",  reduce.random=F,keep=c("Block")))
anova(modalk.phos.per.Pmic.red)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modalk.phos.per.Pmic<-modalk.phos.per.Pmic.full
plot(fitted(modalk.phos.per.Pmic),resid(modalk.phos.per.Pmic),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modalk.phos.per.Pmic)) 
qqline(y=resid(modalk.phos.per.Pmic))
# histogram of residuals
hist(resid(modalk.phos.per.Pmic.red))

shapiro.test(resid(modalk.phos.per.Pmic))
#cld(emmeans(modalk.phos.per.Pmic, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
modalk.phos.per.Pmic_cld <-as.data.frame(multcomp::cld(emmeans(modalk.phos.per.Pmic, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))

modalk.phos.per.Pmic_cld$rhizo<-revalue(modalk.phos.per.Pmic_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
alk.phos.per.Pmic_graph <- ggplot(na.omit(modalk.phos.per.Pmic_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ylab(expression(paste("Specific alk. phosphomonoesterase activity (",nmol, " MUB ","per ",mu, "g ",P[mic]," ",h^-1,")")))+  ggtitle("Alkaline phosphomonoesterase activity under\n cover crops")+
  ggtitle(" ")+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(modalk.phos.per.Pmic.red_cld$upper.CL)+max(modalk.phos.per.Pmic.red_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
alk.phos.per.Pmic_graph
```


``` {r,  include=F}
svg("FigS11b.svg",width=12, height=8)
alk.phos.per.Pmic_graph
dev.off()

png("FigS11b.png",width=32,height=16, units = "cm", res=800)
alk.phos.per.Pmic_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(modalk.phos.per.Pmic)
anova(modalk.phos.per.Pmic) %>%
  kable(caption= "Alkaline phosphomonoesterase activity by Pmic: ANOVA output of fitted model") %>%
  kable_styling()
```


\newpage
## Fig S11c diphos on Pmic basis

``` {r Fig S11c,fig.width = 8, fig.height=4,echo=F}  
dataAll$diphos.per.Pmic<-dataAll$diphos/dataAll$Pmic

ggplot(dataAll, aes(x=CoverCrop:Type, y=diphos.per.Pmic, fill=CoverCrop:Type))+
  geom_boxplot()+Motheme_all+xlab("")+
  geom_text(data = dataAll, aes(x = CoverCrop:Type,  y = diphos.per.Pmic, group=CoverCrop:Type,label = Plot), 
            position=position_dodge(1), hjust=.5)+
  ylab(expression(paste("(",nmol, " MUB ","per ",mu, "g ",P[mic]," ",h^-1,")")))+  ggtitle("Alkaline phosphomonoesterase activity under\n cover crops")+
  ggtitle("Specific phosphodiesterase activity")+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
```

``` {r, include=F}
moddiphos.per.Pmic.full<-
  lmer(log(diphos.per.Pmic)~Type*CoverCrop*Month+
         Block+(1|TypeID)+(1|MonthID), na.action=na.omit, data=dataAll)
anova(moddiphos.per.Pmic.full)
VarCorr(moddiphos.per.Pmic.full)
summary(moddiphos.per.Pmic.full)
vcov(moddiphos.per.Pmic.full)

# Type:Month sig!
moddiphos.per.Pmic.red<-get_model(step(moddiphos.per.Pmic.full, direction="both",  reduce.random=F,keep=c("Block")))

anova(moddiphos.per.Pmic.red)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
moddiphos.per.Pmic<-moddiphos.per.Pmic.full
plot(fitted(moddiphos.per.Pmic),resid(moddiphos.per.Pmic),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(moddiphos.per.Pmic)) 
qqline(y=resid(moddiphos.per.Pmic))
# histogram of residuals
hist(resid(moddiphos.per.Pmic))

shapiro.test(resid(moddiphos.per.Pmic))
#cld(emmeans(moddiphos.per.Pmic, ~Type|Month, adjust="bon"),Letters=letters, type = "response")
moddiphos.per.Pmic_cld <-as.data.frame(multcomp::cld(emmeans(moddiphos.per.Pmic, ~CoverCrop:Type|Month, adjust="bon"),Letters=letters, type = "response"))

moddiphos.per.Pmic_cld$rhizo<-revalue(moddiphos.per.Pmic_cld$Type, c("rhizo"="R", "bulk"="" ))
```

```{r, fig.width = 8, echo=F }
diphos.per.Pmic_graph <- ggplot(na.omit(moddiphos.per.Pmic_cld), aes(x=CoverCrop:Type, y=response, fill=CoverCrop:Type, shape=CoverCrop:Type))+
  scale_shape_manual("Cover crop:soil compartment",values=c(21,21,22, 22, 23,23, 24, 24))+
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1)+
  geom_point(size=6,color="black")+   geom_text(aes(label=rhizo, x=CoverCrop:Type, y=response))+
  #geom_text(aes(label=.group, x=CoverCrop:Type, y=max(upper.CL)), position=position_dodge(0.8),vjust=-0.4, size=rel(6), color="black")+
  ylab(expression(paste("Specific phosphodiesterase activity (",nmol, " MUB ","per ",mu, "g ",P[mic]," ",h^-1,")")))+  ggtitle("Alkaline phosphomonoesterase activity under\n cover crops")+
  ggtitle(" ")+
  #theme(legend.position="none")+
  xlab("")+
  #ylim( NA,max(moddiphos.per.Pmic_cld$upper.CL)+max(moddiphos.per.Pmic_cld$upper.CL)/10)+
  theme_bw()+
  Motheme_all+
  facet_grid(~Month,scales = "free_x")+
  scale_fill_manual("Cover crop:soil compartment",values=c("#000000", "#FFFFFF" ,"#000000", "#FFFFFF" , "#000000", "#FFFFFF" , "#000000","#FFFFFF" ))
diphos.per.Pmic_graph
```

``` {r,  include=F}
svg("FigS11c.svg",width=12, height=8)
diphos.per.Pmic_graph
dev.off()

png("FigS11c.png",width=32,height=16, units = "cm", res=800)
diphos.per.Pmic_graph
dev.off()
```

##### Model fit

```{r, echo=F}
formula(moddiphos.per.Pmic.full)
anova(moddiphos.per.Pmic) %>%
  kable(caption= "Phosphodiesterase activity by Pmic: ANOVA output of fitted model") %>%
  kable_styling()
```

\newpage
# Fig S12: Fungi vs N-acetyl-hexosaminidase

``` {r Fig S12,fig.width = 4, fig.height=4, echo=F}  
#str(dataAll)

ggplot(dataAll,aes(y=NAG,x=log(Fungi), color=Type))+
  geom_point(aes(color=Type))+stat_smooth(method="lm")
```

```{r, include=F}
Fungi_x_NAG_model.full<- lmer(NAG~CoverCrop*Date*Type*log(Fungi)+Block+(1|Plot)+(1|MonthID)+(1|TypeID), na.action=na.exclude, data=dataAll)
anova(Fungi_x_NAG_model.full)
VarCorr(Fungi_x_NAG_model.full)

Fungi_x_NAG_model.red.lm<- lm(NAG~Type*log(Fungi)+Block, na.action=na.exclude, data=dataAll)
anova(Fungi_x_NAG_model.red.lm)

Fungi_x_NAG_model.red1<- lmer(NAG~CoverCrop*Date*Type*log(Fungi)+Block+(1|MonthID), na.action=na.exclude, data=dataAll)
anova(Fungi_x_NAG_model.red1)
VarCorr(Fungi_x_NAG_model.red1)

Fungi_x_NAG_model.red.lm_step<-get_model(step(Fungi_x_NAG_model.red1, direction="both"))
anova(Fungi_x_NAG_model.red.lm_step)
summary(Fungi_x_NAG_model.red.lm_step)

anova(Fungi_x_NAG_model.full, Fungi_x_NAG_model.red.lm, Fungi_x_NAG_model.red.lm_step)

# Fungi_x_NAG_model.red.lm_step has lowest AIC
anova(Fungi_x_NAG_model.red.lm_step)
```

``` {r ,fig.width = 4, fig.height=4, echo=F} 
modFungi_x_NAG<-Fungi_x_NAG_model.red.lm_step
#modFungi<-Fungi_model.full
plot(fitted(modFungi_x_NAG),resid(modFungi_x_NAG),xlab='Fitted Values', ylab='Residuals'); abline(h=0)
# create Q-Q Plot of the residues 
qqnorm(resid(modFungi_x_NAG)) 
qqline(y=resid(modFungi_x_NAG))
# histogram of residuals
hist(resid(modFungi_x_NAG))

shapiro.test(resid(modFungi_x_NAG))
```

##### Model fit 

```{r, echo=F}
formula(modFungi_x_NAG)

anova(modFungi_x_NAG) %>%
  kable(caption= "Relation between fungal abundance and NAG activity") %>%
  kable_styling()
```

``` {r,  echo=F}
# representing data with regression lines
Fungi_x_NAG_graph_Type <-ggplot(dataAll_T1,aes(y=diphos,x=log(Fungi), color=Type, shape=Type, fill=Type))+
  #  ggplot(dataAll,aes(y=NAG,x=log(Fungi), color=Type))+
  geom_point(aes(color=Type, shape=Type),size=2.5)+stat_smooth(method="lm")+
  ggtitle("") +
  scale_color_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_fill_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_shape_manual("Soil compartment",values=c(21,23))+
  xlab(expression(paste("Fungi (log[ng PLFA",  g^-1,"soil])",sep="")))+ labs(fill="")+
  ylab(expression(paste("N-acetyl-hexosaminidase activity (nmol MUB ",g^-1,h^-1,")",sep="")))+ labs(fill="")+
  Motheme_all2
Fungi_x_NAG_graph_Type

Fungi_x_NAG_graph <-ggplot(dataAll,aes(y=NAG,x=Fungi, shape=Type))+
  geom_point(size=3, aes(fill=Type, shape=Type))+
  geom_smooth(method = "lm", se=T, aes(color=Type),linetype="dashed",
              formula = my.formula, size=1.5) +
  stat_poly_eq(formula = my.formula, label.x=0.9,size=6,
               aes(label = paste(..rr.label.., ..p.value.label..,sep = "~~~")), 
               parse = TRUE,
               rr.digits = 2,
               p.digits = 2) + 
  scale_color_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_fill_manual("Soil compartment",values=c("#a6611a","#018571"))+
  scale_shape_manual("Soil compartment",values=c(21,23))+
  ggtitle("Fungal vs N-acetyl-hexosaminidase  activity") +
  xlab(expression(paste("Fungi (log[ng PLFA",  g^-1,"soil])",sep="")))+ labs(fill="")+
  ylab(expression(paste("N-acetyl-hexosaminidase activity (nmol MUB ",g^-1,h^-1,")",sep="")))+ labs(fill="")+
  Motheme_all2
Fungi_x_NAG_graph
```

``` {r,  include=F}
svg("FigS12.svg",width=8,height=8)
Fungi_x_NAG_graph
dev.off()

png("FigS12.png",width=12,height=12, units = "cm", res=800)
Fungi_x_NAG_graph#+  theme(legend.position="none")
dev.off()
```
