---
title: "Plant-mycorrhizal niche mapping analysis and figures"
output: html_notebook
---
Created by: Michael Van Nuland & Kabir Peay (October 2019 - February 2020)
  
Analysis and data visualization associated with "Symbiotic niche mapping reveals nutrient specialization and functional complementarity among ectomycorrhizal fungi".
  
# Load packages, import and manage data 
```{r}

library(lme4)
library(lmerTest)
library(geometry)
library(plotly)
library(fields)
library(reshape2)
library(ggplot2)
library(viridis)
library(RColorBrewer)

pinus.myc.niche.dat <- read.csv(file="~/Desktop/Pinus_Myc_Niche_data.csv")

# Re-defining mycorrhizal treatments based on Thelephora contamination issues found after we measured ECM root colonization

control <- subset(pinus.myc.niche.dat, pinus.myc.niche.dat$MYC=="Control" & 
                    pinus.myc.niche.dat$Thelephora_colonization < 2 & 
                    pinus.myc.niche.dat$Suillus_colonization < 2)

thelephora <- subset(pinus.myc.niche.dat, pinus.myc.niche.dat$MYC=="Thelephora" & pinus.myc.niche.dat$Thelephora_col > 10)

mixed <- subset(pinus.myc.niche.dat, pinus.myc.niche.dat$Suillus_col > 10 & pinus.myc.niche.dat$Thelephora_col > 10) 

pinus.myc.niche.dat.redef <- as.data.frame(rbind(control, thelephora, mixed)) #MYC_recode treatment column reflects new designations
pinus.myc.niche.dat.redef <- na.exclude(pinus.myc.niche.dat.redef)

```
  
  
  
# FIGURE 1

Total plant biomass responses across N and P fertilization gradients and convex hull niche dimensions.

