---
title: "Plant-mycorrhizal niche mapping analysis and figures"
output: html_notebook
---
Created by: Michael Van Nuland & Kabir Peay (October 2019 - February 2020); *REVISED April 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(processx)
library(fields)
library(reshape2)
library(ggplot2)
library(viridis)
library(RColorBrewer)

pinus.myc.niche.dat <- read.csv(file="~/Symbiotic niche mapping analysis and figures/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$Total_colonization < 2)

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

suillus <- subset(pinus.myc.niche.dat, pinus.myc.niche.dat$MYC=="Suillus")

mixed <- subset(pinus.myc.niche.dat, pinus.myc.niche.dat$MYC=="Suillus_Thelephora") 

pinus.myc.niche.dat.redef <- as.data.frame(rbind(control, thelephora, suillus, mixed))
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 * log(Nitrogen) * log(Phosphorus) + (1|REP), data = pinus.myc.niche.dat.redef)
  anova(TotalBiomass.mod, type="III")

# TABLE 1
                                    #Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
#MYC                                0.196   0.065     3 305.01   0.7807   0.50548    
#log(Nitrogen)                     37.576  37.576     1 305.03 448.5069 < 2.2e-16 ***
#log(Phosphorus)                    0.013   0.013     1 305.01   0.1578   0.69147    
#MYC:log(Nitrogen)                  0.421   0.140     3 305.01   1.6741   0.17256    
#MYC:log(Phosphorus)                0.718   0.239     3 305.01   2.8572   0.03728 *  
#log(Nitrogen):log(Phosphorus)      5.476   5.476     1 305.17  65.3592 1.478e-14 ***
#MYC:log(Nitrogen):log(Phosphorus)  0.221   0.074     3 305.12   0.8792   0.45218

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

# Pull out residuals for plotting

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, levels = c("Control", "Thelephora", "Suillus", "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 = 0.4, data = pinus.myc.niche.dat.redef, show.legend = F) +
  geom_smooth(aes(x=log(Nitrogen), y=P.resids, color = Myc_trts), method = "lm", se=F, lwd=0.4, data=pinus.myc.niche.dat.redef, show.legend = F) +
  scale_color_manual(values = c("grey55","#00AFBB", "#E7B800", "#FC4E07")) +
  ylim(-1.5, 1.10) +
  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 = 0.4, data=pinus.myc.niche.dat.redef, show.legend = F) +
  geom_smooth(aes(x=log(Phosphorus), y=N.resids, color = Myc_trts), method ="lm", se=F, lwd=0.4, data=pinus.myc.niche.dat.redef, show.legend = F) +
  scale_color_manual(values = c("grey55","#00AFBB", "#E7B800", "#FC4E07")) +
  ylim(-1.5, 1.10) +
  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)
