library(vegan)
library(ggplot2)
library(mobr)
library(grid)
library(gridExtra)
library(nlme)
library(car)
library(dplyr)

community<-read.csv("McKibben_Henning_2018_plant_community_data.csv")

names(community)
str(community)

community2<- community[,8:71]
names(community2)

rich<-specnumber(community2) # species richness
ENSPIE<-calc_PIE(community2, ENS=TRUE) # evenness w/ richness
PIE<-calc_PIE(community2) # evenness

shan<-diversity(community2)
simp<-diversity(community2, "simp")

braydiss = vegdist(community2, "bray")

diversitymetrics<-cbind(rich,ENSPIE, PIE, shan, simp)
divtable<-cbind(community[,1:7],diversitymetrics)
names(divtable)

lme.simp <- lme(simp~Treatment*Elevation, random=~1|Pairing, data=divtable,  method="REML")
summary(lme.simp)
Anova(lme.simp)

lme.shan <- lme(shan~Treatment*Elevation, random=~1|Pairing, data=divtable,  method="REML")
summary(lme.shan)
Anova(lme.simp)


grouped <- group_by(divtable, Elevation)
summarise(grouped, meanrich=mean(rich), meaneven=mean(PIE))


############################################################
##                   plant richness                       ##
############################################################
## Model selection using lme
lme.rich <- lme(rich~Treatment*Elevation, random=~1|Pairing, data=divtable,  method="ML")
lme.rich1 <- lme(rich~Treatment, random=~1|Pairing, data=divtable,  method="ML")
lme.rich2 <- lme(rich~Elevation, random=~1|Pairing, data=divtable,  method="ML")
lme.rich3 <- lme(rich~1, random=~1|Pairing, data=divtable,  method="ML")
anova(lme.rich,lme.rich1,lme.rich2,lme.rich3)
# best fit model
lme.rich <- lme(rich~Treatment*Elevation, random=~1|Pairing, data=divtable,  method="REML")
summary(lme.rich)
Anova(lme.rich)

mod1<-lm(rich~Treatment*Elevation, data=divtable)
Anova(mod1)

# Richness figure
richness<-ggplot(data=divtable, aes(x=factor(Elevation), y=rich, fill = Treatment))
richfig<- richness + geom_boxplot() +  theme_bw(base_size=14)+ scale_fill_manual(values=c( "dodgerblue2","darkorange"))+ theme(legend.position="none")



#############################################################################
##              Probability of intraspecfic encounter - evenness           ##
#############################################################################
## Model selection using lme
lme.PIE <- lme(PIE~Treatment*Elevation, random=~1|Pairing, data=divtable,  method="ML")
lme.PIE1 <- lme(PIE~Treatment, random=~1|Pairing, data=divtable,  method="ML")
lme.PIE2 <- lme(PIE~Elevation, random=~1|Pairing, data=divtable,  method="ML")
lme.PIE3 <- lme(PIE~1, random=~1|Pairing, data=divtable,  method="ML")
anova(lme.PIE,lme.PIE1,lme.PIE2,lme.PIE3)

# best fit model
lme.PIE <- lme(PIE~Treatment*Elevation, random=~1|Pairing, data=divtable,  method="REML")
summary(lme.PIE)
Anova(lme.PIE)

divtableC <- subset(divtable, Treatment=="C" )
divtableN <- subset(divtable, Treatment=="N" )
mean(divtableC$PIE)
mean(divtableN$PIE)

PIE_even<-ggplot(data=divtable, aes(x=factor(Elevation), y=PIE, fill = Treatment))
PIE_even<- PIE_even + geom_boxplot() +  theme_bw(base_size=14)+ scale_fill_manual(values=c( "dodgerblue2","darkorange"))+ theme(legend.position="none")

p1 <- ggplot_gtable(ggplot_build(richfig))
p2 <- ggplot_gtable(ggplot_build(PIE_even))


maxWidth = unit.pmax(p1$widths[2:3], p2$widths[2:3])

p1$widths[2:3] <- maxWidth
p2$widths[2:3] <- maxWidth