```{r}

TotalBiomass.mod <- lmer(log(TotalBiomass) ~ MYC_recode * log(Nitrogen) * log(Phosphorus) + (1|REP), data = pinus.myc.niche.dat.redef)
anova(TotalBiomass.mod)

# TABLE 1
                                          #Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
#MYC_recode                                0.0541  0.0270     2 221.02   0.2826   0.75410    
#log(Nitrogen)                            25.5485 25.5485     1 221.03 267.0666 < 2.2e-16 ***
#log(Phosphorus)                           0.1439  0.1439     1 221.01   1.5045   0.22128    
#MYC_recode:log(Nitrogen)                  0.5364  0.2682     2 221.02   2.8034   0.06276 .  
#MYC_recode:log(Phosphorus)                0.1236  0.0618     2 221.02   0.6458   0.52521    
#log(Nitrogen):log(Phosphorus)             4.8656  4.8656     1 221.39  50.8615 1.391e-11 ***
#MYC_recode:log(Nitrogen):log(Phosphorus)  0.0524  0.0262     2 221.16   0.2737   0.76082


############

TotalBiomass.mod.N.resids <- lm(log(TotalBiomass) ~ log(Nitrogen), data = pinus.myc.niche.dat.redef)
  N.resids <- resid(TotalBiomass.mod.N.resids)

TotalBiomass.mod.P.resids <- lm(log(TotalBiomass) ~ log(Phosphorus), data = pinus.myc.niche.dat.redef)
  P.resids <- resid(TotalBiomass.mod.P.resids)

cbind(pinus.myc.niche.dat.redef, N.resids, P.resids)

Myc_trts <- factor(pinus.myc.niche.dat.redef$MYC_recode, levels = c("Control", "Thelephora", "Suillus_Thelephora"))


###  Plotting linear regressions

ggplot(aes(x=log(Nitrogen), y=P.resids), data = pinus.myc.niche.dat.redef) +
  geom_point(position = position_jitterdodge(dodge.width = 0.3, jitter.width = 0.1), aes(x=log(Nitrogen), y=P.resids, color = Myc_trts), size = 1, data = pinus.myc.niche.dat.redef, show.legend = T) +
  geom_smooth(aes(x=log(Nitrogen), y=P.resids, color = Myc_trts), method = "lm", se=F, lwd=0.75, data=pinus.myc.niche.dat.redef, show.legend = F) +
  scale_color_manual(values = c("#bab0ac","#4e79a7", "#f28e2b")) +
  xlab("Log Nitrogen (mg/kg soil)") +
  ylab("Total Plant Biomass (residuals)") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())


ggplot(aes(x=log(Phosphorus), y=N.resids), data = pinus.myc.niche.dat.redef) +
  geom_point(position = position_jitterdodge(dodge.width = 0.3, jitter.width = 0.1), aes(x=log(Phosphorus), y=N.resids, color = Myc_trts), size = 1, data=pinus.myc.niche.dat.redef, show.legend = T) +
  geom_smooth(aes(x=log(Phosphorus), y=N.resids, color = Myc_trts), method ="lm", se=F, lwd=0.75, data=pinus.myc.niche.dat.redef, show.legend = F) +
  scale_color_manual(values = c("#bab0ac","#4e79a7", "#f28e2b")) +
  xlab("Log Phosphorus (mg/kg soil)") +
  ylab("Total Plant Biomass (residuals)") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())



###   Constructing plant biomass convex hulls and plotting response surfaces

# Create & smooth 7x7 matrices of total biomass responses across N and P gradients

control.mat <- acast(control, Nitrogen~Phosphorus, value.var = "TotalBiomass", mean, drop=TRUE)
thel.mat <- acast(thelephora,Nitrogen~Phosphorus,value.var = "TotalBiomass",mean, drop=TRUE)
mixed.mat <- acast(mixed, Nitrogen~Phosphorus, value.var = "TotalBiomass", mean, drop=TRUE)

control.mat.smooth <- image.smooth(control.mat, theta=1) #smooths over missing matrix values
thel.mat.smooth<-image.smooth(thel.mat, theta=1)
mixed.mat.smooth<-image.smooth(mixed.mat, theta=1)

# Create convex hulls and calculate volumes (single measures)

control.hull <- convhulln(melt(control.mat.smooth$z), output.options=TRUE)
  control.hull$vol  #7315.871
  
thel.hull <- convhulln(melt(thel.mat.smooth$z), output.options=TRUE)
  thel.hull$vol     #7490.693

mixed.hull <- convhulln(melt(mixed.mat.smooth$z), output.options=TRUE)
  mixed.hull$vol    #6448.734

# Response surface plotting

axz <- list(title="Total biomass (g)", range = c(185, 1100)) #Setting Z-axes to same limits across plots

# Control
plot_ly( z = ~control.mat.smooth$z) %>% add_surface(opacity = 0.95) %>% 
  layout(
  title="Control",
  scene=list(
    aspectmode = 'cube',
    xaxis=list(title="Phosphorus", autorange = "reversed"),
    yaxis=list(title="Nitrogen", autorange = "reversed"), 
    zaxis=axz))

# Thelephora
plot_ly( z = ~thel.mat.smooth$z) %>% add_surface(opacity = 0.95) %>% layout(
  title="Thelephora",
  scene=list(
    aspectmode = 'cube',
    xaxis=list(title="Phosphorus", autorange = "reversed"),
    yaxis=list(title="Nitrogen", autorange = "reversed"), 
    zaxis=axz))

# Mixed (Thelephora + Suillus)
plot_ly( z = ~mixed.mat.smooth$z) %>% add_surface(opacity = 0.95) %>% layout(
  title="Thelephora + Suillus",
  scene=list(
    aspectmode = 'cube',
    xaxis=list(title="Phosphorus", autorange = "reversed"),
    yaxis=list(title="Nitrogen", autorange = "reversed"), 
    zaxis=axz))



###  Random subsampling and bootstrapping approach to create distributions and calculate error estimates for convex hull volumes from each mycorrhizal treatment

## CONTROL

set.seed(1324)

control.mat.mean.hull.vols <- NA

for(i in 1:1000){
  
  control_matrix_list <- list()
  volume_output_control <- NA
  
  for(j in 1:100){
    original.matrix<-control
    random.number<-rnorm(length(original.matrix$Nitrogen),0,1)
    sample.order<-order(random.number)
    resampled.matrix<-original.matrix[sample.order[1:round(length(sample.order)*0.5)],]
    sampled.matrix <- acast(resampled.matrix, Nitrogen~Phosphorus, value.var = "TotalBiomass", mean, drop=TRUE)
    sampled.mat.smooth<- image.smooth(sampled.matrix, theta=1) #smooths over missing matrix values
    sampled.hull <- convhulln(melt(sampled.mat.smooth$z), output.options=TRUE)
    volume_output_control[j]<-sampled.hull$vol
    control_matrix_list[[j]] <- sampled.mat.smooth$z
  }
  control.ids.to.remove <- sapply(control_matrix_list, function(i) length(i) < 49) 
  control_matrix_list <- control_matrix_list[!control.ids.to.remove] 
  control.mat.mean <- Reduce("+", control_matrix_list) / length(control_matrix_list) # Calculate means by x/y elements
  control.mat.mean.hull <- convhulln(melt(control.mat.mean), output.options=TRUE)
  control.mat.mean.hull.vols[i] <- control.mat.mean.hull$vol
}

mean(control.mat.mean.hull.vols, na.rm=T)
  sd(control.mat.mean.hull.vols, na.rm=T)


## THELEPHORA

set.seed(1324)
  
thelephora.mat.mean.hull.vols <- NA

for(i in 1:1000){
  
  thelephora_matrix_list <- list()
  volume_output_thelephora <- NA
  
  for(j in 1:100){
    original.matrix<-thelephora
    random.number<-rnorm(length(original.matrix$Nitrogen),0,1)
    sample.order<-order(random.number)
    resampled.matrix<-original.matrix[sample.order[1:round(length(sample.order)*0.5)],]
    sampled.matrix <- acast(resampled.matrix, Nitrogen~Phosphorus, value.var = "TotalBiomass", mean, drop=TRUE)
    sampled.mat.smooth<- image.smooth(sampled.matrix, theta=1) #smooths over missing matrix values
    sampled.hull <- convhulln(melt(sampled.mat.smooth$z), output.options=TRUE)
    volume_output_thelephora[j]<-sampled.hull$vol
    thelephora_matrix_list[[j]] <- sampled.mat.smooth$z
  }
  thelephora.ids.to.remove <- sapply(thelephora_matrix_list, function(i) length(i) < 49) 
  thelephora_matrix_list <- thelephora_matrix_list[!thelephora.ids.to.remove]
  thelephora.mat.mean <- Reduce("+", thelephora_matrix_list) / length(thelephora_matrix_list) # Calculate means by x/y elements
  thelephora.mat.mean.hull <- convhulln(melt(thelephora.mat.mean), output.options=TRUE)
  thelephora.mat.mean.hull.vols[i] <- thelephora.mat.mean.hull$vol
}

mean(thelephora.mat.mean.hull.vols, na.rm=T)
  sd(thelephora.mat.mean.hull.vols, na.rm=T)


## MIXED (Thelephora + Suillus)

set.seed(1324)  
  
mixed.mat.mean.hull.vols <- NA

for(i in 1:1000){
  
  mixed_matrix_list <- list()
  volume_output_mixed <- NA

  for(j in 1:100){
    original.matrix<-mixed
    random.number<-rnorm(length(original.matrix$Nitrogen),0,1)
    sample.order<-order(random.number)
    resampled.matrix<-original.matrix[sample.order[1:round(length(sample.order)*0.5)],]
    sampled.matrix <- acast(resampled.matrix, Nitrogen~Phosphorus, value.var = "TotalBiomass", mean, drop=TRUE)
    sampled.mat.smooth<- image.smooth(sampled.matrix, theta=1) #smooths over missing matrix values
    sampled.hull <- convhulln(melt(sampled.mat.smooth$z), output.options=TRUE)
    volume_output_mixed[j]<-sampled.hull$vol
    mixed_matrix_list[[j]] <- sampled.mat.smooth$z
  }
  mixed.ids.to.remove <- sapply(mixed_matrix_list, function(i) length(i) < 49) 
  mixed_matrix_list <- mixed_matrix_list[!mixed.ids.to.remove]
  mixed.mat.mean <- Reduce("+", mixed_matrix_list) / length(mixed_matrix_list) # Calculate means by x/y elements
  mixed.mat.mean.hull <- convhulln(melt(mixed.mat.mean), output.options=TRUE)
  mixed.mat.mean.hull.vols[i] <- mixed.mat.mean.hull$vol
}

mean(mixed.mat.mean.hull.vols, na.rm=T)
  sd(mixed.mat.mean.hull.vols, na.rm=T)

##

t.test(thelephora.mat.mean.hull.vols, control.mat.mean.hull.vols, conf.level=0.95)
t.test(mixed.mat.mean.hull.vols, control.mat.mean.hull.vols, conf.level=0.95)
t.test(mixed.mat.mean.hull.vols, thelephora.mat.mean.hull.vols, conf.level=0.95)

##
  
# Combine the density distributions for plotting
vol.control.thel.mix.dat <- data.frame(dens = c(control.mat.mean.hull.vols, thelephora.mat.mean.hull.vols, mixed.mat.mean.hull.vols), lines = rep(c("control", "thel", "mix"), each = 1000))
  
# Control
ggplot() + 
  geom_histogram(data = subset(vol.control.thel.mix.dat, lines=="control"), aes(x = dens), alpha = 1, color = "#bab0ac", fill="#bab0ac", binwidth=50) +
  geom_vline(xintercept = mean(control.mat.mean.hull.vols), linetype="dashed") +
  ylab("Count") +
  xlab("Convex Hull Volume") +
  ylim(0,175) +
  xlim(5500, 7750) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          strip.background = element_blank(), strip.text = element_blank())

# Thelephora
ggplot() + 
  geom_histogram(data = subset(vol.control.thel.mix.dat, lines=="thel"), aes(x = dens), alpha = 1, color = "#4e79a7", fill="#4e79a7", binwidth=50) +
  geom_vline(xintercept = mean(thelephora.mat.mean.hull.vols), linetype="dashed") +
  ylab("Count") +
  xlab("Convex Hull Volume") +
  ylim(0,175) +
  xlim(5500, 7750) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          strip.background = element_blank(), strip.text = element_blank())

# Mixed (Thelephora + Suillus)
ggplot() + 
  geom_histogram(data = subset(vol.control.thel.mix.dat, lines=="mix"), aes(x = dens), alpha = 1, color = "#f28e2b", fill="#f28e2b", binwidth=50) +
  geom_vline(xintercept = mean(mixed.mat.mean.hull.vols), linetype="dashed") +
  ylab("Count") +
  xlab("Convex Hull Volume") +
  ylim(0,175) +
  xlim(5500, 7750) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          strip.background = element_blank(), strip.text = element_blank())

```