suil.mat <- acast(suillus, 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)
suil.mat.smooth <- image.smooth(suil.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     #7467.71

suil.hull <- convhulln(melt(suil.mat.smooth$z), output.options=TRUE)
  suil.hull$vol     #7535.753
  
mixed.hull <- convhulln(melt(mixed.mat.smooth$z), output.options=TRUE)
  mixed.hull$vol    #5231.45
  
  
# Response surface plotting

axz <- list(title="Total biomass (g)", range = c(180, 1180)) #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))

# Suillus
plot_ly( z = ~suil.mat.smooth$z) %>% add_surface(opacity = 0.95) %>% layout(
  title="Suillus",
  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)

  
## SUILLUS

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

for(i in 1:1000){
  
  suillus_matrix_list <- list()
  volume_output_suillus <- NA
  
  for(j in 1:100){
    original.matrix<-suillus
    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_suillus[j]<-sampled.hull$vol
    suillus_matrix_list[[j]] <- sampled.mat.smooth$z
  }
  suillus.ids.to.remove <- sapply(suillus_matrix_list, function(i) length(i) < 49) 
  suillus_matrix_list <- suillus_matrix_list[!suillus.ids.to.remove]
  suillus.mat.mean <- Reduce("+", suillus_matrix_list) / length(suillus_matrix_list) # Calculate means by x/y elements
  suillus.mat.mean.hull <- convhulln(melt(suillus.mat.mean), output.options=TRUE)
  suillus.mat.mean.hull.vols[i] <- suillus.mat.mean.hull$vol
}

mean(suillus.mat.mean.hull.vols, na.rm=T)
  sd(suillus.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(thelephora.mat.mean.hull.vols, suillus.mat.mean.hull.vols, conf.level=0.95)

t.test(suillus.mat.mean.hull.vols, control.mat.mean.hull.vols, conf.level=0.95)
t.test(suillus.mat.mean.hull.vols, thelephora.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)
t.test(mixed.mat.mean.hull.vols, suillus.mat.mean.hull.vols, conf.level=0.95)

##
  
# Combine the density distributions for plotting
vol.control.thel.suil.mix.dat <- data.frame(dens = c(control.mat.mean.hull.vols, thelephora.mat.mean.hull.vols, suillus.mat.mean.hull.vols, mixed.mat.mean.hull.vols), lines = rep(c("control", "thel", "suil", "mix"), each = 1000))
  
# Combined plot of bootstrapped convex hulls

ggplot() + 
  geom_histogram(data = subset(vol.control.thel.suil.mix.dat, lines=="control"), aes(x = dens), alpha = 0.75, color = "grey55", fill="grey55", binwidth=50) +
  geom_vline(xintercept = mean(control.mat.mean.hull.vols), linetype="dashed") +
  geom_histogram(data = subset(vol.control.thel.suil.mix.dat, lines=="thel"), aes(x = dens), alpha = 0.75, color = "#00AFBB", fill="#00AFBB", binwidth=50) +
  geom_vline(xintercept = mean(thelephora.mat.mean.hull.vols), linetype="dashed") +
  geom_histogram(data = subset(vol.control.thel.suil.mix.dat, lines=="suil"), aes(x = dens), alpha = 0.75, color = "#E7B800", fill="#E7B800", binwidth=50) +
  geom_vline(xintercept = mean(suillus.mat.mean.hull.vols), linetype="dashed") +
  geom_histogram(data = subset(vol.control.thel.suil.mix.dat, lines=="mix"), aes(x = dens), alpha = 0.75, color = "#FC4E07", fill="#FC4E07", binwidth=50) +
  geom_vline(xintercept = mean(mixed.mat.mean.hull.vols), linetype="dashed") +
  ylab("Count") +
  xlab("Bootstrapped Convex Hull Volume") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 210)) +
  xlim(4200, 7500) +
  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.suil <- 100*((suil.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.suil.thel <- 100*((suil.mat.smooth$z-thel.mat.smooth$z)/thel.mat.smooth$z)
change.mat.mix.thel <- 100*((mixed.mat.smooth$z-thel.mat.smooth$z)/thel.mat.smooth$z)
change.mat.mix.suil <- 100*((mixed.mat.smooth$z-suil.mat.smooth$z)/suil.mat.smooth$z)


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

# 0 - min(change.mat.thel) / (max(change.mat.thel) - min(change.mat.thel)) = example for calculating mid-point color values in % change plots


# Thelephora vs. Control

plot_ly( z = ~change.mat.thel) %>% add_surface(opacity = 1, colorscale = list(
  c(0, (0 - min(change.mat.thel) / (max(change.mat.thel) - min(change.mat.thel))), 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))

# Suillus vs. Control

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

# Mixed (Thelephora+Suillus) vs. Control

plot_ly( z = ~change.mat.mixed) %>% add_surface(opacity = 1, colorscale = list(
  c(0, (0 - min(change.mat.mixed) / (max(change.mat.mixed) - min(change.mat.mixed))), 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))

# Suillus vs. Thelephora

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

# Mixed vs. Thelephora

plot_ly( z = ~change.mat.mix.thel) %>% add_surface(opacity = 1, colorscale = list(
  c(0,  (0 - min(change.mat.mix.thel) / (max(change.mat.mix.thel) - min(change.mat.mix.thel))), 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))

# Mixed vs. Suillus

plot_ly( z = ~change.mat.mix.suil) %>% add_surface(opacity = 1, colorscale = list(
  c(0,  (0 - min(change.mat.mix.suil) / (max(change.mat.mix.suil) - min(change.mat.mix.suil))), 1), c("blue","white", "red"))) %>% layout(
    title="Mixed - Suillus",
    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}

Myc.col.mod <- lmer(Total_colonization ~ MYC* log(Nitrogen) * log(Phosphorus) + (1|REP), data=subset(pinus.myc.niche.dat.redef, pinus.myc.niche.dat.redef$MYC!="Control"))
  anova(Myc.col.mod)

# TABLE 2
                                  #Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
#MYC                               4524.1 2262.04     2 262.02  9.9937 6.567e-05 ***
#log(Nitrogen)                     2785.4 2785.38     1 262.00 12.3058 0.0005312 ***
#log(Phosphorus)                    524.5  524.48     1 262.00  2.3172 0.1291594    
#MYC:log(Nitrogen)                 1161.3  580.63     2 262.01  2.5652 0.0788329 .  
#MYC:log(Phosphorus)               1429.2  714.62     2 262.01  3.1572 0.0441682 *  
#log(Nitrogen):log(Phosphorus)       19.2   19.18     1 262.00  0.0847 0.7712259    
#MYC:log(Nitrogen):log(Phosphorus)  396.6  198.30     2 262.01  0.8761 0.4176332 


Thelephora.col.mod <- lm(Total_colonization ~ log(Nitrogen) * log(Phosphorus) + (1|REP), data = thelephora)
  anova(Thelephora.col.mod)
                              #Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
#log(Nitrogen)                 3430.7  3430.7     1    90 10.0516 0.002081 **
#log(Phosphorus)                327.5   327.5     1    90  0.9594 0.329968   
#log(Nitrogen):log(Phosphorus)  358.6   358.6     1    90  1.0508 0.308079
  
Suillus.col.mod <- lmer(Total_colonization ~ log(Nitrogen) * log(Phosphorus) + (1|REP), data = suillus)
  anova(Suillus.col.mod)
                              #Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
#log(Nitrogen)                 139.28  139.28     1 84.000  0.8167 0.3687  
#log(Phosphorus)               902.53  902.53     1 84.006  5.2920 0.0239 *
#log(Nitrogen):log(Phosphorus)  38.49   38.49     1 84.042  0.2257 0.6360
  
Mixed.col.mod <- lmer(Total_colonization ~ log(Nitrogen) * log(Phosphorus) + (1|REP), data = mixed)
  anova(Mixed.col.mod)
                              #Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
#log(Nitrogen)                 484.02  484.02     1 89.003  3.0548 0.08395 .
#log(Phosphorus)               757.95  757.95     1 89.000  4.7836 0.03135 *
#log(Nitrogen):log(Phosphorus)  25.46   25.46     1 89.000  0.1607 0.68948

  
# Mycorrhizal colonization ~ Nitrogen 

ggplot(aes(x=log(Nitrogen), y=Total_colonization), data = pinus.myc.niche.dat.redef) +
  geom_point(position = position_jitterdodge(dodge.width = 0.3, jitter.width = 0.1), aes(x=log(Nitrogen), y=Total_colonization, color = Myc_trts), size = 0.55, data = pinus.myc.niche.dat.redef, show.legend = T) +
  geom_smooth(aes(x=log(Nitrogen), y=Total_colonization, color = Myc_trts), method = "lm", se=F, lwd=0.5, data=pinus.myc.niche.dat.redef, show.legend = F) +
  scale_color_manual(values = c("grey55","#00AFBB", "#E7B800", "#FC4E07")) +
  ylim(3,95) +
  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(aes(x=log(Phosphorus), y=Total_colonization), data = pinus.myc.niche.dat.redef) +
  geom_point(position = position_jitterdodge(dodge.width = 0.3, jitter.width = 0.1), aes(x=log(Phosphorus), y=Total_colonization, color = Myc_trts), size = 0.55, data = pinus.myc.niche.dat.redef, show.legend = F) +
  geom_smooth(aes(x=log(Phosphorus), y=Total_colonization, color = Myc_trts), method = "lm", se=F, lwd=0.5, data=pinus.myc.niche.dat.redef, show.legend = F) +
  scale_color_manual(values = c("grey55","#00AFBB", "#E7B800", "#FC4E07")) +
  ylim(3,95) +
  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.col.mat <- acast(thelephora, Nitrogen~Phosphorus, value.var = "Total_colonization", mean, drop=TRUE)
suil.col.mat <- acast(suillus, Nitrogen~Phosphorus, value.var="Total_colonization", mean, drop=TRUE)
mixed.col.mat <- acast(mixed, Nitrogen~Phosphorus, value.var="Total_colonization", mean, drop=TRUE)

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


# Response surface plotting

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

# Suillus
plot_ly( z = ~suil.col.mat.smooth$z) %>% add_surface(opacity = 0.95,  colors = viridis::magma(50)) %>% layout(
  title="Suillus",
  scene=list(
    aspectmode = 'cube',
    xaxis=list(title="Phosphorus", autorange = "reversed"),
    yaxis=list(title="Nitrogen", autorange = "reversed"), 
    zaxis=list(title="Myc Colonization (%)", range = c(40, 85))))

# Mixed
plot_ly( z = ~mixed.col.mat.smooth$z) %>% add_surface(opacity = 0.95, colors = viridis::magma(50)) %>% layout(
  title="Mixed",
  scene=list(
    aspectmode = 'cube',
    xaxis=list(title="Phosphorus", autorange = "reversed"),
    yaxis=list(title="Nitrogen", autorange = "reversed"), 
    zaxis=list(title="Myx Colonization (%)", range = c(40, 85))))


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

## Thelephora

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 = "Total_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)


## Suillus

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  <- suillus
    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 = "Total_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.ids.to.remove <- sapply(suillus.col_matrix_list, function(i) length(i) < 49) 
  suillus.col_matrix_list <- suillus.col_matrix_list[!suillus.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)


## Mixed

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

for(i in 1:1000){
  
  mixed.col_matrix_list <- list()
  volume_output_mixed.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 = "Total_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_mixed.col[j]<-sampled.hull$vol
    mixed.col_matrix_list[[j]] <- sampled.mat.smooth$z
  }
  mixed.col.ids.to.remove <- sapply(mixed.col_matrix_list, function(i) length(i) < 49) 
  mixed.col_matrix_list <- mixed.col_matrix_list[!mixed.col.ids.to.remove]
  mixed.col.mat.mean <- Reduce("+", mixed.col_matrix_list) / length(mixed.col_matrix_list)
  mixed.col.mat.mean.hull <- convhulln(melt(mixed.col.mat.mean), output.options=TRUE)
  mixed.col.mat.mean.hull.vols[i] <- mixed.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)
mean(suillus.col.mat.mean.hull.vols, na.rm=T); sd(suillus.col.mat.mean.hull.vols, na.rm=T)
mean(mixed.col.mat.mean.hull.vols, na.rm=T); sd(mixed.col.mat.mean.hull.vols, na.rm=T)

##

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

##

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

#

ggplot() + 
  geom_histogram(data = subset(myc.vol.thel.suil.mix.dat, lines=="Thelephora"), aes(x = dens), alpha = 0.75, color = "#00AFBB", fill="#00AFBB", binwidth=6) +
  geom_histogram(data = subset(myc.vol.thel.suil.mix.dat, lines=="Suillus"), aes(x = dens), alpha = 0.75, color = "#E7B800", fill="#E7B800", binwidth=6) +
  geom_histogram(data = subset(myc.vol.thel.suil.mix.dat, lines=="Mixed"), aes(x = dens), alpha = 0.75, color = "#FC4E07", fill="#FC4E07", binwidth=6) +
  ylab("Count") +
  xlab("Convex Hull Volume") +
  geom_vline(xintercept = mean(thelephora.col.mat.mean.hull.vols), linetype="dashed") +
  geom_vline(xintercept = mean(suillus.col.mat.mean.hull.vols), linetype="dashed") +
  geom_vline(xintercept = mean(mixed.col.mat.mean.hull.vols), linetype="dashed") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 200)) +
  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.MycCol.mod <- lmer(log(TotalBiomass) ~ MYC * Total_colonization + (1|REP), data = subset(pinus.myc.niche.dat.redef, pinus.myc.niche.dat.redef$MYC!="Control"))
  anova(TotalBiomass.MycCol.mod)

# TABLE 3
                        #Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
#MYC                    0.4744  0.2372     2 268.62  0.9822    0.3758    
#Total_colonization     5.6237  5.6237     1 268.08 23.2895 2.342e-06 ***
#MYC:Total_colonization 0.2184  0.1092     2 268.53  0.4523    0.6366 


# Plotting (linear regression of biomass ~ colonization)

ggplot() +
  geom_point(aes(x=Total_colonization, y=log(TotalBiomass), color = MYC), data = subset(pinus.myc.niche.dat.redef, pinus.myc.niche.dat.redef$MYC!="Control"), size = 1.25, show.legend = F) +
  #geom_smooth(aes(x=Total_colonization, y=log(TotalBiomass), color=MYC), data = subset(pinus.myc.niche.dat.redef, pinus.myc.niche.dat.redef$MYC!="Control"), method = "lm", se=F, lwd=0.75, show.legend = F) +
  geom_smooth(aes(x=Total_colonization, y=log(TotalBiomass)), color="black", data = subset(pinus.myc.niche.dat.redef, pinus.myc.niche.dat.redef$MYC!="Control"), method = "lm", se=F, lwd=1.25, show.legend = F) +
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07")) +
  xlim(3,95) +
  ylab("Log Total Plant Biomass (mg)") +
  xlab("Mycorrhizal Colonization (% root length)") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

```

