Core microbiota of Zamia species
Zamia Coralloid root Core Cyanobacteria
Contributors
Juan Carlos Villarreal A.
Adriel M. Sierra
We use the two dataset generated with two amplicon markers to identify the core microbiota.
First we use a subset dataset with only ASVs whose phylum corresponds to Cyanobacteria in the 16S dataset.
knitr::opts_knit$set(root.dir = "/Users/adrielsierra/Documents/GitHub/cycads_rbcLX_roots")
root16ps = readRDS("./Data/Zamia_root_bacteria16s_ps_mai23.rds")
ps_cyano = subset_taxa(root16ps, Phylum=="Cyanobacteria")Second, we use the cyanobacterial community identified with the amplicon marker rbcL-X.
knitr::opts_knit$set(root.dir = "/Users/adrielsierra/Documents/GitHub/cycads_rbcLX_roots")
rbclx.ps = readRDS("./Data/RBCLX_ps_rootZamia.rds")
rbclx.ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 175 taxa and 25 samples ]
## sample_data() Sample Data: [ 25 samples by 3 sample variables ]
## tax_table() Taxonomy Table: [ 175 taxa by 6 taxonomic ranks ]16S
Zamia pseudoparasitica
cyano.zpseudo = subset_samples(ps_cyano, Species %in% c("Zamia pseudoparasitica"))
rank.pseudo16 = rankabundance(otu_table(cyano.zpseudo))
rankabunplot(rank.pseudo16, labels=" ",scale="logabun", lty=2, pch= 20)## Warning in xy.coords(x, y, xlabel, ylabel, log): 20 y values <= 0 omitted from
## logarithmic plot
Zamia nana
cyano.znana = subset_samples(ps_cyano, Species %in% c("Zamia nana"))
rank.nana16 = rankabundance(otu_table(cyano.znana))
rankabunplot(rank.nana16, labels=" ",scale="logabun", lty=2, pch= 20)## Warning in xy.coords(x, y, xlabel, ylabel, log): 76 y values <= 0 omitted from
## logarithmic plot
rbcL-X
Zamia pseudoparasitica
ps_pseudo = subset_samples(rbclx.ps, Species_name %in% c("Zamia pseudoparasitica"))
rank.pseudo = rankabundance(t(otu_table(ps_pseudo)))
rankabunplot(rank.pseudo, labels=" ",scale="logabun", lty=2, pch= 20)## Warning in xy.coords(x, y, xlabel, ylabel, log): 34 y values <= 0 omitted from
## logarithmic plot
Zamia nana
ps_nana = subset_samples(rbclx.ps, Species_name %in% c("Zamia nana"))
rank.nana = rankabundance(t(otu_table(ps_nana)))
rankabunplot(rank.nana, labels=" ",scale="logabun", lty=2, pch= 20)## Warning in xy.coords(x, y, xlabel, ylabel, log): 104 y values <= 0 omitted from
## logarithmic plot
Core microbiome.
Indicator species analysis, the labdsv package (version 2.0-1) was used to identify which ASVs were relatively more abundant and predominantly found to one of the specific Zamia species, in contrast to the other host species.
The script is adapted from https://bocasbiome.github.io/wf1.html
The Indicator species analysis was done for both marker. A core microbiota was only observed when we analyzed the dataset generated from the rbcL-X marker. Therefore, here is only presented for the rbcL-X.
zamia.host = data.frame(t(otu_table(rbclx.ps)))
cyano.ASV = tibble::remove_rownames(zamia.host)
host.group = data.frame(sample_data(rbclx.ps)) %>%
select(Species_name)
host.group$Status = host.group$Species
host.group = host.group[, c(2,1)]
host.group$Status = str_replace(host.group$Status, "Zamia nana", "1")
host.group$Status = str_replace(host.group$Status, "Zamia pseudoparasitica", "2")
host.group = tibble::rownames_to_column(host.group, "Label")
host.group$Status = as.integer(host.group$Status)
host.group$Species = as.character(host.group$Species)Here we calculate the indicator values.
set.seed(1280)
iva = indval(cyano.ASV, host.group$Status)
gr = iva$maxcls[iva$pval <= 0.01]
iv = iva$indcls[iva$pval <= 0.01]
pv = iva$pval[iva$pval <= 0.01]
fr = apply(cyano.ASV > 0, 2, sum)[iva$pval <= 0.01]
indval.out = data.frame(group = gr, indval = iv, pval = pv, freq = fr)
indval.out = indval.out[order(indval.out$group, -indval.out$indval),]
indval.out## group indval pval freq
## ASV_2 2 0.74343 0.004 21indval.out$prob.corrected = p.adjust(indval.out$pval, "bonferroni")
indval.out## group indval pval freq prob.corrected
## ASV_2 2 0.74343 0.004 21 0.004core microbiome using threshold frequency
For the threshold frequency base method, we used the microbiome R package (Lahti and Shetty 2012). Due to the lack of consensus about a fixed threshold (Risely 2020). We tested for the presence of a ‘core’ microbiota as defined by any taxon with a prevalence higher than 50%, 75% or 90% (Jorge et al. 2020; Neu et al. 2021). We pruned out the phyloseq object to retain those samples with more than 50 reads. We detected the core microbiota with the function (core_members) and estimated the total abundance of specific members of the core microbiota in each sample (sample_sums).
To allow for comparisons of the core microbiomes, we considered samples that we have amplicon sequencing data for the two markers and from the same locality.
ps_cyano_sub = subset_samples(ps_cyano, Locality !="El Cope")
ps_cyano_prun = prune_taxa(taxa_sums(ps_cyano_sub)!=0, ps_cyano_sub)First, we pruned samples with less than 50 reads, and transform incidence of ASVs in the samples to compositional relative abundance.
which(!rowSums(otu_table(rbclx.ps)) > 50)## ASV_145 ASV_146 ASV_147 ASV_148 ASV_149 ASV_150 ASV_151 ASV_152 ASV_153 ASV_154
## 145 146 147 148 149 150 151 152 153 154
## ASV_155 ASV_156 ASV_157 ASV_158 ASV_159 ASV_160 ASV_161 ASV_162 ASV_163 ASV_164
## 155 156 157 158 159 160 161 162 163 164
## ASV_166 ASV_167 ASV_168 ASV_170 ASV_171 ASV_172 ASV_173 ASV_176 ASV_177 ASV_180
## 165 166 167 168 169 170 171 172 173 174
## ASV_183
## 175total = median(sample_sums(rbclx.ps))
standf = function(x, t=total) round(t * (x / sum(x)))
rbclx.ps = transform_sample_counts(rbclx.ps, standf)
ps2prun.ra = microbiome::transform(rbclx.ps, "compositional")
ps2prun.ra2 = prune_taxa(taxa_sums(ps2prun.ra) > 0, ps2prun.ra)
taxonomy = as.data.frame(tax_table(ps2prun.ra2))We first extract the taxonomy table. Subset this taxonomy table to include only core ASVs with the different threshold. Total core abundance in each sample (sum of abundances of the core members).
core50_name = core_members(ps2prun.ra2, detection = 0.0001, prevalence = 50/100, include.lowest = FALSE)
core50_seq = core_members(ps2prun.ra2, detection = 0.0001, prevalence = 50/100, include.lowest = FALSE)
core.abundance.50 = sample_sums(core(ps2prun.ra2, detection = 0.0001, prevalence = 50/100))
DT::datatable(as.data.frame(core.abundance.50))core_taxa_id50 = subset(taxonomy, rownames(taxonomy) %in% core50_name)
DT::datatable(core_taxa_id50)core75_name = core_members(ps2prun.ra2, detection = 0.0001, prevalence = 75/100, include.lowest = FALSE)
core75_seq = core_members(ps2prun.ra2, detection = 0.0001, prevalence = 75/100, include.lowest = FALSE)
core.abundance.75 = sample_sums(core(ps2prun.ra2, detection = 0.0001, prevalence = 75/100))
DT::datatable(as.data.frame(core.abundance.75))core_taxa_id75 = subset(taxonomy, rownames(taxonomy) %in% core75_name)
DT::datatable(core_taxa_id75)core90_name = core_members(ps2prun.ra2, detection = 0.0001, prevalence = 90/100, include.lowest = FALSE)
core90_seq = core_members(ps2prun.ra2, detection = 0.0001, prevalence = 90/100, include.lowest = FALSE)
core.abundance.90 = sample_sums(core(ps2prun.ra2, detection = 0.0001, prevalence = 90/100))
DT::datatable(as.data.frame(core.abundance.90))core_taxa_id90 = subset(taxonomy, rownames(taxonomy) %in% core90_name)
DT::datatable(core_taxa_id90)Visualitation of the core microbiome with Core line plots considering the compositional (relative) abundances.
det = c(0, 0.1, 0.5, 2, 5, 20)/100
prevalences = seq(.05, 1, .05)
plot_core(ps2prun.ra2, prevalences = prevalences,
detections = det, plot.type = "lineplot") +
xlab("Relative Abundance (%)") +
theme_bw()## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
prevalences = seq(.05, 1, .05)
detections = round(10^seq(log10(1e-3), log10(.2), length = 10), digits = 4)
p.core = plot_core(ps2prun.ra2,
plot.type = "heatmap",
colours = rev(brewer.pal(5, "Spectral")),
prevalences = prevalences,
detections = detections,
min.prevalence = .5) + theme_classic() +
xlab("Detection Threshold (Relative Abundance (%))")
p.core 