EEID Adult Newts Temperature Experiments: Microbiomes (16S & ITS)

Packages

library(tidyverse)
library(ggplot2)
library(ggthemes)
library(RColorBrewer)
library(lme4)
library(qiime2R)
library(ecodist)
library(vegan)

16S Bacterial Communities

Data Compiling

Data Compiling

Metadata merging of antifungal function data

Metadata_AntiFungal_PCO_Alpha = read_tsv("16S_Metadata_AlphaDiv_AntiFungalPred_TimePre.txt")

Alpha diversity

Alpha diversity metrics were calculated in QIIME2 with the diversity core-metrics-phylogenetic command with a rarefaction depth of 5000 reads

Alpha Models

sOTU Richness models

# dist check
hist(Metadata_AntiFungal_PCO_Alpha$observed_otus)

# main model
mod_sOTU_1 <- glmer(observed_otus ~ Temperature + (1|CollectionLocation), family = poisson(link = "log"),data = Metadata_AntiFungal_PCO_Alpha)

car::Anova(mod_sOTU_1, type=c("III"), test.statistic=c("Chisq"))

# model fit
qqnorm(resid(mod_sOTU_1))
qqline(resid(mod_sOTU_1))

# post hoc test
emmeans::emmeans(mod_sOTU_1, list(pairwise ~ Temperature), adjust = "tukey")

## $`emmeans of Temperature`
##  Temperature emmean     SE  df asymp.LCL asymp.UCL
##         13.8   4.56 0.0439 Inf      4.47      4.64
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $` of Temperature`
##  contrast  estimate SE df z.ratio p.value
##  (nothing)   nonEst NA NA NA      NA     
## 
## Note: contrasts are still on the log scale

Phylogenetic diversity models

# dist check
hist(Metadata_AntiFungal_PCO_Alpha$faith_pd)

# main model
mod_PD_1 <-glmer(faith_pd ~ Temperature + (1|CollectionLocation), family = Gamma, data = Metadata_AntiFungal_PCO_Alpha)

car::Anova(mod_PD_1, type=c("III"), test.statistic=c("Chisq"))

qqnorm(resid(mod_PD_1))
qqline(resid(mod_PD_1))

# Post hoc test
emmeans::emmeans(mod_PD_1, list(pairwise ~ Temperature), adjust = "tukey")

## $`emmeans of Temperature`
##  Temperature emmean      SE  df asymp.LCL asymp.UCL
##         13.8  0.123 0.00771 Inf     0.108     0.138
## 
## Results are given on the inverse (not the response) scale. 
## Confidence level used: 0.95 
## 
## $` of Temperature`
##  contrast  estimate SE df z.ratio p.value
##  (nothing)   nonEst NA NA NA      NA     
## 
## Note: contrasts are still on the inverse scale

Alpha diversity across temperatures plot

trt_order <- c("6", "14", "22")

Metadata_AntiFungal_PCO_Alpha %>% 
  ggplot(aes(x = factor(Temperature, level = trt_order), y = observed_otus, color = factor(Temperature, level = trt_order))) +
  geom_violin()+
  geom_point(size = 3, alpha = 0.6) +
  scale_color_manual(values = c("#087CAF", "#E7B800", "#DB6E0B"),name = "Temperature (C)") +
  theme_classic() +
  xlab("Temperature (C)") +
  ylab("sOTU Richness") +
  ggtitle("Richness of adult newt bacterial communities across temperatures")

Beta diversiy

Beta diversity models

Beta diversity metrics were calculated in QIIME2 with the diversity core-metrics-phylogenetic command with a rarefaction depth of 5000 reads

Weighted Unifrac models

#read in dist matrix
AdNewt_WUF<-read_qza("16S_weighted_unifrac_distance_matrix.qza")

AdNewt_WUF_mx <- as.matrix(AdNewt_WUF$data)

# filtering dist matrix and metadata to same samples
AdNewtT0_distList2<-colnames(AdNewt_WUF_mx)

Metadata_AntiFungal_PCO_Alpha_T0_BD2<- Metadata_AntiFungal_PCO_Alpha %>%
  filter(SampleID %in% AdNewtT0_distList2)

metaList<-as.vector(Metadata_AntiFungal_PCO_Alpha_T0_BD2$SampleID)
AdNewt_WUF_mx2 <- usedist::dist_subset(AdNewt_WUF_mx,metaList)
AdNewt_WUF_mx3 <-as.matrix(AdNewt_WUF_mx2)

# main model
adonis(AdNewt_WUF_mx3 ~ Temperature, data = Metadata_AntiFungal_PCO_Alpha_T0_BD2)

## 
## Call:
## adonis(formula = AdNewt_WUF_mx3 ~ Temperature, data = Metadata_AntiFungal_PCO_Alpha_T0_BD2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##              Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Temperature   1    0.6124 0.61238   13.85 0.09975  0.001 ***
## Residuals   125    5.5269 0.04422         0.90025           
## Total       126    6.1393                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise models

# install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")

# pairwise post hoc
WUFmodel <- pairwiseAdonis::pairwise.adonis(AdNewt_WUF_mx3, Metadata_AntiFungal_PCO_Alpha_T0_BD2$Temperature, p.adjust.m = "bonferroni")

Unweighted Unifrac models

# read in dist matrix

AdNewt_unWUF<-read_qza("16S_unweighted_unifrac_distance_matrix.qza")
AdNewt_unWUF_mx <- as.matrix(AdNewt_unWUF$data)

# filtering dist matrix and metadata to same samples

AdNewtT0_distList3<-colnames(AdNewt_unWUF_mx)
Metadata_AntiFungal_PCO_Alpha_T0_BD3<- Metadata_AntiFungal_PCO_Alpha %>% 
 filter(SampleID %in% AdNewtT0_distList3)

metaList3<-as.vector(Metadata_AntiFungal_PCO_Alpha_T0_BD3$SampleID)
AdNewt_unWUF_mx2 <- usedist::dist_subset(AdNewt_unWUF_mx,metaList3)
AdNewt_unWUF_mx3 <-as.matrix(AdNewt_unWUF_mx2)

# main model
adonis(AdNewt_unWUF_mx3 ~ Temperature, data = Metadata_AntiFungal_PCO_Alpha_T0_BD2)

## 
## Call:
## adonis(formula = AdNewt_unWUF_mx3 ~ Temperature, data = Metadata_AntiFungal_PCO_Alpha_T0_BD2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##              Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Temperature   1    4.8776  4.8776  32.301 0.20534  0.001 ***
## Residuals   125   18.8757  0.1510         0.79466           
## Total       126   23.7533                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise models

# post hoc
pairwiseAdonis::pairwise.adonis(as.dist(AdNewt_unWUF_mx3), Metadata_AntiFungal_PCO_Alpha_T0_BD2$Temperature, p.adjust.m = "bonferroni")

PCoA Plots

trt_order <- c("6", "14", "22")

uWUF_16S_PCoA <- ape::pcoa(AdNewt_unWUF_mx3)

PCoA_Axes16S = as.data.frame(uWUF_16S_PCoA$vectors)
PCoA_Axes16S =tibble::rownames_to_column(PCoA_Axes16S, "SampleID")

PCoA_Axes16S %>%
  dplyr::select(SampleID, Axis.1, Axis.2) %>%
  left_join(Metadata_AntiFungal_PCO_Alpha) %>%
  filter(TimeWeek=="0")%>%
  ggplot(aes(x=Axis.1, y=Axis.2, color = as.factor(Temperature),size = Propor_TotalAntiFungal_BSALonly)) +
  geom_point(aes(color = as.factor(Temperature), alpha=0.5))+
   scale_color_manual(values = c("#087CAF", "#E7B800", "#DB6E0B"),name = "Temperature (C)") +
  stat_ellipse(type = "t",linetype = 2)+
   theme_classic() +
  xlab("PCo1 - 4.5%") +
  ylab("PCo2 - 2.8%")

Taxa Plots

Taxonomic Composition

Grouped sOTU tables with one column per temperature were calculated in QIIME2 with the feature-table group command using a individual sample OTU table rarified to 5000 reads and converted to a tab-delimited file using biom convert

TempTable2 = read_tsv("16S_TempGrouped_rarefied_table_r5000.txt")

Taxa = read_tsv("16S_Taxonomy.txt") %>% parse_taxonomy() %>%
  rownames_to_column("OTUID")

TempTable2Taxa = TempTable2 %>%
  left_join(.,Taxa, by = "OTUID")

TempTable2Taxa_Order = TempTable2Taxa %>%
  group_by(Order) %>%
  summarize(AvgAbund6 = sum(t6),AvgAbund14 = sum(t14),AvgAbund22 = sum(t22), totalAbund = sum(t14+t6+t22))

## select rows that are greater than 300 in total abundance --> move to new table
TempTable2Taxa_Order_Top = TempTable2Taxa_Order %>%
  filter(totalAbund>=300) %>%
  select( -totalAbund)

TempTable2Taxa_Order_Top <- as.data.frame(TempTable2Taxa_Order_Top)

## select rows that are less than 300 and move to new table --> sum all to new table
TempTable2Taxa_Order_Low = TempTable2Taxa_Order %>%
  filter(totalAbund<300) %>%
  summarize(AvgAbund6 = sum(AvgAbund6),AvgAbund14 = sum(AvgAbund14),AvgAbund22 = sum(AvgAbund22)) %>%
  mutate(Order = "Other")
