# Xenopus microbiome development biodiversity analysis
# Comparaison des communautés entre tetards prométamorphiques, prémétamorphiques, métamorphiques et adultes.
# Chargement des packages
library(picante)
## Loading required package: ape
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
## Loading required package: nlme
library(ggplot2)
library(phytools)
## Loading required package: maps
##
## Attaching package: 'phytools'
## The following object is masked from 'package:vegan':
##
## scores
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
##
## collapse
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(grid)
library(PMCMRplus)
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
library(multcompView)
## Folders, Themes, colors
source("prelude.R")
# Setting up directories
data_dir <-paste0(data_dir_path,"xpdev")
output_dir <- paste0(output_dir_path,"Supp_Figure_4")
# Working on abundance data from 100 rarefactions
setwd(data_dir)
allcom<-read.table("xpdev_fichier_allcom.txt", header=TRUE, row.names=1)
# Getting rid of OTUs with zero abundance
allcom<-allcom %>% select_if((function(col) is.numeric(col) && sum(col) > 0))
commRS<-decostand(allcom, method="total")
comm<-commRS
#Importing metadata
metadata<-read.csv("nmt3_2_metadata_xp_dev.csv",sep=";",header=TRUE,row.names=1)
metadata$LifeStage<-as.factor(metadata$LifeStage)
#Importing phylogeny
phy<-read.tree("nmt3_sequences_mpr.tree")
#Rooting the tree using midpoint rooting
phy<-midpoint.root(phy)
#Getting rid of taxa that are missing in the studied samples
phy<-prune.sample(comm, phy)
#Combining the input data
combined <- match.phylo.comm(phy, comm)
phy <- combined$phy
comm <- combined$comm
metadata <- metadata[rownames(comm), ]
#Checking the index
all.equal(rownames(comm), rownames(metadata))
## [1] TRUE
#Checking the normalisation
rowSums(comm)
## XT1A2G XT1A4G XT1A5G XT1B1G XT1B2G XT1B2Gbis XT1B4G XT1B5G
## 1 1 1 1 1 1 1 1
## XTLS1 XTLS10 XTLS11 XTLS12 XTLS13 XTLS14 XTLS15 XTLS16
## 1 1 1 1 1 1 1 1
## XTLS17 XTLS2 XTLS3 XTLS4 XTLS5 XTLS6 XTLS7 XTLS8
## 1 1 1 1 1 1 1 1
## XTLS9 XTLSM1 XTLSM2 XTLSM3 XTLTW2M XTLTW3F XTLTW3M XTLTW4F
## 1 1 1 1 1 1 1 1
## XTLTW5F
## 1
#######################################
# Analysis of taxonomic diversity #
#######################################
# Defining the output dir
setwd(output_dir)
# For the boxplot with ggplot2
# OTU richness
df<-data.frame(Nb.OTU= specnumber(comm),Sample=metadata$LifeStage)
#Phylogenetic diversity
comm.pd <- pd(comm, phy)
df2<-data.frame(Faith.PD= comm.pd$PD,Sample=metadata$LifeStage)
#Significance
kw_sr<-kruskal.test(specnumber(comm) ~ metadata$LifeStage, data= comm)
kw_pd<-kruskal.test(comm.pd$PD ~ metadata$LifeStage, data= comm)
#Pairwise comparisons using Dunn-Bonferroni test implemented in PMCMRplus
pp_sr<-kwAllPairsDunnTest(x=specnumber(comm), metadata$LifeStage,p.adjust.method="bonferroni")
## Warning in kwAllPairsDunnTest.default(x = specnumber(comm),
## metadata$LifeStage, : Ties are present. z-quantiles were corrected for ties.
pp_pd<-kwAllPairsDunnTest(x=comm.pd$PD, metadata$LifeStage,p.adjust.method="bonferroni")
# Including results of significance to plots
mymat_sr <-tri.to.squ(pp_sr$p.value)
mymat_pd <-tri.to.squ(pp_pd$p.value)
# Assigning letter code for significance
myletters_sr <- multcompLetters(mymat_sr,compare="<=",threshold=0.05,Letters=letters)
myletters_pd <- multcompLetters(mymat_pd,compare="<=",threshold=0.05,Letters=letters)
# Put the letters in a dataframe
myletters_sr_df <- data.frame(Sample=names(myletters_sr$Letters),letter = myletters_sr$Letters )
myletters_pd_df <- data.frame(Sample=names(myletters_pd$Letters),letter = myletters_pd$Letters )
ggplot(df,aes(x=Sample,y=Nb.OTU))+ geom_boxplot(outlier.colour=NA)+ geom_dotplot(binaxis="y",binwidth=.5,stackdir="center",dotsize=15,fill="grey")+geom_text(data = myletters_sr_df, aes(label = letter, y = 365 ),fontface="bold",size=6) + scale_x_discrete(limits=c("Premetamorph","Prometamorph","Metamorph","Froglet","Adult"), labels=c("Premeta.","Prometa.","Meta.","Froglet","Adult"))+theme_nb()
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
ggsave("Plot_stage_1_2.pdf",width=26,height=20,units="cm")
ggplot(df2,aes(x=Sample,y=Faith.PD))+geom_boxplot(outlier.colour=NA)+geom_dotplot(binaxis="y",binwidth=.5,stackdir="center",dotsize=0.75,fill="grey")+geom_text(data = myletters_pd_df, aes(label = letter, y = 19.5 ),fontface="bold",size=6)+scale_x_discrete(limits=c("Premetamorph","Prometamorph","Metamorph","Froglet","Adult"), labels=c("Premeta.","Prometa.","Meta.","Froglet","Adult"))+theme_nb()
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
ggsave("Plot_stage_2_2.pdf",width=26,height=20,units="cm")
#Diversity using the Shannon index
comm.H<-diversity(comm, index="shannon", base=exp(1))
dfH<-data.frame(H= comm.H,Sample=metadata$LifeStage)
plotH<-ggplot(dfH,aes(x=Sample,y=H))+scale_x_discrete(limits=c("Premetamorph","Prometamorph","Metamorph","Froglet","Adult"), labels=c("Premeta.","Prometa.","Meta.","Froglet","Adult"))+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey")+geom_boxplot(outlier.colour=NA,alpha=0.1)+theme_nb()+labs(y="Shannon diversity",x="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plotH
#Diversity using the Simpson index
comm.D<-diversity(comm, index="simpson")
dfD<-data.frame(D= comm.D,Sample=metadata$LifeStage)
plotD<-ggplot(dfD,aes(x=Sample,y=D))+scale_x_discrete(limits=c("Premetamorph","Prometamorph","Metamorph","Froglet","Adult"), labels=c("Premeta.","Prometa.","Meta.","Froglet","Adult"))+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey")+geom_boxplot(outlier.colour=NA,alpha=0.1)+theme_nb()+labs(y="Simpson diversity",x="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plotD
#Diversity using the inverse Simpson index
comm.I<-diversity(comm, index="invsimpson")
dfI<-data.frame(I= comm.I,Sample=metadata$LifeStage)
plotI<-ggplot(dfI,aes(x=Sample,y=I))+scale_x_discrete(limits=c("Premetamorph","Prometamorph","Metamorph","Froglet","Adult"), labels=c("Premeta.","Prometa.","Meta.","Froglet","Adult"))+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey")+geom_boxplot(outlier.colour=NA,alpha=0.1)+theme_nb()+labs(y="Inverse Simpson",x="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plotI
#Diversity using the equitability index
comm.J <- comm.H/log(specnumber(comm), base=exp(1))
dfJ<-data.frame(J= comm.J,Sample=metadata$LifeStage)
plotJ<-ggplot(dfJ,aes(x=Sample,y=J))+scale_x_discrete(limits=c("Premetamorph","Prometamorph","Metamorph","Froglet","Adult"), labels=c("Premeta.","Prometa.","Meta.","Froglet","Adult"))+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey")+geom_boxplot(outlier.colour=NA,alpha=0.1)+theme_nb()+labs(y="Pielou evenness",x="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plotJ
grid.arrange(plotH, plotD, plotI, plotJ,nrow=4,ncol=1)
#pdf("plotHDIJ_Stage_2.pdf",width=10,height=10)
#require(gridExtra)
#print(grid.arrange(plotH, plotD, plotI, plotJ,nrow=4,ncol=1))
#dev.off()
#Comparison of diversity indices
rankindex(comm,comm,c("euc","man","bray","jac","kul"))
## euc man bray jac kul
## -0.21234567 0.08143398 0.08143398 0.08143398 0.08143398
# Computing dissimilarity using Bray-Curtis distance (weighted)
comm.bc.dist <- vegdist(comm, method = "bray")
# Clustering to analyse samples - Bray
comm.bc.clust <- hclust(comm.bc.dist, method = "average")
pdf("plot_bc_clust_stage.pdf",width=8,height=6)
print(plot(comm.bc.clust, ylab = "Bray-Curtis dissimilarity"))
## NULL
dev.off()
## quartz_off_screen
## 2
#Dissimilarity using Sorenson (unweighted)
comm.so.dist<-vegdist(comm,method="bray",binary="TRUE")
# Clustering to analyse samples - Sorenson
comm.so.clust <- hclust(comm.so.dist, method = "average")
pdf("plot_so_clust_stage.pdf",width=8,height=6)
print(plot(comm.so.clust, ylab = "Sorenson dissimilarity"))
## NULL
dev.off()
## quartz_off_screen
## 2
################################
# Ordination NMDS Bray-Curtis #
################################
# Ordination using NMDS
comm.bc.mds <- metaMDS(comm, dist = "bray")
## Run 0 stress 0.1690183
## Run 1 stress 0.1836237
## Run 2 stress 0.1691755
## ... Procrustes: rmse 0.005939721 max resid 0.02358807
## Run 3 stress 0.2204997
## Run 4 stress 0.3857212
## Run 5 stress 0.1690453
## ... Procrustes: rmse 0.00196101 max resid 0.008065467
## ... Similar to previous best
## Run 6 stress 0.203513
## Run 7 stress 0.1691752
## ... Procrustes: rmse 0.005822825 max resid 0.02318964
## Run 8 stress 0.1981167
## Run 9 stress 0.1690183
## ... Procrustes: rmse 2.266568e-05 max resid 6.632704e-05
## ... Similar to previous best
## Run 10 stress 0.3928898
## Run 11 stress 0.1935153
## Run 12 stress 0.1691974
## ... Procrustes: rmse 0.006055161 max resid 0.0229117
## Run 13 stress 0.2280351
## Run 14 stress 0.1882471
## Run 15 stress 0.2142176
## Run 16 stress 0.1836299
## Run 17 stress 0.202762
## Run 18 stress 0.2005621
## Run 19 stress 0.194421
## Run 20 stress 0.1691985
## ... Procrustes: rmse 0.006279185 max resid 0.02374309
## *** Solution reached
# Checking the NMDS
stressplot(comm.bc.mds)
#pdf("plot_nmds_stressplot.pdf",width=8,height=6)
#print(stressplot(comm.bc.mds))
#dev.off()
# Plotting NMDS ordination using Bray-Curtis distance
mds.fig2 <- ordiplot(comm.bc.mds, type = "none",xlim=c(-2.5,3.5),ylim=c(-2,2))
points(mds.fig2, "species", pch = 1,cex=0.5, col = "grey")
ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Premetamorph",col="burlywood")
points(mds.fig2, "sites", pch = 19, col = "burlywood", select = metadata$LifeStage == "Premetamorph")
ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Prometamorph",col="aquamarine4")
points(mds.fig2, "sites", pch = 19, col = "aquamarine4", select = metadata$LifeStage == "Prometamorph")
ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Metamorph",col="green")
points(mds.fig2, "sites", pch = 19, col = "green", select = metadata$LifeStage == "Metamorph")
ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Froglet",col="gold")
points(mds.fig2, "sites", pch = 19, col = "gold", select = metadata$LifeStage == "Froglet")
ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Adult",col="darkorange")
points(mds.fig2, "sites", pch = 19, col = "darkorange", select = metadata$LifeStage == "Adult")
ordiellipse(comm.bc.mds, metadata$LifeStage, conf = 0.95, label = FALSE,lty = 'dotted')
#pdf("plot3_nmds_2.pdf",width=8,height=6)
#print(mds.fig2 <- ordiplot(comm.bc.mds, type = "none",xlim=c(-2.5,3.5),ylim=c(-2,2)))
#print(points(mds.fig2, "species", pch = 1,cex=0.5, col = "grey"))
#print(ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Premetamorph",col="burlywood"))
#print(points(mds.fig2, "sites", pch = 19, col = "burlywood", select = metadata$LifeStage == "Premetamorph"))
#print(ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Prometamorph",col="aquamarine4"))
#print(points(mds.fig2, "sites", pch = 19, col = "aquamarine4", select = metadata$LifeStage == "Prometamorph"))
#print(ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Metamorph",col="green"))
#print(points(mds.fig2, "sites", pch = 19, col = "green", select = metadata$LifeStage == "Metamorph"))
#print(ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Froglet",col="gold"))
#print(points(mds.fig2, "sites", pch = 19, col = "gold", select = metadata$LifeStage == "Froglet"))
#print(ordispider(comm.bc.mds,groups=metadata$LifeStage,show.groups="Adult",col="darkorange"))
#print(points(mds.fig2, "sites", pch = 19, col = "darkorange", select = metadata$LifeStage == "Adult"))
#print(ordiellipse(comm.bc.mds, metadata$LifeStage, conf = 0.95, label = FALSE,lty = 'dotted'))
#dev.off()
################################
# Ordination beta-dispersion #
################################
# Plotting the ordination results of beta-dispersion PCoA (PERMDISP)
comm.betadisp=betadisper(comm.bc.dist, metadata$LifeStage, type="median")
comm.betadisp
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = comm.bc.dist, group = metadata$LifeStage, type =
## "median")
##
## No. of Positive Eigenvalues: 26
## No. of Negative Eigenvalues: 6
##
## Average distance to median:
## Adult Froglet Metamorph Premetamorph Prometamorph
## 0.3611 0.5265 0.4041 0.2061 0.2938
##
## Eigenvalues for PCoA axes:
## (Showing 8 of 32 eigenvalues)
## PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
## 3.2375 2.0213 1.2734 1.0513 0.7138 0.5901 0.4745 0.3849
plot(comm.betadisp)
#pdf("plot1_betadispa_2.pdf",width=8,height=6)
#print(mds.fig3a<-plot(comm.betadisp))
#dev.off()
# Another plot of PCoA for beta dispersion
#pdf("plot1_betadisp.pdf",width=8,height=6)
#print(mds.fig3 <- ordiplot(comm.betadisp, type = "none"))
#print(ordispider(comm.betadisp,groups=metadata$LifeStage,show.groups="Premetamorph",col="burlywood"))
#print(points(mds.fig3, "sites", pch = 19, col = "burlywood", select = metadata$LifeStage == "Premetamorph"))
#print(ordispider(comm.betadisp,groups=metadata$LifeStage,show.groups="Prometamorph",col="aquamarine4"))
#print(points(mds.fig3, "sites", pch = 19, col = "aquamarine4", select = metadata$LifeStage == "Prometamorph"))
#print(ordispider(comm.betadisp,groups=metadata$LifeStage,show.groups="Metamorph",col="green"))
#print(points(mds.fig3, "sites", pch = 19, col = "green", select = metadata$LifeStage == "Metamorph"))
#print(ordispider(comm.betadisp,groups=metadata$LifeStage,show.groups="Froglet",col="gold"))
#print(points(mds.fig3, "sites", pch = 19, col = "gold", select = metadata$LifeStage == "Froglet"))
#print(ordispider(comm.betadisp,groups=metadata$LifeStage,show.groups="Adult",col="darkorange"))
#print(points(mds.fig3, "sites", pch = 19, col = "darkorange", select = metadata$LifeStage == "Adult"))
#print(ordiellipse(comm.betadisp, metadata$LifeStage, conf = 0.95, label = FALSE,lty = 'dotted'))
#dev.off()
#######################################################
# Analysis of community structure using MDP and MNTD #
#######################################################
#Phylogenetic distance to measure the community phylogenetic structure without abundance
phy.dist <- cophenetic(phy)
comm.sesmpd <- ses.mpd(comm, phy.dist, null.model = "richness", abundance.weighted = FALSE, runs = 999)
comm.sesmntd <- ses.mntd(comm, phy.dist, null.model = "richness", abundance.weighted = FALSE, runs = 999)
#Inspecting obtained values
#comm.sesmpd
#comm.sesmntd
#Significance
kw_sesmpd<-kruskal.test(comm.sesmpd$mpd.obs.z ~ metadata$LifeStage, data= comm)
kw_sesmntd<-kruskal.test(comm.sesmntd$mntd.obs.z ~ metadata$LifeStage, data= comm)
#Pairwise comparisons using Dunn-Bonferroni test implemented in PMCMRplus
pp_sesmpd<-kwAllPairsDunnTest(x=comm.sesmpd$mpd.obs.z, metadata$LifeStage,p.adjust.method="bonferroni")
pp_sesmntd<-kwAllPairsDunnTest(x=comm.sesmntd$mntd.obs.z, metadata$LifeStage,p.adjust.method="bonferroni")
# Including results of significance to plots
mymat_sesmpd <-tri.to.squ(pp_sesmpd$p.value)
mymat_sesmntd <-tri.to.squ(pp_sesmntd$p.value)
# Assigning letter code for significance
myletters_sesmpd <- multcompLetters(mymat_sesmpd,compare="<=",threshold=0.05,Letters=letters)
myletters_sesmntd <- multcompLetters(mymat_sesmntd,compare="<=",threshold=0.05,Letters=letters)
# Put the letters in a dataframe
myletters_sesmpd_df <- data.frame(Sample=names(myletters_sesmpd$Letters),letter = myletters_sesmpd$Letters )
myletters_sesmntd_df <- data.frame(Sample=names(myletters_sesmntd$Letters),letter = myletters_sesmntd$Letters )
#Preparing data for plotting SES(MPD) et SES(MNTD)
df3<-data.frame(SES_MPD = comm.sesmpd$mpd.obs.z, Sample = metadata$LifeStage)
df4<-data.frame(SES_MNTD = comm.sesmntd$mntd.obs.z, Sample = metadata$LifeStage)
plot3<-ggplot(df3,aes(x=Sample,y=SES_MPD))+geom_boxplot()+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey",alpha=0.9)+scale_x_discrete(limits=c("Premetamorph","Prometamorph","Metamorph","Froglet","Adult"), labels=c("Premeta.","Prometa.","Meta.","Froglet","Adult"))+labs(x="",y="SES MPD")+geom_text(data = myletters_sesmpd_df, aes(label = letter, y = 0 ), colour="black", size=5,fontface="bold")+theme_nb()
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plot4<-ggplot(df4,aes(x=Sample,y=SES_MNTD))+geom_boxplot()+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey",alpha=0.9)+scale_x_discrete(limits=c("Premetamorph","Prometamorph","Metamorph","Froglet","Adult"), labels=c("Premeta.","Prometa.","Meta.","Froglet","Adult"))+labs(x="",y="SES MPD")+geom_text(data = myletters_sesmntd_df, aes(label = letter, y = 0 ), colour="black", size=5,fontface="bold")+theme_nb()
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
#Plot for Supp Figure
grid.arrange(plot3, plot4, ncol=2)
#require(gridExtra)
#pdf("plot_sesmpdmntd_12_stage.pdf",width=20,height=8)
#print(grid.arrange(plot3, plot4, ncol=2))
#dev.off()
#####################################
# Beta-diversité phylogénétique #
#####################################
# beta diversity: Taxonomic (Bray-Curtis) dissimilarity explained (Permanova)
adonis(comm.bc.dist ~ LifeStage, data = metadata)
##
## Call:
## adonis(formula = comm.bc.dist ~ LifeStage, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## LifeStage 4 6.6035 1.65089 9.5018 0.5758 0.001 ***
## Residuals 28 4.8648 0.17374 0.4242
## Total 32 11.4684 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Significance of beta-dispersion
comm.betadisp.perm=permutest(comm.betadisp, group= metadata$LifeStage, type="median", permutations=999, pairwise=TRUE)
comm.betadisp.perm
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 4 0.38703 0.096758 3.1512 999 0.021 *
## Residuals 28 0.85975 0.030705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Adult Froglet Metamorph Premetamorph Prometamorph
## Adult 1.0000e-03 7.7500e-01 1.0000e-03 0.074
## Froglet 5.7799e-05 4.9900e-01 1.0000e-03 0.002
## Metamorph 7.7731e-01 5.0130e-01 1.8500e-01 0.565
## Premetamorph 7.6427e-04 9.1688e-06 2.2212e-01 0.090
## Prometamorph 7.6114e-02 2.8073e-04 5.4855e-01 9.6521e-02
comm.betadisp.HSD<-TukeyHSD(comm.betadisp)
plot(comm.betadisp.HSD)
#Test for variance homogeneity
anova(comm.betadisp)
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 4 0.38703 0.096758 3.1512 0.02939 *
## Residuals 28 0.85975 0.030705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# beta diversity: compute phylogenetic dissimilarity using MPD
comm.mpd.dist<-comdist(comm,phy.dist,abundance.weighted=TRUE)
adonis(comm.mpd.dist ~ LifeStage, data = metadata)
##
## Call:
## adonis(formula = comm.mpd.dist ~ LifeStage, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## LifeStage 4 1.2278 0.306947 3.6682 0.34384 0.001 ***
## Residuals 28 2.3430 0.083678 0.65616
## Total 32 3.5708 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# beta diversity: compute phylogenetic dissimilarity using MNTD
comm.mntd.dist <- comdistnt(comm, phy.dist, abundance.weighted = TRUE)
adonis(comm.mntd.dist ~ LifeStage, data = metadata)
##
## Call:
## adonis(formula = comm.mntd.dist ~ LifeStage, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## LifeStage 4 0.0140614 0.0035153 90.649 0.92831 0.001 ***
## Residuals 28 0.0010858 0.0000388 0.07169
## Total 32 0.0151472 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# calculate Mantel correlation for taxonomic Bray-Curtis vs. phylogenetic
# MPD diversity
mantel(comm.bc.dist, comm.mpd.dist)
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = comm.bc.dist, ydis = comm.mpd.dist)
##
## Mantel statistic r: 0.6143
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.082 0.101 0.127 0.144
## Permutation: free
## Number of permutations: 999
#MNTD diversity
mantel(comm.bc.dist, comm.mntd.dist)
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = comm.bc.dist, ydis = comm.mntd.dist)
##
## Mantel statistic r: 0.4593
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0667 0.0846 0.0987 0.1151
## Permutation: free
## Number of permutations: 999