Zamia Coralloid root bacterial community

Contributors

Juan Carlos Villarreal A.

Adriel M. Sierra

To study the overall coralloid root bacterial community of Zamia species, we used the processed phyloseq object with the whole bacterial community.

The dataset include a total of fourty four coralloid roots of Z. pseudoparasitica from two localities, and ten of Zamia nana from one locality

knitr::opts_knit$set(root.dir = "/Users/adrielsierra/Documents/GitHub/cycads_rbcLX_roots")
root16ps = readRDS("./Data/Zamia_root_bacteria16s_ps_mai23.rds")
root16ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1561 taxa and 53 samples ]
## sample_data() Sample Data:       [ 53 samples by 5 sample variables ]
## tax_table()   Taxonomy Table:    [ 1561 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1561 tips and 1559 internal nodes ]

Compute and visualize alpha diversity

Alpha diversity was assessed using four indexes: i) the Chao index which shows the observed richness, ii) the exponential of Shannon which weights the ASVs by their frequency, iii) the Simpson’s diversity index which considers the presence and relative abundances of ASVs present in a sample, and iv) the Pielou index which measures diversity along with species richness.

richness.df = estimate_richness(root16ps, split = TRUE, measures = c("Observed","Shannon", "Simpson"))

richness.df$Pielou = (richness.df$Shannon)/log(specnumber(otu_table(root16ps)))

Now, we add extra column with read count by sample, the Zamia host species and locality respective to each sample

richness.df$readcount = sample_sums(root16ps)

richness.df$Host = as.factor(root16ps@sam_data$Species)

richness.df$Locality = as.factor(root16ps@sam_data$Locality)

Alpha diversity table : Alpha diversity (richness indexes) for the complete bacterial community in coralloid roots of the two species Z. nana and Z. pseudoparasitica.

 DT::datatable(richness.df) %>%
formatRound(columns = c("Shannon", "Simpson", "Pielou"), digits=4)

The assumptions of normality was tested with the Shapiro-Wilk Test

shapiro.test(richness.df$Observed)
## 
##  Shapiro-Wilk normality test
## 
## data:  richness.df$Observed
## W = 0.94207, p-value = 0.01244
shapiro.test(richness.df$Shannon)
## 
##  Shapiro-Wilk normality test
## 
## data:  richness.df$Shannon
## W = 0.93602, p-value = 0.007031
shapiro.test(richness.df$Simpson)
## 
##  Shapiro-Wilk normality test
## 
## data:  richness.df$Simpson
## W = 0.74905, p-value = 3.715e-08
shapiro.test(richness.df$Pielou)
## 
##  Shapiro-Wilk normality test
## 
## data:  richness.df$Pielou
## W = 0.81468, p-value = 1.338e-06

The homogeneity of variance were tested with the Levene test

test_names <- c( "Observed", "Shannon", "Simpson", "Pielou")
for (test in test_names) {
  test_formula <- paste(test,"~","Host")
  print(paste("Testing homogeneity of variance:", test))
  print(leveneTest(as.formula(test_formula), data=richness.df))
}
## [1] "Testing homogeneity of variance: Observed"
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.2698 0.2651
##       51               
## [1] "Testing homogeneity of variance: Shannon"
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  5.0997 0.02824 *
##       51                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Testing homogeneity of variance: Simpson"
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.4914 0.06744 .
##       51                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Testing homogeneity of variance: Pielou"
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.6555 0.06162 .
##       50                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Statistical test for significance between host species and study localities

Here, the significant differences between species and markers were assessed using the Kruskal-Wallis nonparametric test. Post hoc pairwise comparisons were performed with the exact sum of Wilcoxon ranks.

for (test in test_names) {
  test_formula <- paste(test,"~","Host")
  print(paste("Alpha diversity:", test))
  print(kruskal.test(as.formula(test_formula), data=richness.df))
}
## [1] "Alpha diversity: Observed"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Observed by Host
## Kruskal-Wallis chi-squared = 0.81082, df = 1, p-value = 0.3679
## 
## [1] "Alpha diversity: Shannon"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by Host
## Kruskal-Wallis chi-squared = 2.1572, df = 1, p-value = 0.1419
## 
## [1] "Alpha diversity: Simpson"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Simpson by Host
## Kruskal-Wallis chi-squared = 2.0881, df = 1, p-value = 0.1485
## 
## [1] "Alpha diversity: Pielou"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Pielou by Host
## Kruskal-Wallis chi-squared = 2.434, df = 1, p-value = 0.1187
for (test in test_names) {
  test_formula <- paste(test,"~","Locality")
  print(paste("Alpha diversity:", test))
  print(kruskal.test(as.formula(test_formula), data=richness.df))
}
## [1] "Alpha diversity: Observed"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Observed by Locality
## Kruskal-Wallis chi-squared = 0.91613, df = 2, p-value = 0.6325
## 
## [1] "Alpha diversity: Shannon"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Shannon by Locality
## Kruskal-Wallis chi-squared = 4.3259, df = 2, p-value = 0.115
## 
## [1] "Alpha diversity: Simpson"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Simpson by Locality
## Kruskal-Wallis chi-squared = 7.5496, df = 2, p-value = 0.02294
## 
## [1] "Alpha diversity: Pielou"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Pielou by Locality
## Kruskal-Wallis chi-squared = 4.1152, df = 2, p-value = 0.1278

Alpha diversity Figures : Alpha diversity (richness indexes) for the complete bacterial community in coralloid roots of the two species Z. nana and Z. pseudoparasitica.

Figsupp_A = 
  ggplot(richness.df, aes(x=Host, y=Observed, fill=Locality)) + 
  ggtitle("Alpha Diversity", subtitle="Complete bacterial community")+
  scale_fill_brewer(palette="Pastel2")+
  geom_boxplot(width = 0.3)+
  #scale_fill_grey()+
  theme_classic()+
  theme(text = element_text(size=12))+
  labs(x=" ", y = "Observed species")

Figsupp_A

Figsupp_B = 
  ggplot(richness.df, aes(x=Host, y=Shannon, fill=Locality)) + 
  #ggtitle("Alpha Diversity", subtitle="Complete bacterial community")+
  scale_fill_brewer(palette="Pastel2")+
  geom_boxplot(width = 0.3)+
  #scale_fill_grey()+
  theme_classic()+ theme(text = element_text(size=12))+
  labs(x=" ", y = "Shannon")

Figsupp_B

Figsupp_C = 
  ggplot(richness.df, aes(x=Host, y=Simpson, fill=Locality)) + 
  #ggtitle("Alpha Diversity", subtitle="Complete bacterial community")+
  scale_fill_brewer(palette="Pastel2")+
  geom_boxplot(width = 0.3)+
  #scale_fill_grey()+
  theme_classic()+ theme(text = element_text(size=12))+
  labs(x=" ", y = "Simpson")
  
Figsupp_C

Figsupp_D = 
  ggplot(richness.df, aes(x=Host, y=Pielou, fill=Locality)) + 
  #ggtitle("Alpha Diversity", subtitle="Complete bacterial community")+
  scale_fill_brewer(palette="Pastel2")+
  geom_boxplot(width = 0.3)+
  #scale_fill_grey()+
  theme_classic()+ theme(text = element_text(size=12))+
  labs(x=" ", y = "Pielou")

Figsupp_D
## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).