Chapter 12 Does the inoculum sample maintain the microbial community found in acclimation?

sample_metadata <- sample_metadata %>%
  mutate(treatment=time_point) %>% 
    mutate(treatment=case_when(
    treatment %in% c("Inoculum1", "Inoculum2") ~ "Inoculum",
    TRUE ~ treatment
  ))

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Warm_control" & treatment %in% c("Acclimation","Inoculum"))

samples_to_keep <- sample_metadata %>%
  filter(type == "Warm_control" & treatment %in% c("Acclimation","Inoculum")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Warm_control" & treatment %in% c("Acclimation","Inoculum"))

12.1 Alpha diversity

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., alpha_div_meta, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
  filter(type == "Warm_control" & treatment %in% c("Acclimation","Inoculum")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  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("#d57d2c", "#d5b52c")) +
  scale_fill_manual(values=c("#d57d2c", "#d5b52c")) +
  facet_wrap( ~ metric, scales = "free")+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.text = element_text(size = 12),
    axis.title.y = element_text(size = 14),
    axis.title.x = element_blank(),
    strip.text = element_text(size = 16, lineheight = 0.6)
      ) +
  ylab("Alpha diversity")

Richness

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

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

      AIC       BIC    logLik -2*log(L)  df.resid 
    228.0     232.7    -110.0     220.0        20 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.37872 -0.53180  0.06467  0.46904  1.76837 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0        0       
Number of obs: 24, groups:  individual, 9

Fixed effects:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)         4.4569     0.0834  53.441   <2e-16 ***
treatmentInoculum   0.1855     0.1049   1.769   0.0769 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
trtmntInclm -0.795
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Neutral

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

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  185.3367 189.7009 -88.66837

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev: 0.001923638 12.18199

Fixed effects:  neutral ~ treatment 
                     Value Std.Error DF   t-value p-value
(Intercept)       44.09058  4.060663 14 10.857976  0.0000
treatmentInoculum  1.80350  5.136377 14  0.351122  0.7307
 Correlation: 
                  (Intr)
treatmentInoculum -0.791

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.66610957 -0.48314824 -0.03902931  0.62479437  2.02597977 

Number of Observations: 24
Number of Groups: 9

Phylogenetic

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

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  80.84515 85.20932 -36.42257

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:    0.722741 0.9632749

Fixed effects:  phylogenetic ~ treatment 
                      Value Std.Error DF   t-value p-value
(Intercept)        6.481542 0.4014215 14 16.146474  0.0000
treatmentInoculum -0.579687 0.4129656 14 -1.403716  0.1822
 Correlation: 
                  (Intr)
treatmentInoculum -0.622

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-2.329639154 -0.568062343 -0.001994187  0.552203432  1.529230884 

Number of Observations: 24
Number of Groups: 9

12.2 Beta diversity

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$treatment) %>% 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.00836 0.0083600 1.4409    999  0.253
Residuals 22 0.12764 0.0058019

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.2657351	0.06396199	1.503319	0.107
Residual	22	3.8888430	0.93603801	NA	NA
Total	23	4.1545781	1.00000000	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$treatment) %>% 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.000661 0.0006614 0.0736    999  0.815
Residuals 22 0.197804 0.0089911

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3164816	0.07596593	1.808646	0.079
Residual	22	3.8496178	0.92403407	NA	NA
Total	23	4.1660995	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$treatment) %>% 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.00204 0.0020405 0.2622    999  0.675
Residuals 22 0.17118 0.0077807

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.0412056	0.056244	1.31111	0.242
Residual	22	0.6914166	0.943756	NA	NA
Total	23	0.7326222	1.000000	NA	NA

NMDS

richness %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 3) +
  scale_color_manual(name="Time point",values=c("#d57d2c", "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.6, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 3) +
  scale_color_manual(name="Time point",values=c("#d57d2c", "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.6, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

phylo %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 3) +
  scale_color_manual(name="Time point",values=c("#d57d2c", "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.6, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )