Chapter 10 Does the antibiotic administration change the microbial community?

10.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point %in% c("Acclimation", "Antibiotics") )

label_map <- c(
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity")

alpha_div%>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic")))%>%
  #time_point = factor(time_point, levels = c("Wild", "Acclimation"))) %>%
  filter(time_point %in% c("Acclimation", "Antibiotics")) %>% 
  ggplot(aes(y = value, x = time_point, color=species, fill=species)) +
  geom_boxplot(width = 0.5,alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  geom_jitter(width = 0.2,show.legend = FALSE) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#008080', "#d57d2c")) +
  facet_grid(metric ~ species, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.text = element_text(size = 14),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 18, lineheight = 0.6),
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ time_point*species+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(5.1645)  ( log )
Formula: richness ~ time_point * species + (1 | individual)
   Data: alpha_div_meta

      AIC       BIC    logLik -2*log(L)  df.resid 
    434.8     446.1    -211.4     422.8        43 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.61978 -0.56608  0.06074  0.51883  2.39922 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.2193   0.4683  
Number of obs: 49, groups:  individual, 27

Fixed effects:
                                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                                     4.4514     0.2179  20.431  < 2e-16 ***
time_pointAntibiotics                          -1.5109     0.2404  -6.286 3.25e-10 ***
speciesPodarcis_muralis                        -0.7055     0.2719  -2.595  0.00946 ** 
time_pointAntibiotics:speciesPodarcis_muralis   0.3586     0.3037   1.181  0.23768    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) tm_pnA spcsP_
tm_pntAntbt -0.452              
spcsPdrcs_m -0.802  0.367       
tm_pntAn:P_  0.358 -0.790 -0.461

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ time_point*species, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  342.5018 353.3418 -165.2509

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.282619 7.923552

Fixed effects:  neutral ~ time_point * species 
                                                  Value Std.Error DF   t-value p-value
(Intercept)                                    44.09058  2.858870 25 15.422378   0e+00
time_pointAntibiotics                         -32.38019  3.866707 20 -8.374099   0e+00
speciesPodarcis_muralis                       -24.97345  3.534826 25 -7.064972   0e+00
time_pointAntibiotics:speciesPodarcis_muralis  20.91701  4.793363 20  4.363745   3e-04
 Correlation: 
                                              (Intr) tm_pnA spcsP_
time_pointAntibiotics                         -0.631              
speciesPodarcis_muralis                       -0.809  0.510       
time_pointAntibiotics:speciesPodarcis_muralis  0.509 -0.807 -0.632

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8258256 -0.5685215 -0.1867832  0.5823711  2.0651455 

Number of Observations: 49
Number of Groups: 27

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*species, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  203.9883 214.8282 -95.99413

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.7199475 1.694643

Fixed effects:  phylogenetic ~ time_point * species 
                                                  Value Std.Error DF   t-value p-value
(Intercept)                                    6.481542 0.6137444 25 10.560653  0.0000
time_pointAntibiotics                         -1.469110 0.8271432 20 -1.776126  0.0909
speciesPodarcis_muralis                       -1.101403 0.7588453 25 -1.451419  0.1591
time_pointAntibiotics:speciesPodarcis_muralis  0.704699 1.0254458 20  0.687212  0.4998
 Correlation: 
                                              (Intr) tm_pnA spcsP_
time_pointAntibiotics                         -0.629              
speciesPodarcis_muralis                       -0.809  0.508       
time_pointAntibiotics:speciesPodarcis_muralis  0.507 -0.807 -0.629

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8593632 -0.5146695 -0.0514960  0.4893122  1.8026133 

Number of Observations: 49
Number of Groups: 27

10.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.)


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     1 0.03530 0.035300 10.073    999  0.002 **
Residuals 47 0.16471 0.003505                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ time_point*species,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	1.9348458	0.1002347	6.065418	0.001
species	1	2.1356198	0.1106358	6.694811	0.001
time_point:species	1	0.8778553	0.0454773	2.751930	0.001
Residual	45	14.3548318	0.7436522	NA	NA
Total	48	19.3031527	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.)


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     1 0.051657 0.051657 9.9224    999  0.003 **
Residuals 47 0.244689 0.005206                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ time_point*species,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	2.0429592	0.11061666	7.254623	0.001
species	1	2.8764490	0.15574623	10.214376	0.001
time_point:species	1	0.8770558	0.04748846	3.114457	0.005
Residual	45	12.6723552	0.68614864	NA	NA
Total	48	18.4688192	1.00000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)


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     1 0.67439 0.67439 55.932    999  0.001 ***
Residuals 47 0.56670 0.01206                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ time_point*species,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	1.8783556	0.23032747	16.358708	0.001
species	1	0.7551942	0.09260332	6.577030	0.001
time_point:species	1	0.3545686	0.04347786	3.087959	0.024
Residual	45	5.1670341	0.63359135	NA	NA
Total	48	8.1551525	1.00000000	NA	NA

dbRDA

#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
    scale_fill_manual(values=c('#00808050', "#d57d2c50"))+
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()