# FIGURE 2

Differences in plant growth responses between mycorrhizal treatments, expressed as % of Control.

```{r}

###  Calculating growth differences

change.mat.thel <- 100*((thel.mat.smooth$z-control.mat.smooth$z)/control.mat.smooth$z)
change.mat.mixed <- 100*((mixed.mat.smooth$z-control.mat.smooth$z)/control.mat.smooth$z)
change.mat.myco <- 100*((mixed.mat.smooth$z-thel.mat.smooth$z)/thel.mat.smooth$z)

axz2 <- list(title="Biomass difference (% of Control)", range = c(-37, 70)) #Setting Z-axes to same limits across plots


# Thelephora vs. Control

0 - min(change.mat.thel) / (max(change.mat.thel) - min(change.mat.thel)) # calculate midpoint value

plot_ly( z = ~change.mat.thel) %>% add_surface(opacity = 1, colorscale = list(
  c(0, 0.2336282, 1), c("blue", "white", "red"))) %>% layout(
    title="Thelephora - Control",
    scene=list(
      aspectmode = 'cube',
      xaxis=list(title="Phosphorus", autorange = "reversed"),
      yaxis=list(title="Nitrogen", autorange = "reversed"), 
      zaxis=axz2))


# Mixed (Thelephora+Suillus) vs. Control

0 - min(change.mat.mixed) / (max(change.mat.mixed) - min(change.mat.mixed)) # calculate midpoint value

plot_ly( z = ~change.mat.mixed) %>% add_surface(opacity = 1, colorscale = list(
  c(0, 0.4166209, 1), c("blue","white", "red"))) %>% layout(
  title="Mixed - Control",
  scene=list(
    aspectmode = 'cube',
    xaxis=list(title="Phosphorus", autorange = "reversed"),
    yaxis=list(title="Nitrogen", autorange = "reversed"), 
    zaxis=axz2))


# Mixed (Thelephora+Suillus) vs. Thelephora

0 - min(change.mat.myco) / (max(change.mat.myco) - min(change.mat.myco)) # calculate midpoint value

plot_ly( z = ~change.mat.myco) %>% add_surface(opacity = 1, colorscale = list(
  c(0,  0.679591, 1), c("blue","white", "red"))) %>% layout(
    title="Mixed - Thelephora",
    scene=list(
      aspectmode = 'cube',
      xaxis=list(title="Phosphorus", autorange = "reversed"),
      yaxis=list(title="Nitrogen", autorange = "reversed"), 
      zaxis=axz2))

```
  
  
  
