Comparable early-stage decomposition but contrasting underlying drivers between surface and cave habitats along an elevational gradient

Analyses of 16S microbiome dataset

Importing the data and creating a phyloseq object

# load following packages
library(vegan)

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.6-2

library(phyloseq) 
library(microbiome)

## Loading required package: ggplot2

## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2020 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>

## 
## Attaching package: 'microbiome'

## The following object is masked from 'package:ggplot2':
## 
##     alpha

## The following object is masked from 'package:vegan':
## 
##     diversity

## The following object is masked from 'package:base':
## 
##     transform

library(ape)
library(knitr)
library(cluster)
library(viridis)

## Loading required package: viridisLite

library(ggplot2)
library(grid)
library(gridExtra)
library(picante)

## Loading required package: nlme

library(adephylo)

## Loading required package: ade4

library(igraph)

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:ape':
## 
##     degree, edges, mst, ring

## The following object is masked from 'package:microbiome':
## 
##     diversity

## The following object is masked from 'package:vegan':
## 
##     diversity

## The following object is masked from 'package:permute':
## 
##     permute

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:igraph':
## 
##     as_data_frame, groups, union

## The following object is masked from 'package:nlme':
## 
##     collapse

## The following object is masked from 'package:gridExtra':
## 
##     combine

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(microeco)
library(file2meco)
library(pairwiseAdonis)

# Bring in the ASV table and transform it 
abund_table <- read.csv2("ASV_table_Tenerife.csv", row.names = 1, check.names = FALSE)
#Transpose the OTU table
abund_table <- t(abund_table)
#Bring in the Taxanomy table
taxa <- read.csv2("Tax_table_Tenerife.csv", row.names = 1, check.names = FALSE)
#Bring in the metadata table
meta_table <- read.csv2("Meta_table_Tenerife.csv", row.names = 1, check.names = FALSE)
meta_table <- meta_table[rownames(abund_table),]
#bring in the bacterial phylogeny
OTU_tree <- read_tree("Tenerife.tre")

#prepare the tables for phyloseq
OTU1 <- otu_table(as.matrix(abund_table), taxa_are_rows = FALSE)
TAX1 <- tax_table(as.matrix(taxa))
SAM <- sample_data(meta_table)
OTU_tree <- compute.brlen(OTU_tree, method = "Grafen")

#create a phyloseq object with all the data tables
physeq <- merge_phyloseq(phyloseq(OTU1, TAX1), SAM, OTU_tree)

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

#Different samples in the dataset are sequenced in different depths so we rarefied our phyloseq object. Here we used the sequence number in the sample with the lowest sequencing depth.

physeq_R<- rarefy_even_depth(physeq, sample.size = min(sample_sums(physeq)),
                  rngseed = TRUE, replace = TRUE, trimOTUs = TRUE, verbose = TRUE)

## `set.seed(TRUE)` was used to initialize repeatable random subsampling.

## Please record this for your records so others can reproduce.

## Try `set.seed(TRUE); .Random.seed` for the full vector

## ...

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## 1993OTUs were removed because they are no longer 
## present in any sample after random subsampling

## ...

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

Now we can calculate alpha diversities

#if you use "all" as the index you will get a table with all the alpha diversity indexes that microbiome package calculate. But here we specifed jut to get chao1 richness estimate and the Shannon's diversity index.
Div <- alpha(physeq_R, index = c("observed", "chao1", "diversity_shannon"))

## Observed richness

## Other forms of richness

## Diversity

## Evenness

## Dominance

## Rarity

write.csv2(Div, file = "Alpha_tenerife.csv")

We calculated the Faith’s phylogenetic diversity seperately using Picante package

asvtab <- as.data.frame(physeq_R@otu_table)
tree <- physeq_R@phy_tree
df.pd <- pd(asvtab, tree,include.root=T)
write.csv2(df.pd, file = "PD_tenerife.csv")

Plotting phylum level bargraphs

#prune taxa that are less than 2%
prune.dat <- prune_taxa(taxa_sums(physeq_R) > 2, physeq_R)

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

#get dominanat phyla
top6 <- names(sort(taxa_sums(prune.dat), decreasing=TRUE)[1:200])
dat.aglo <- tax_glom(physeq_R, taxrank = "Phylum")

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

dat.trans = transform_sample_counts(dat.aglo, function(x) x/sum(x))

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

#all the samples
prune.dat.two = prune_taxa(taxa_sums(dat.trans) > 0.001, dat.trans)

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'RNeXML'

#only top taxa
prune.dat.two = prune_taxa(top6, dat.trans)

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'RNeXML'