tiff("Henning_McKibben_Figure1.tiff", width = 4, height = 6, units = 'in', res = 300, compression = 'none')
grid.arrange(p1,p2, ncol = 1)
dev.off()



mod1<-lm(invsimp~Site*Treatment, data=divtable)
Anova(mod1)
summary(mod1)

mod2<-lm(shan~Site*Treatment, data=divtable)
Anova(mod2)
summary(mod2)
boxplot(divtable$shan~divtable$Treatment*divtable$Site)

mod3<-lm(even~Site*Treatment, data=divtable)
Anova(mod3)
summary(mod3)
boxplot(divtable$even~divtable$Treatment*divtable$Site)

mod4<-lm(simp~Site*Treatment, data=divtable)
Anova(mod4)
summary(mod4)
boxplot(divtable$simp~divtable$Treatment*divtable$Site)


###########################################################
#     COMMUNITY ANALYSIS 
###########################################################
library(vegan)

# abundance weighted Bray Curtis
totalperm<-adonis(community2 ~ Elevation*Treatment*Pairing,data=community, permutations=99) #significant Elevation*Treatment and Elevation*Pairing interaction, so let's look at relationships at each site

# presence/absence Jaccard's distance
community2J<- community2 
community2J[community2J>0] <-1
modcommunity2J<-adonis(community2J ~ Elevation*Treatment*Pairing,data=community, permutations=99)

#################################################
#        almont low (weather station) - ELEVATION - 2480m
#################################################

Alow<-subset(community, Site=="Alm_low")
names(Alow)
Alow2<-Alow[,c("UNKNOWN2", "AGOGLA",  "ARTTRI", "BALSAG" ,"CAPBUR"  ,"CARALB", "DRAAUR",  "ELYELY" , "FESTHU",  "IPOAGG", "MUHMON" , "ERISPE"  ,"POA_SP" , "POLPAT",  "CHRVIS" , "SYMROT", "TAROFF")]

popcols<-c(rep("dodgerblue2", 10), rep("darkorange", 10))
LowAlmdis<-vegdist(Alow2, method="bray")
Alow
Treatment<-c("P", "P","P", "P","P","P", "P","P", "P","P","A","A","A","A","A","A","A","A","A","A")
Alow$Treatment<-c("P", "P","P", "P","P","P", "P","P", "P","P","A","A","A","A","A","A","A","A","A","A")
mod<-betadisper(LowAlmdis, group=Alow$Treatment)
anova(mod)

# bray curtis
Alowperm<-adonis(Alow2 ~ Treatment*Pairing  ,data=Alow, permutations=99)

# jaccard's
Alow2J<- Alow2 
Alow2J[Alow2J>0] <-1
modAlow2J<-adonis(Alow2J ~ Elevation*Treatment*Pairing,data=Alow, permutations=99)


AlLow.mds <- metaMDS(Alow2, trymax = 1000)
data.scores <- as.data.frame(scores(AlLow.mds))  #Using the scores function from vegan to extract the site scores and convert to a data.frame
data.scores$site <- rownames(data.scores)  # create a column of site names, from the rownames of data.scores
data.scores$Treatment <- Treatment  #  add the Treatment variable created earlier
head(data.scores)  #look at the data

species.scores <- as.data.frame(scores(AlLow.mds, "species"))  #Using the scores function from vegan to extract the species scores and convert to a data.frame
species.scores$species <- rownames(species.scores)  # create a column of species, from the rownames of species.scores
head(species.scores)  #look at the data