# FIGURE 3  
  
Mycorrhizal colonization responses across N and P fertilization gradients and convex hull niche dimensions.

```{r}

Thelephora.col.alone.mod <- lmer(Thelephora_colonization ~ log(Nitrogen) * log(Phosphorus) + (1|REP), data = thelephora)
Thelephora.col.mixed.mod <- lmer(Thelephora_colonization ~ log(Nitrogen) * log(Phosphorus) + (1|REP), data = mixed)
Suillus.col.mixed.mod <- lmer(Suillus_colonization ~ log(Nitrogen) * log(Phosphorus) + (1|REP), data = mixed)

anova(Thelephora.col.alone.mod)
anova(Thelephora.col.mixed.mod)
anova(Suillus.col.mixed.mod)

# TABLE 2 

# Thelephora (alone)
                              #Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
#log(Nitrogen)                 3198.0  3198.0     1    88 11.4719 0.001058 **
#log(Phosphorus)                 50.0    50.0     1    88  0.1794 0.672954   
#log(Nitrogen):log(Phosphorus)  311.3   311.3     1    88  1.1168 0.293502

# Thelephora (+ S)
                              #Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
#log(Nitrogen)                 613.04  613.04     1 91.006  3.3760 0.06942 .
#log(Phosphorus)                 8.89    8.89     1 91.000  0.0490 0.82535  
#log(Nitrogen):log(Phosphorus)  38.35   38.35     1 91.218  0.2112 0.64693

# Suillus (+ T)
                              #Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
#log(Nitrogen)                 123.45  123.45     1 91.010  0.8381 0.36237  
#log(Phosphorus)               598.74  598.74     1 91.001  4.0649 0.04673 *
#log(Nitrogen):log(Phosphorus) 138.94  138.94     1 91.362  0.9433 0.33400


# Mycorrhizal colonization ~ Nitrogen 

ggplot() +
  geom_point(position = position_jitter(width = 0.1), aes(x=log(Nitrogen), y=Suillus_colonization), color = "#f28e2b", size = 1, data = suillus.mixed.dat, show.legend = F) +
  geom_smooth(aes(x=log(Nitrogen), y=Suillus_colonization), method = "lm", color = "#f28e2b", se=F, lwd=0.75, data = suillus.mixed.dat, show.legend = F) +
  geom_point(position = position_jitter(width = 0.1), aes(x=log(Nitrogen), y=Thelephora_colonization), color = "#4e79a7", fill = "white", pch = 21, size = 1, data = thelephora.mixed.dat, show.legend = F) +
  geom_smooth(aes(x=log(Nitrogen), y=Thelephora_colonization), color = "#4e79a7", method = "lm", linetype = "dashed", se=F, lwd=0.75, data = thelephora.mixed.dat, show.legend = F) +
  geom_point(position = position_jitter(width = 0.1), aes(x=log(Nitrogen), y=Thelephora_colonization), color = "#4e79a7", size = 1, data = thelephora.alone.dat, show.legend = F) +
  geom_smooth(aes(x=log(Nitrogen), y=Thelephora_colonization), color = "#4e79a7", method = "lm", se=F, lwd=0.75, data = thelephora.alone.dat, show.legend = F) +
  xlab("Log Nitrogen (mg/kg soil)") +
  ylab("Mycorrhizal Colonization (%)") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())


# Mycorrhizal colonization ~ Phosphorus

ggplot() +
  geom_point(position = position_jitter(width = 0.1), aes(x=log(Phosphorus), y=Suillus_colonization), color = "#f28e2b", size = 1, data = suillus.mixed.dat, show.legend = F) +
  geom_smooth(aes(x=log(Phosphorus), y=Suillus_colonization), method = "lm", color = "#f28e2b", se=F, lwd=0.75, data = suillus.mixed.dat, show.legend = F) +
  geom_point(position = position_jitter(width = 0.1), aes(x=log(Phosphorus), y=Thelephora_colonization), color = "#4e79a7", fill = "white", pch = 21, size = 1, data = thelephora.mixed.dat, show.legend = F) +
  geom_smooth(aes(x=log(Phosphorus), y=Thelephora_colonization), color = "#4e79a7", method = "lm", linetype = "dashed", se=F, lwd=0.75, data = thelephora.mixed.dat, show.legend = F) +
  geom_point(position = position_jitter(width = 0.1), aes(x=log(Phosphorus), y=Thelephora_colonization), color = "#4e79a7", size = 1, data = thelephora.alone.dat, show.legend = F) +
  geom_smooth(aes(x=log(Phosphorus), y=Thelephora_colonization), color = "#4e79a7", method = "lm", se=F, lwd=0.75, data = thelephora.alone.dat, show.legend = F) +
  xlab("Log Phosphorus (mg/kg soil)") +
  ylab("Mycorrhizal Colonization (%)") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())



###   Constructing mycorrhizal colonization convex hulls and plotting response surfaces

# Create & smooth 7x7 matrices of mycorrhizal responses across N and P gradients

thel.colon.mat <- acast(thelephora, Nitrogen~Phosphorus, value.var = "Thelephora_colonization",mean, drop=TRUE)
mixed.thel.mat <- acast(mixed, Nitrogen~Phosphorus, value.var="Thelephora_colonization",mean,drop=TRUE)
mixed.suillus.mat <- acast(mixed, Nitrogen~Phosphorus, value.var="Suillus_colonization",mean,drop=TRUE)

thel.colon.smooth<-image.smooth(thel.colon.mat, theta=1) #smooths over missing matrix values
mixed.thel.smooth<-image.smooth(mixed.thel.mat, theta=1)
mixed.suillus.smooth<-image.smooth(mixed.suillus.mat, theta=1)


# Response surface plotting

# Thelephora (alone)
plot_ly( z = ~thel.colon.smooth$z) %>% add_surface(opacity = 0.95, colors = viridis::magma(50)) %>% 
  layout(
    title="Thelephora (alone)",
    scene=list(
      aspectmode = 'cube',
      xaxis=list(title="Phosphorus", autorange = "reversed"),
      yaxis=list(title="Nitrogen", autorange = "reversed"), 
      zaxis=list(title="Root Colonization (%)", range = c(30, 80))))

# Thelephora (+S)
plot_ly( z = ~mixed.thel.smooth$z) %>% add_surface(opacity = 0.95,  colors = viridis::magma(50)) %>% layout(
  title="Thelephora (+S)",
  scene=list(
    aspectmode = 'cube',
    xaxis=list(title="Phosphorus", autorange = "reversed"),
    yaxis=list(title="Nitrogen", autorange = "reversed"), 
    zaxis=list(title="Root Colonization (%)", range = c(30, 80))))

# Suillus (+T)
plot_ly( z = ~mixed.suillus.smooth$z) %>% add_surface(opacity = 0.95, colors = viridis::magma(50)) %>% layout(
  title="Suillus (+T)",
  scene=list(
    aspectmode = 'cube',
    xaxis=list(title="Phosphorus", autorange = "reversed"),
    yaxis=list(title="Nitrogen", autorange = "reversed"), 
    zaxis=list(title="Root Colonization (%)", range = c(15, 35))))


###  Random subsampling and bootstrapping approach to create distributions and calculate error estimates for convex hull volumes from each mycorrhizal treatment

## Thelephora (alone)

thelephora.col.mat.mean.hull.vols <- NA

for(i in 1:1000){
  
  thelephora.col_matrix_list <- list()
  volume_output_thelephora.col <- NA
  
  for(j in 1:100){
    original.matrix  <- thelephora
    random.number<-rnorm(length(original.matrix$Nitrogen),0,1)
    sample.order<-order(random.number)
    resampled.matrix<-original.matrix[sample.order[1:round(length(sample.order)*0.5)],]
    sampled.matrix <- acast(resampled.matrix, Nitrogen~Phosphorus, value.var = "Thelephora_colonization", mean, drop=TRUE)
    sampled.mat.smooth<- image.smooth(sampled.matrix, theta=1) #smooths over missing matrix values
    sampled.hull <- convhulln(melt(sampled.mat.smooth$z), output.options=TRUE)
    volume_output_thelephora.col[j]<-sampled.hull$vol
    thelephora.col_matrix_list[[j]] <- sampled.mat.smooth$z
  }
  thelephora.ids.to.remove <- sapply(thelephora.col_matrix_list, function(i) length(i) < 49) 
  thelephora.col_matrix_list <- thelephora.col_matrix_list[!thelephora.ids.to.remove]
  thelephora.col.mat.mean <- Reduce("+", thelephora.col_matrix_list) / length(thelephora.col_matrix_list)
  thelephora.col.mat.mean.hull <- convhulln(melt(thelephora.col.mat.mean), output.options=TRUE)
  thelephora.col.mat.mean.hull.vols[i] <- thelephora.col.mat.mean.hull$vol
}

mean(thelephora.col.mat.mean.hull.vols, na.rm=T)
  sd(thelephora.col.mat.mean.hull.vols, na.rm=T)


## Thelephora (+S)

thelephora.colS.mat.mean.hull.vols <- NA

for(i in 1:1000){
  
  thelephora.colS_matrix_list <- list()
  volume_output_thelephora.colS <- NA
  
  for(j in 1:100){
    original.matrix  <- mixed
    random.number<-rnorm(length(original.matrix$Nitrogen),0,1)
    sample.order<-order(random.number)
    resampled.matrix<-original.matrix[sample.order[1:round(length(sample.order)*0.5)],]
    sampled.matrix <- acast(resampled.matrix, Nitrogen~Phosphorus, value.var = "Thelephora_colonization", mean, drop=TRUE)
    sampled.mat.smooth<- image.smooth(sampled.matrix, theta=1) #smooths over missing matrix values
    sampled.hull <- convhulln(melt(sampled.mat.smooth$z), output.options=TRUE)
    volume_output_thelephora.colS[j]<-sampled.hull$vol
    thelephora.colS_matrix_list[[j]] <- sampled.mat.smooth$z
  }
  thelephora.ids.to.remove <- sapply(thelephora.colS_matrix_list, function(i) length(i) < 49) 
  thelephora.colS_matrix_list <- thelephora.colS_matrix_list[!thelephora.ids.to.remove] 
  thelephora.colS.mat.mean <- Reduce("+", thelephora.colS_matrix_list) / length(thelephora.colS_matrix_list) 
  thelephora.colS.mat.mean.hull <- convhulln(melt(thelephora.colS.mat.mean), output.options=TRUE)
  thelephora.colS.mat.mean.hull.vols[i] <- thelephora.colS.mat.mean.hull$vol
}

mean(thelephora.colS.mat.mean.hull.vols, na.rm=T)
  sd(thelephora.colS.mat.mean.hull.vols, na.rm=T)


## Suillus (+T)

suillus.col.mat.mean.hull.vols <- NA

for(i in 1:1000){
  
  suillus.col_matrix_list <- list()
  volume_output_suillus.col <- NA
  
  for(j in 1:100){
    original.matrix  <- mixed
    random.number<-rnorm(length(original.matrix$Nitrogen),0,1)
    sample.order<-order(random.number)
    resampled.matrix<-original.matrix[sample.order[1:round(length(sample.order)*0.5)],]
    sampled.matrix <- acast(resampled.matrix, Nitrogen~Phosphorus, value.var = "Suillus_colonization", mean, drop=TRUE)
    sampled.mat.smooth<- image.smooth(sampled.matrix, theta=1) #smooths over missing matrix values
    sampled.hull <- convhulln(melt(sampled.mat.smooth$z), output.options=TRUE)
    volume_output_suillus.col[j]<-sampled.hull$vol
    suillus.col_matrix_list[[j]] <- sampled.mat.smooth$z
  }
  suillus.col.ids.to.remove <- sapply(suillus.col_matrix_list, function(i) length(i) < 49) 
  suillus.col_matrix_list <- suillus.col_matrix_list[!suillus.col.ids.to.remove]
  suillus.col.mat.mean <- Reduce("+", suillus.col_matrix_list) / length(suillus.col_matrix_list)
  suillus.col.mat.mean.hull <- convhulln(melt(suillus.col.mat.mean), output.options=TRUE)
  suillus.col.mat.mean.hull.vols[i] <- suillus.col.mat.mean.hull$vol
}

mean(suillus.col.mat.mean.hull.vols, na.rm=T)
sd(suillus.col.mat.mean.hull.vols, na.rm=T)


##

t.test(thelephora.colS.mat.mean.hull.vols, thelephora.col.mat.mean.hull.vols, conf.level=0.95)
t.test(mixed.mat.mean.hull.vols, control.mat.mean.hull.vols, conf.level=0.95)
t.test(mixed.mat.mean.hull.vols, thelephora.mat.mean.hull.vols, conf.level=0.95)

##


# Combine density data for plotting
myc.vol.thel.mix.sui.dat <- data.frame(dens = c(thelephora.col.mat.mean.hull.vols, thelephora.colS.mat.mean.hull.vols, suillus.col.mat.mean.hull.vols), lines = rep(c("Thel alone", "Thel mixed", "Suillus mixed"), each = 1000))

# Thelephora (alone)
ggplot() + 
  geom_histogram(data = subset(myc.vol.thel.mix.sui.dat, lines=="Thel alone"), aes(x = dens), alpha = 1, color = "#4e79a7", fill="#4e79a7", binwidth=10) +
  ylab("Count") +
  xlab("Convex Hull Volume") +
  geom_vline(xintercept = mean(thelephora.col.mat.mean.hull.vols), linetype="dashed") +
  xlim(200, 575) +
  ylim(0,375) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.background = element_blank(), strip.text = element_blank())


# Thelephora (+S)
ggplot() + 
  geom_histogram(data = subset(myc.vol.thel.mix.sui.dat, lines=="Thel mixed"), aes(x = dens), alpha = 1, color = "#4e79a7", fill="white", binwidth=10) +
  ylab("Count") +
  xlab("Convex Hull Volume") +
  geom_vline(xintercept = mean(thelephora.colS.mat.mean.hull.vols), linetype="dashed") +
  xlim(200, 575) +
  ylim(0,375) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.background = element_blank(), strip.text = element_blank())


# Suillus (+T)
ggplot() + 
  geom_histogram(data = subset(myc.vol.thel.mix.sui.dat, lines=="Suillus mixed"), aes(x = dens), alpha = 1, color = "#f28e2b", fill="#f28e2b", binwidth=10) +
  ylab("Count") +
  xlab("Convex Hull Volume") +
  geom_vline(xintercept = mean(suillus.col.mat.mean.hull.vols), linetype="dashed") +
  xlim(200, 575) +
  ylim(0,375) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.background = element_blank(), strip.text = element_blank())

```