dat.dataframe = psmelt(prune.dat.two)
dat.agr = aggregate(Abundance~Cave+Habitat+Phylum, data=dat.dataframe, FUN=mean)

#plot the bar graph. Here we also faucet the graph according to Surface and Cave
ggplot(dat.agr, aes(x=Cave, y=Abundance, fill=Phylum)) + geom_bar(stat="identity") + facet_grid(~Habitat, scale="free")+ theme_classic() + scale_fill_manual(values=c("#a2a2b0","#7979a6","#111112", "#5b5bc9","#1616c9", "#1f2cd1", "#ed47bb", "#818282", "#ed4945", "#0c6e6e", "#a30b07", "#bd921c", "#0af20a", "#038703", "#636b6b"))

Beta diversity of microbiomes

#This can be calculated using multiple distance matrixes but here we used Bray-Curtis and weighted Unifrac distances. In following examples we will only show analyses with Bray-Curtis distances. For visualizing we used NMDS plots.

ord <- ordinate(physeq_R, "NMDS", "bray")

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1647741 
## Run 1 stress 0.2026344 
## Run 2 stress 0.1523125 
## ... New best solution
## ... Procrustes: rmse 0.1126109  max resid 0.2661109 
## Run 3 stress 0.1523125 
## ... Procrustes: rmse 5.4421e-06  max resid 1.411551e-05 
## ... Similar to previous best
## Run 4 stress 0.1937815 
## Run 5 stress 0.1879122 
## Run 6 stress 0.1523141 
## ... Procrustes: rmse 0.0003679265  max resid 0.001774492 
## ... Similar to previous best
## Run 7 stress 0.1690842 
## Run 8 stress 0.1541766 
## Run 9 stress 0.1465648 
## ... New best solution
## ... Procrustes: rmse 0.08988216  max resid 0.4764822 
## Run 10 stress 0.1539707 
## Run 11 stress 0.1912181 
## Run 12 stress 0.1642052 
## Run 13 stress 0.1898784 
## Run 14 stress 0.146565 
## ... Procrustes: rmse 7.668128e-05  max resid 0.0003949434 
## ... Similar to previous best
## Run 15 stress 0.1539707 
## Run 16 stress 0.1713799 
## Run 17 stress 0.1465649 
## ... Procrustes: rmse 0.0001048992  max resid 0.000425372 
## ... Similar to previous best
## Run 18 stress 0.1465648 
## ... New best solution
## ... Procrustes: rmse 5.803446e-05  max resid 0.0003247318 
## ... Similar to previous best
## Run 19 stress 0.1717456 
## Run 20 stress 0.1539707 
## *** Solution reached

p2 <- plot_ordination(physeq_R, ord, type="samples",color="Cave",shape = "Habitat",  axes = c(1,2))
p2 + geom_point(size=7)+ scale_color_manual(values=c("#990000","#ed5e21", "#edb64e", "#a4eda1", "#f041c4", "#87dee8", "#1f60ed", "#042b80")) + theme_classic() +  scale_shape_manual(values=c(15, 17))+ stat_ellipse(aes(group = Habitat), linetype = 2)

We used the envfit function in Vegan package to investigate the influence of enviroenmental parameters on beta diversities

# We only chose Elevation, pH, water, SOM and DOC for these analyses. We also conducted the envfit analyses separately for two habitat types (caves and surfaces)
env <- physeq_R@sam_data[,4:8]
en <- envfit(ord, env, permutations = 999, na.rm = TRUE)
en

## 
## ***VECTORS
## 
##                  NMDS1    NMDS2     r2 Pr(>r)    
## Elevation      0.52560 -0.85073 0.5456  0.001 ***
## Water_content -0.73870 -0.67403 0.1453  0.034 *  
## pH            -0.81362  0.58140 0.7060  0.001 ***
## SOM            0.75350  0.65745 0.0880  0.156    
## DOC            0.11530  0.99333 0.0216  0.628    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

We will use adonis2 to conduct PERMANOVAs to see the microbial community level differences in soil zones (habitat) and between different caves.

#first calculate the similarity of communities. Here I use Bray-Curtis.
bDist <- distance(physeq_R, method = "bray")

#Then we used the distance matrix to investigate the influence of soil zone (Habitat) and the cave on microbial community structures.
AD <- adonis2(bDist~  physeq_R@sam_data$Cave + physeq_R@sam_data$Habitat, by = "margin", permutations = 10000)

AD

## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = bDist ~ physeq_R@sam_data$Cave + physeq_R@sam_data$Habitat, permutations = 10000, by = "margin")
##                           Df SumOfSqs      R2      F    Pr(>F)    
## physeq_R@sam_data$Cave     7   6.1672 0.30952 2.5553 9.999e-05 ***
## physeq_R@sam_data$Habitat  1   1.2852 0.06450 3.7276 9.999e-05 ***
## Residual                  36  12.4122 0.62295                     
## Total                     44  19.9248 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Finally we investigated the correlations between the relative abundance of bacterial genera with k and S seperately in the two habitat types.

# first we subsample the data based on the habitat
Sub<- subset_samples(physeq_R, Habitat=="Surface")

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'RNeXML'

#Then we  convert our phyloseq object to fit microeco package
data <- phyloseq2meco(Sub)

## 10918 taxa are removed from the otu_table, as the abundance is 0 ...

#calculate relative abundances
data$cal_abund()

## The row number of tax_table is not equal to that of otu_table ...

## Automatically applying tidy_dataset() function to trim the data ...

## microtable class:
## sample_table have 23 rows and 20 columns
## otu_table have 13260 rows and 23 columns
## tax_table have 13260 rows and 7 columns
## phylo_tree have 13260 tips

## The result is stored in object$taxa_abund ...

#calculate correlations
t1 <- trans_env$new(dataset = data, add_data = data$sample_table[, 19:20])
t1$cal_cor (use_data = "Genus", p_adjust_method = "fdr")

## The correlation result is stored in object$res_cor ...

t1$plot_cor()

#remove any taxa without significant associations
t1$plot_cor(filter_feature = c(""))

Analyses of alpha diversities and soil properties

###Enviroenmental variable analyses (full dataset)

Env_full<-read.csv2("Env_var.csv", header = TRUE, row.names = 1)

Following packages are needed for the analyses

library(dplyr)
library(ggpubr)

## 
## Attaching package: 'ggpubr'

## The following object is masked from 'package:ape':
## 
##     rotate

library(FSA)

## Registered S3 methods overwritten by 'FSA':
##   method       from
##   confint.boot car 
##   hist.boot    car

## ## FSA v0.9.3. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✔ tibble  3.1.7     ✔ purrr   0.3.4
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ microbiome::alpha()     masks ggplot2::alpha()
## ✖ tibble::as_data_frame() masks dplyr::as_data_frame(), igraph::as_data_frame()
## ✖ dplyr::collapse()       masks nlme::collapse()
## ✖ dplyr::combine()        masks gridExtra::combine()
## ✖ purrr::compose()        masks igraph::compose()
## ✖ tidyr::crossing()       masks igraph::crossing()
## ✖ dplyr::filter()         masks stats::filter()
## ✖ dplyr::groups()         masks igraph::groups()
## ✖ dplyr::lag()            masks stats::lag()
## ✖ purrr::simplify()       masks igraph::simplify()

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:FSA':
## 
##     bootCase

## The following object is masked from 'package:dplyr':
## 
##     recode

library("ggplot2")
library("gplots")

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

library("viridis")
library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

Then we examined the corlinearity between environmental variables to see what variables are significantly correlated with each other using ggpairs function in ggally package. Make sure to specify the columns of the variables. NOTE: We did not include conductivity in this.

cor <- ggpairs(Env_full, columns = 5:18)
cor

Based on the correlation plots we can see that many environmental variables are correlated with each other. pH, Water, SOM and DOC do not correlate with one another so we use those variables for downstream analyses.

Before running analyses we investigated the distribution of data. Example is given with pH and DOC.

hist(Env_full$pH)

hist(Env_full$DOC)

If the distribution is normal (like in pH) then we ran a GLM with Gaussian distribution where link function set to identity and if the distribution is exponential (like in DOC) we change the link function to log. NOTE: “Habitat” indicte the soil zone (surface or within cave).

model1 <- glm(pH~ Cave + Habitat + Cave*Habitat, family = gaussian(link = "identity"), data = Env_full)

Anova(model1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: pH
##              LR Chisq Df Pr(>Chisq)    
## Cave          137.645  7  < 2.2e-16 ***
## Habitat        32.051  1  1.502e-08 ***
## Cave:Habitat   44.653  7  1.597e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

NOTE: We need to log transform DOC, but there are two soil samples where DOC is 0. In that case log transformation turn this into infinity and model give an error. So we will add a 0.01 to all DOC values to fix this issue. This wont have a big problem for our statistical outcome.

model2 <- glm(DOC+0.01~ Cave + Habitat + Cave*Habitat, family = gaussian(link = "log"), data = Env_full)

Anova(model2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: DOC + 0.01
##              LR Chisq Df Pr(>Chisq)    
## Cave           12.287  7    0.09151 .  
## Habitat         0.723  1    0.39531    
## Cave:Habitat   32.252  7  3.648e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

now lets plot pH and DOC

#pH
ggplot(Env_full, aes(x=Cave, y=pH)) + geom_boxplot (aes(color = Habitat, fill = Habitat), outlier.size = 3) + scale_color_manual(values=c("#140303","#220bb8"))+ scale_fill_manual(values=c("#474545","#6754e8")) +theme_classic()

#DOC
ggplot(Env_full, aes(x=Cave, y=DOC)) + geom_boxplot (aes(color = Habitat, fill = Habitat), outlier.size = 3) + scale_color_manual(values=c("#140303","#220bb8"))+ scale_fill_manual(values=c("#474545","#6754e8")) +theme_classic()

Enviroenmental variable associations with elevation and decomposition parameters (k and S)

For these analyses we will use separate datasheets for Surface and Cave

Cave_Env <- read.csv2("Cave_Env.csv", header = TRUE, row.names = 1)
Surf_Env <- read.csv2("Surf_Env.csv", header = TRUE, row.names = 1)

Following codes will focus on Surface environmental variables but the steps are the same. We conducted both linear regressions and polynomial regressions and then compare the models. Example is on pH.

#Generate the squire of the independent variable
EL2 <- Surf_Env$Elevation*Surf_Env$Elevation

#linear regression
model1 <- lm(pH~ Elevation, data = Surf_Env)

#polynomial regression
model2 <- lm(pH~ Elevation + EL2, data = Surf_Env)

#get model summaries
summary(model1)

## 
## Call:
## lm(formula = pH ~ Elevation, data = Surf_Env)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6931 -0.2703  0.2106  0.3677  1.1818 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.9495360  0.2388016  29.102  < 2e-16 ***
## Elevation   -0.0005513  0.0001636  -3.369  0.00277 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6681 on 22 degrees of freedom
## Multiple R-squared:  0.3403, Adjusted R-squared:  0.3103 
## F-statistic: 11.35 on 1 and 22 DF,  p-value: 0.00277

summary(model2)

## 
## Call:
## lm(formula = pH ~ Elevation + EL2, data = Surf_Env)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3018 -0.1537 -0.0046  0.1792  0.9764 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.447e+00  2.429e-01  26.542   <2e-16 ***
## Elevation    1.181e-03  5.173e-04   2.284   0.0329 *  
## EL2         -7.387e-07  2.131e-07  -3.467   0.0023 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5454 on 21 degrees of freedom
## Multiple R-squared:  0.5805, Adjusted R-squared:  0.5405 
## F-statistic: 14.53 on 2 and 21 DF,  p-value: 0.0001094

#get AIC values
AIC(model1)

## [1] 52.66426

AIC(model2)

## [1] 43.801

#Test whether model 2 is better than model 1
anova (model1, model2)

## Analysis of Variance Table
## 
## Model 1: pH ~ Elevation
## Model 2: pH ~ Elevation + EL2
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1     22 9.8210                                
## 2     21 6.2456  1    3.5753 12.021 0.002304 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here we can see that model 2 is better than the model 1 both based on the anova test and the AIC values. Now lets plot both regressions.

#Linear 
ggplot(Surf_Env, aes(x=Elevation, y=pH, color=Cave, shape = Habitat)) + geom_smooth(method=lm, se = TRUE, fill="light gray", colour="black", aes(group = Habitat)) + theme_classic() + scale_color_manual(values=c("#990000","#ed5e21", "#edb64e", "#a4eda1", "#f041c4", "#87dee8", "#1f60ed", "#042b80" ))+ geom_jitter(aes(colour = Cave), size = 3.5, width = 0.5) + scale_shape_manual(values=c(17,15,1))

## `geom_smooth()` using formula 'y ~ x'

#Polynomial
ggplot(Surf_Env, aes(x=Elevation, y=pH, color=Cave, shape = Habitat)) + geom_smooth(method=lm,formula =y~poly(x,2), se = TRUE, fill="light gray", color = "black", aes(group = Habitat)) + theme_classic() + scale_color_manual(values=c("#990000","#ed5e21", "#edb64e", "#a4eda1", "#f041c4", "#87dee8", "#1f60ed", "#042b80" ))+ geom_jitter(aes(colour = Cave), size = 3.5, width = 0.5) + scale_shape_manual(values=c(17,15,1))

We repeated the same analyses to investigate associations between environmental variables and decomposition parameters (k and S).

Association microbial alpha diversities with enviroenemntal variables and decomposition rates.

Associations between microbial alpha diversities and abiotic and decomposition parameters were conducted using the same methods as above.