TempTable2Taxa_Order_Low <- as.data.frame(TempTable2Taxa_Order_Low)

## merge to tables and make it long format
TempTable2Taxa_Order_Summarize <- full_join(TempTable2Taxa_Order_Top,TempTable2Taxa_Order_Low)


TempTable2Taxa_Order_Summarize_long = TempTable2Taxa_Order_Summarize %>%
  gather(key = Temperature, value = Abundance, -Order)

TempTable2Taxa_Order_Summarize_long_prop = TempTable2Taxa_Order_Summarize_long %>% 
  group_by(Temperature) %>% 
   mutate(total = sum(Abundance), rel.freq = Abundance / total)

trt_order2 <- c("AvgAbund6", "AvgAbund14", "AvgAbund22")

TempTable2Taxa_Order_Summarize_long_prop %>%
  ggplot(aes(x = factor(Temperature, level = trt_order2), y = rel.freq, fill = Order)) +
  geom_bar(stat = "identity")+
  scale_fill_stata() +
  ylab("Relative Abundance")+
  xlab("Temperature") +
  scale_x_discrete(labels=c("AvgAbund6" = "6", "AvgAbund22" = "22", "AvgAbund14" = "14"))

ITS Fungal Communities

Alpha Diversity

Alpha diversity metrics were calculated in QIIME2 with the diversity core-metrics-phylogenetic command with a rarefaction depth of 5000 reads

ITSMeta_Alpha = read_tsv("ITS_Metadata_AlphaDiv_TimePre.txt")

Alpha Diversity models

Note: included Collection location as a random effect in alpha diversity models.

#dist check
hist(ITSMeta_Alpha$observed_otus)

#main model
mod_sOTU_ITS1 <- glmer(observed_otus ~ Temperature + (1|CollectionLocation), family = Gamma,data = ITSMeta_Alpha)

car::Anova(mod_sOTU_ITS1, type=c("III"), test.statistic=c("Chisq"))

qqnorm(resid(mod_sOTU_ITS1))
qqline(resid(mod_sOTU_ITS1))

# post hoc test
emmeans::emmeans(mod_sOTU_ITS1, list(pairwise ~ Temperature), adjust = "tukey")

## $`emmeans of Temperature`
##  Temperature emmean      SE  df asymp.LCL asymp.UCL
##         14.7 0.0347 0.00554 Inf    0.0239    0.0456
## 
## Results are given on the inverse (not the response) scale. 
## Confidence level used: 0.95 
## 
## $` of Temperature`
##  contrast  estimate SE df z.ratio p.value
##  (nothing)   nonEst NA NA NA      NA     
## 
## Note: contrasts are still on the inverse scale

Alpha Diversity plot

Temptrt_order <- c("6", "14", "22")

ITSMeta_Alpha %>% 
  filter(observed_otus<200) %>%
  ggplot(aes(x = factor(Temperature, level = Temptrt_order), y = observed_otus, color = factor(Temperature, level = Temptrt_order))) +
  geom_violin()+
  geom_point(size = 3 , alpha = 0.6) +
  scale_color_manual(values = c("#087CAF", "#E7B800", "#DB6E0B"),name = "Temperature (C)") +
  theme_classic() +
  xlab("Temperature (C)") +
  ylab("ITS sOTU Richness") +
  ggtitle("Richness of adult newt fungal communities across temperatures")

Beta diversity

Beta diversity

Beta diversity metrics were calculated in QIIME2 with the diversity core-metrics-phylogenetic command with a rarefaction depth of 5000 reads

#read in dist matrix
AdNewt_JacITSmx<-read_qza("ITS_jaccard_distance_matrix.qza")

AdNewt_JacITSmx2 <- as.matrix(AdNewt_JacITSmx$data)

#filtering distance matrix and metadata

AdNewtT0_distList2<-colnames(AdNewt_JacITSmx2)
ITSMeta_ObsOTUS_AntiF_0b<- ITSMeta_Alpha %>%
  filter(SampleID %in% AdNewtT0_distList2)

metaList<-as.vector(ITSMeta_ObsOTUS_AntiF_0b$SampleID)

AdNewt_Jac_mx2c <- usedist::dist_subset(AdNewt_JacITSmx2,metaList)
AdNewt_JacITSmx3 <-as.matrix(AdNewt_JacITSmx2)

# main model
adonis(AdNewt_JacITSmx3 ~ Temperature, data = ITSMeta_Alpha)