Treatment.a <- data.scores[data.scores$Treatment == "A", ][chull(data.scores[data.scores$Treatment == 
                                                                       "A", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment A
Treatment.p <- data.scores[data.scores$Treatment == "P", ][chull(data.scores[data.scores$Treatment == 
                                                                       "P", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment B

hull.data <- rbind(Treatment.a, Treatment.p)  #combine Treatment.a and Treatment.b
hull.data

ALLOW.MDS<-ggplot() + 
  geom_polygon(data=hull.data,aes(x=NMDS1,y=NMDS2,fill=Treatment,group=Treatment),alpha=0.30) + # add the convex hulls
  geom_text(data=species.scores,aes(x=NMDS1,y=NMDS2,label=species)) +  # add the species labels
  geom_point(data=data.scores,aes(x=NMDS1,y=NMDS2,shape=Treatment,colour=Treatment),size=4) + # add the point markers
  scale_colour_manual(values=c("A" = "darkorange", "P" = "dodgerblue2")) +
  theme_bw() + 
  theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        axis.title.x = element_text(size=18), # remove x-axis labels
        axis.title.y = element_text(size=18), # remove y-axis labels
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())
#################################################
#        almont obs  ,    Elevation = 2740m
#################################################

AObs<-subset(community, Site=="Alm_obs")
names(AObs)
AObs2<-AObs[,c("ALLSEP","AST_SP","CARALB","DELNUT","DRAAUR" , "ERIGLA","ERISPE", "FESTHU" ,"ADELEW","PHLOX","HELQUI" ,"MAHAQU" ,"MUHMON","BALSAG","PHASER","POA_SP", "POTGRA","POLPAT", "ANTROS","CHRVIS" , "ARTTRI", "SEDLAN","SYMROT" ,"SENCRA","TAROFF","UNKNOWN1")]

# bray-curtis
AObsperm<-adonis(AObs2 ~ Treatment*Pairing  ,data=AObs, permutations=99)
# jaccard's
AObs2J<- AObs2 
AObs2J[AObs2J>0] <-1
modAObs2J<-adonis(AObs2J ~ Elevation*Treatment*Pairing,data=AObs, permutations=99)



sim <-  simper(AObs2, group= AObs$Treatment  )
summary(sim)


popcols<-c(rep("dodgerblue2", 10), rep("darkorange", 10))

LowAlmdis<-vegdist(AObs2, method="bray")
AObs
AObs$Treatment<-c("P", "P","P", "P","P","P", "P","P", "P","P","A","A","A","A","A","A","A","A","A","A")
mod<-betadisper(LowAlmdis, group=AObs$Treatment)
anova(mod)

AlmOb.mds <- metaMDS(AObs2, trymax = 1000)
data.scores <- as.data.frame(scores(AlmOb.mds))  #Using the scores function from vegan to extract the site scores and convert to a data.frame
data.scores$site <- rownames(data.scores)  # create a column of site names, from the rownames of data.scores
data.scores$Treatment <- Treatment  #  add the Treatment variable created earlier
head(data.scores)  #look at the data

species.scores <- as.data.frame(scores(AlmOb.mds, "species"))  #Using the scores function from vegan to extract the species scores and convert to a data.frame
species.scores$species <- rownames(species.scores)  # create a column of species, from the rownames of species.scores
head(species.scores)  #look at the data

Treatment.a <- data.scores[data.scores$Treatment == "A", ][chull(data.scores[data.scores$Treatment == 
                                                                       "A", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment A
Treatment.p <- data.scores[data.scores$Treatment == "P", ][chull(data.scores[data.scores$Treatment == 
                                                                       "P", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment B

hull.data <- rbind(Treatment.a, Treatment.p)  #combine Treatment.a and Treatment.b
hull.data

ALOBS.MDS<-ggplot() + 
  geom_polygon(data=hull.data,aes(x=NMDS1,y=NMDS2,fill=Treatment,group=Treatment),alpha=0.30) + # add the convex hulls
  geom_text(data=species.scores,aes(x=NMDS1,y=NMDS2,label=species)) +  # add the species labels
  geom_point(data=data.scores,aes(x=NMDS1,y=NMDS2,shape=Treatment,colour=Treatment),size=4) + # add the point markers
  scale_colour_manual(values=c("A" = "darkorange", "P" = "dodgerblue2")) +
  theme_bw() + 
  theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        axis.title.x = element_text(size=18), # remove x-axis labels
        axis.title.y = element_text(size=18), # remove y-axis labels
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())
#
#################################################
#                Elkton    Elevation = 3200
#################################################

elk<-subset(community, Site=="Elk")
names(elk)
elk2<-elk[,c("AGOGLA","Annual","BROINE","CLATON","UNKNOWN1", "DANTER","DELBAR","ELYELY","ADELEW" , "ERYGRA","HELQUI","HYDROP", "IPOAGG","LATLEU","LIGPOR", "MERFUS","ERISPE","POA_SP","POTGRA","ROSWOO","SENCRA","HELMUL","THAFEN","VIOADU","VICAME","FRAVIR")]

colnames(elk2)[colnames(elk2)=="Annual"] <- "ANDSEP"

elkperm<-adonis(elk2 ~ Treatment*Pairing  ,data=elk, permutations=99)

# jaccard's
elk2J<- elk2 
elk2J[elk2J>0] <-1
modelk2J<-adonis(elk2J ~ Elevation*Treatment*Pairing,data=elk, permutations=99)


elk2dis<-vegdist(elk2, method="bray")
elk
elk$Treatment<-c("P", "P","P", "P","P","P", "P","P", "P","P","A","A","A","A","A","A","A","A","A","A")
mod<-betadisper(elk2dis, group=elk$Treatment)
anova(mod)

Elk.mds <- metaMDS(elk2, trymax = 1000)
data.scores <- as.data.frame(scores(Elk.mds))  #Using the scores function from vegan to extract the site scores and convert to a data.frame
data.scores$site <- rownames(data.scores)  # create a column of site names, from the rownames of data.scores
data.scores$Treatment <- Treatment  #  add the Treatment variable created earlier
head(data.scores)  #look at the data

species.scores <- as.data.frame(scores(Elk.mds, "species"))  #Using the scores function from vegan to extract the species scores and convert to a data.frame
species.scores$species <- rownames(species.scores)  # create a column of species, from the rownames of species.scores
head(species.scores)  #look at the data

Treatment.a <- data.scores[data.scores$Treatment == "A", ][chull(data.scores[data.scores$Treatment == 
                                                                       "A", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment A
Treatment.p <- data.scores[data.scores$Treatment == "P", ][chull(data.scores[data.scores$Treatment == 
                                                                       "P", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment B

hull.data <- rbind(Treatment.a, Treatment.p)  #combine Treatment.a and Treatment.b
hull.data

Elk.MDS<-ggplot() + 
  geom_polygon(data=hull.data,aes(x=NMDS1,y=NMDS2,fill=Treatment,group=Treatment),alpha=0.30) + # add the convex hulls
  geom_text(data=species.scores,aes(x=NMDS1,y=NMDS2,label=species)) +  # add the species labels
  geom_point(data=data.scores,aes(x=NMDS1,y=NMDS2,shape=Treatment,colour=Treatment),size=4) + # add the point markers
  scale_colour_manual(values=c("A" = "darkorange", "P" = "dodgerblue2")) +
  theme_bw() + 
  theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        axis.title.x = element_text(size=18), # remove x-axis labels
        axis.title.y = element_text(size=18), # remove y-axis labels
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())
#

#################################################
#                Painter Boy Mine - ELEVATION = 3392
#################################################

PB<-subset(community, Site=="PB")
names(PB)
PB2<-PB[,c("ADELEW","AGOGLA","Annual","AQUCOL","BOESTR","BROINE","CARALB", "COLLIN", "DRAAUR","ERISPE","ELYELY", "ERYGRA", "CHAANG","FRAVIR" ,"GALSEP" , "HELQUI","IRIGRA","LIGPOR", "MERFUS","POA_SP", "POTGRA","UNKNOWN1" ,"SENCRA", "TAROFF", "THAFEN","VIOADU")]

colnames(PB2)[colnames(PB2)=="Annual"] <- "ANDSEP"

# bray-curtis
PBperm<-adonis(PB2 ~ Treatment*Pairing  ,data=PB, permutations=99)

# jaccard's
PB2J<- PB2 
PB2J[PB2J>0] <-1
modPB2J<-adonis(PB2J ~ Elevation*Treatment*Pairing,data=PB, permutations=99)



PB2dis<-vegdist(PB2, method="bray")
PB
PB$Treatment<-c("P", "P","P", "P","P","P", "P","P", "P","P","A","A","A","A","A","A","A","A","A","A")
mod<-betadisper(PB2dis, group=PB$Treatment)
anova(mod)

pb.mds <- metaMDS(PB2, trymax = 1000)
data.scores <- as.data.frame(scores(pb.mds))  #Using the scores function from vegan to extract the site scores and convert to a data.frame
data.scores$site <- rownames(data.scores)  # create a column of site names, from the rownames of data.scores
data.scores$Treatment <- Treatment  #  add the Treatment variable created earlier
head(data.scores)  #look at the data

species.scores <- as.data.frame(scores(pb.mds, "species"))  #Using the scores function from vegan to extract the species scores and convert to a data.frame
species.scores$species <- rownames(species.scores)  # create a column of species, from the rownames of species.scores
head(species.scores)  #look at the data

Treatment.a <- data.scores[data.scores$Treatment == "A", ][chull(data.scores[data.scores$Treatment == 
                                                                       "A", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment A
Treatment.p <- data.scores[data.scores$Treatment == "P", ][chull(data.scores[data.scores$Treatment == 
                                                                       "P", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment B

hull.data <- rbind(Treatment.a, Treatment.p)  #combine Treatment.a and Treatment.b
hull.data

pb.MDS<-ggplot() + 
  geom_polygon(data=hull.data,aes(x=NMDS1,y=NMDS2,fill=Treatment,group=Treatment),alpha=0.30) + # add the convex hulls
  geom_text(data=species.scores,aes(x=NMDS1,y=NMDS2,label=species)) +  # add the species labels
  geom_point(data=data.scores,aes(x=NMDS1,y=NMDS2,shape=Treatment,colour=Treatment),size=4) + # add the point markers
  scale_colour_manual(values=c("A" = "darkorange", "P" = "dodgerblue2")) +
  theme_bw() + 
  theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        axis.title.x = element_text(size=18), # remove x-axis labels
        axis.title.y = element_text(size=18), # remove y-axis labels
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())
#

#################################################
#                Cinnamon Mountain - ELEVATION = 3460 m
#################################################
cin<-subset(community, Site=="Cin_obs")
names(cin)
cin2<-cin[,c("ACHMIL","AGOGLA","Annual","ARCUVA", "ARNMOL","BOESTR","CARALB","COLLIN","DRAAUR","ERECON","ERIGLA","ERYGRA","FRAVIR","JUNDRU","LIGPOR", "POAALP","POALEP","POTGRA","PSEMON","SENCRA","CLATON","VIOADU")]
cin3<-cin2[complete.cases(cin2), ]
colnames(cin3)[colnames(cin3)=="Annual"] <- "ANDSEP"

#bray-curtis
cinperm<-adonis(cin2 ~ Treatment*Pairing  ,data=cin, permutations=99)
# jaccards
cin2J<- cin2
cin2J[cin2J>0] <-1
modcin2J<-adonis(cin2J ~ Elevation*Treatment*Pairing,data=cin, permutations=99)



cin3dis<-vegdist(cin3, method="bray")
cin
cin$Treatment<-c("P", "P","P", "P","P","P", "P","P", "P","P","A","A","A","A","A","A","A","A","A","A")
mod<-betadisper(cin3dis, group=cin$Treatment)
anova(mod)

cin.mds <- metaMDS(cin3, trymax = 1000)
data.scores <- as.data.frame(scores(cin.mds))  #Using the scores function from vegan to extract the site scores and convert to a data.frame
data.scores$site <- rownames(data.scores)  # create a column of site names, from the rownames of data.scores
data.scores$Treatment <- Treatment  #  add the Treatment variable created earlier
head(data.scores)  #look at the data

species.scores <- as.data.frame(scores(cin.mds, "species"))  #Using the scores function from vegan to extract the species scores and convert to a data.frame
species.scores$species <- rownames(species.scores)  # create a column of species, from the rownames of species.scores
head(species.scores)  #look at the data

Treatment.a <- data.scores[data.scores$Treatment == "A", ][chull(data.scores[data.scores$Treatment == 
                                                                       "A", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment A
Treatment.p <- data.scores[data.scores$Treatment == "P", ][chull(data.scores[data.scores$Treatment == 
                                                                       "P", c("NMDS1", "NMDS2")]), ]  # hull values for Treatment B

hull.data <- rbind(Treatment.a, Treatment.p)  #combine Treatment.a and Treatment.b
hull.data

Cinnie.MDS<-ggplot() + 
  geom_polygon(data=hull.data,aes(x=NMDS1,y=NMDS2,fill=Treatment,group=Treatment),alpha=0.30) + # add the convex hulls
  geom_text(data=species.scores,aes(x=NMDS1,y=NMDS2,label=species)) +  # add the species labels
  geom_point(data=data.scores,aes(x=NMDS1,y=NMDS2,shape=Treatment,colour=Treatment),size=4) + # add the point markers
  scale_colour_manual(values=c("A" = "darkorange", "P" = "dodgerblue2")) +
  theme_bw() + 
  theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        axis.title.x = element_text(size=18), # remove x-axis labels
        axis.title.y = element_text(size=18), # remove y-axis labels
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())

######################################################################################################
###  FUNGAL COLONIZATION DATA
######################################################################################################
roots<-read.csv("McKibben_Henning_2018_colonization_data.csv")
roots

grouped <- group_by(roots, Elevation)
summarise(grouped, meanMyco=mean(Myco), meanDSE=mean(DSE))
(84.8-67.8)/84.8

## Model selection using lme
lme.Myco <- lme(Myco~Treatment*Elevation, random=~1|Pairing, data=roots,  method="ML")
lme.Myco1 <- lme(Myco~Treatment, random=~1|Pairing, data=roots,  method="ML")
lme.Myco2 <- lme(Myco~Elevation, random=~1|Pairing, data=roots,  method="ML")
lme.Myco3 <- lme(Myco~1, random=~1|Pairing, data=roots,  method="ML")
anova(lme.Myco,lme.Myco1,lme.Myco2,lme.Myco3)
# best fit model
lme.Myco <- lme(Myco~Treatment, random=~1|Pairing, data=roots,  method="REML")
summary(lme.Myco)
Anova(lme.Myco)

## Model selection using lme
lme.DSE <- lme(DSE~Treatment*Elevation, random=~1|Pairing, data=roots,  method="ML")
lme.DSE1 <- lme(DSE~Treatment, random=~1|Pairing, data=roots,  method="ML")
lme.DSE2 <- lme(DSE~Elevation, random=~1|Pairing, data=roots,  method="ML")
lme.DSE3 <- lme(DSE~1, random=~1|Pairing, data=roots,  method="ML")
anova(lme.DSE,lme.DSE1,lme.DSE2,lme.DSE3)
# best fit model
lme.DSE <- lme(DSE~Elevation, random=~1|Pairing, data=roots,  method="REML")
summary(lme.DSE)
Anova(lme.DSE)


AMF<-ggplot(data=roots, aes(x=factor(Elevation), y=Myco, fill = Treatment))
AMFfig<- AMF + geom_boxplot() + ylim(0,100) + theme_bw(base_size=14)+ scale_fill_manual(values=c( "darkorange","dodgerblue2"))+ theme(legend.position="none")

DSE<-ggplot(data=roots, aes(x=factor(Elevation), y=DSE, fill = Treatment))
DSEfig<-DSE + geom_boxplot()+ ylim(0,100) +  theme_bw(base_size=14)+ scale_fill_manual(values=c("darkorange", "dodgerblue2"))+ theme(legend.position="none")

p1 <- ggplot_gtable(ggplot_build(AMFfig))
p2 <- ggplot_gtable(ggplot_build(DSEfig))

maxWidth = unit.pmax(p1$widths[2:3], p2$widths[2:3])

p1$widths[2:3] <- maxWidth
p2$widths[2:3] <- maxWidth

tiff("Henning_McKibben_Figure3.tiff", width = 4, height = 6, units = 'in', res = 300, compression = 'none')
grid.arrange(p1,p2, ncol = 1)
dev.off()