# FIGURE 4

Correlation between plant biomass and mycorrhizal colonization

```{r}
TotalBiomass.thel.alone.mod <- lmer(log(TotalBiomass) ~ Thelephora_colonization + (1|REP), data = thelephora)
TotalBiomass.thel.mixed.mod <- lmer(log(TotalBiomass) ~ Thelephora_colonization + (1|REP), data = mixed)
TotalBiomass.suillus.mixed.mod <- lmer(log(TotalBiomass) ~ Suillus_colonization + (1|REP), data = mixed)
  
# TABLE 3
summary(lm(TotalBiomass ~ Thelephora_colonization, data = thelephora))
summary(lm(TotalBiomass ~ Thelephora_colonization, data = mixed))
summary(lm(TotalBiomass ~ Suillus_colonization, data = mixed))


# Plotting (linear regression of biomass ~ colonization)

Myc_trts2 <- factor(pinus.myc.niche.dat$MYC_recode, levels = c("Thelephora", "Suillus_Thelephora")) # re-arranging the order of these groups for myc colonization plotting

ggplot(aes(x=Thelephora_colonization, y=log(TotalBiomass)), data = pinus.myc.niche.dat) +
  geom_point(aes(x=Suillus_colonization, y=log(TotalBiomass)), color = "#f28e2b", size = 1.5, data = subset(pinus.myc.niche.dat, Suillus_colonization > 0), show.legend = F, na.rm = T) +
  geom_smooth(aes(x=Suillus_colonization, y=log(TotalBiomass)), color = "#f28e2b", method = "lm", se=F, lwd=1, data = subset(pinus.myc.niche.dat, Suillus_colonization > 0), show.legend = F, na.rm = T) +
  geom_point(aes(x=Thelephora_colonization, y=log(TotalBiomass), color = Myc_trts2, fill = Myc_trts2), pch = 21, size = 1.5, data = pinus.myc.niche.dat, show.legend = F, na.rm = T) +
  geom_smooth(aes(x=Thelephora_colonization, y=log(TotalBiomass), color = Myc_trts2, linetype = Myc_trts2), method = "lm", se=F, lwd=1, data = pinus.myc.niche.dat, show.legend = F, na.rm = T) +
  scale_color_manual(values = c("#4e79a7", "#4e79a7")) +
  scale_fill_manual(values = c("#4e79a7", "white")) +
  ylab("Log Total Plant Biomass (mg)") +
  xlab("Mycorrhizal Colonization (%)") +
  xlim(15, 85) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

```



# Other analyses
```{r}

# Correlation between Thelephora and Suillus root colonization levels

Suillus.Thelephora.corr.mod <- lm(Suillus_colonization ~ Thelephora_colonization, data = subset(pinus.myc.niche.dat, MYC == "Suillus_Thelephora" & Suillus_colonization > 0)) # Subsetting to only treatments where both fungi were introduced at same time

anova(Suillus.Thelephora.corr.mod)
summary(Suillus.Thelephora.corr.mod)  #adj R-squared = 0.60

plot(Suillus_colonization ~ Thelephora_colonization, data = subset(pinus.myc.niche.dat, MYC == "Suillus_Thelephora" & Suillus_colonization > 0))

```