## 
## Call:
## adonis(formula = AdNewt_JacITSmx3 ~ Temperature, data = ITSMeta_Alpha) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Temperature  1     2.125 2.12529  4.7978 0.05226  0.001 ***
## Residuals   87    38.539 0.44297         0.94774           
## Total       88    40.664                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# dispersion model
AdNewt_JacITSmx3_disp<-betadisper(as.dist(AdNewt_JacITSmx3), ITSMeta_Alpha$Temperature, type = c("median"), bias.adjust = FALSE, sqrt.dist = FALSE, add = FALSE)

anova(AdNewt_JacITSmx3_disp)

Pairwise models

#devtools::install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")

model <- pairwiseAdonis::pairwise.adonis(AdNewt_JacITSmx3, ITSMeta_Alpha$Temperature, p.adjust.m = "bonferroni")

Beta diversity plot

#read in dist matrix

PCoA_ITS <- ape::pcoa(AdNewt_JacITSmx3)

PCoA_AxesITS = as.data.frame(PCoA_ITS$vectors)
PCoA_AxesITS =tibble::rownames_to_column(PCoA_AxesITS, "SampleID")

PCoA_AxesITS %>%
  dplyr::select(SampleID, Axis.1, Axis.2) %>%
  left_join(ITSMeta_Alpha) %>%
  filter(TimeWeek=="0")%>%
  ggplot(aes(x=Axis.1, y=Axis.2, color = as.factor(Temperature))) +
  geom_point(aes(color = as.factor(Temperature), alpha=0.5))+
   scale_color_manual(values = c("#087CAF", "#E7B800", "#DB6E0B"),name = "Temperature (C)") +
  stat_ellipse(type = "t",linetype = 2)+
   theme_classic() +
  xlab("PCo1 - 4.5%") +
  ylab("PCo2 - 2.8%")

Taxa Plots

ITS Taxonomic Composition

Grouped sOTU tables with one column per temperature were calculated in QIIME2 with the feature-table group command using a individual sample OTU table rarified to 1000 reads and converted to a tab-delimited file using biom convert

TempTable2ITS = read_tsv("ITS_TempGrouped_rarefied_table_r1000.txt")

TaxaITS = read_tsv("ITS_Taxonomy.txt") %>% parse_taxonomy() %>%
  rownames_to_column("OTUID")

TempTable2ITSTaxa = TempTable2ITS %>%
  left_join(.,TaxaITS, by = "OTUID")

TempTable2ITSTaxa_Order = TempTable2ITSTaxa %>%
  group_by(Order) %>%
  summarize(AvgAbund6 = sum(t6),AvgAbund14 = sum(t14),AvgAbund22 = sum(t22), totalAbund = sum(t14+t6+t22))

## select rows that are greater than 100 in total abundance --> move to new table
TempTable2ITSTaxa_Order_Top = TempTable2ITSTaxa_Order %>%
  filter(totalAbund>=100) %>%
  select( -totalAbund)

TempTable2ITSTaxa_Order_Top <- as.data.frame(TempTable2ITSTaxa_Order_Top)

## select rows that are less than 100 and move to new table --> sum all to new table
TempTable2ITSTaxa_Order_Low = TempTable2ITSTaxa_Order %>%
  filter(totalAbund<100) %>%
  summarize(AvgAbund6 = sum(AvgAbund6),AvgAbund14 = sum(AvgAbund14),AvgAbund22 = sum(AvgAbund22)) %>%
  mutate(Order = "Other")
TempTable2ITSTaxa_Order_Low <- as.data.frame(TempTable2ITSTaxa_Order_Low)

## merge to tables and make it long format
TempTable2ITSTaxa_Order_Summarize <- full_join(TempTable2ITSTaxa_Order_Top,TempTable2ITSTaxa_Order_Low)


TempTable2ITSTaxa_Order_Summarize_long = TempTable2ITSTaxa_Order_Summarize %>%
  gather(key = Temperature, value = Abundance, -Order)

TempTable2ITSTaxa_Order_Summarize_long_prop = TempTable2ITSTaxa_Order_Summarize_long %>% 
  group_by(Temperature) %>% 
   mutate(total = sum(Abundance), rel.freq = Abundance / total) %>% 
  replace_na(list(Order = "Unknown Fungi"))

trt_order2 <- c("AvgAbund6", "AvgAbund14", "AvgAbund22")

as.data.frame(TempTable2ITSTaxa_Order_Summarize_long_prop) %>%
  ggplot(aes(x = factor(Temperature, level = trt_order2), y = rel.freq, fill = Order)) +
  geom_bar(stat = "identity")+
  scale_fill_stata() +
  ylab("Relative Abundance")+
  xlab("Temperature") +
  scale_x_discrete(labels=c("AvgAbund6" = "6", "AvgAbund22" = "22", "AvgAbund14" = "14")) +
  theme_classic()