Chapter 6 Bacterial community composition

6.1 Load MG data

load("data/data_mg.Rdata")

6.2 Hill numbers - beta diversity

# Neutral
beta_q1n <- 
  mag_weighted %>% 
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, metric = "S", out = "pair")

# Phylogenetic
beta_q1p <- 
  mag_weighted %>% 
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree, metric = "S", out = "pair")

# Functional KEGG
dist_kegg <-
  genome_kegg_paths %>%
  column_to_rownames(var = "genome") %>%
  traits2dist(., method = "gower")

beta_q1f_kegg <-
  mag_weighted %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist_kegg, metric = "S", out = "pair")

save(beta_q1n, 
     beta_q1p, 
     beta_q1f_gifts,
     file = "data/beta_div.Rdata")

6.2.1 Compare individuals from same trial and pen between sampling days

The resulting stats correspond to Supplementary Table S3.

load("data/beta_div.Rdata")

metadata_2 <-
  chicken_metadata %>%
  dplyr::rename(second = 'animal_code',
                sampling_time_2 = 'sampling_time',
                pen_2 = 'pen',
                trial_2 = 'trial')

# Neutral 
betadiv_n_dis <-
  beta_q1n %>% 
  rename(animal_code = 'first') %>%
  left_join(chicken_metadata, by = 'animal_code') %>%
  left_join(metadata_2, by = 'second') %>%
  mutate(diff = case_when(
    sampling_time == '7' & sampling_time_2 == '21' ~ '7_21',
    sampling_time == '21' & sampling_time_2 == '35' ~ '21_35',
    sampling_time == '7' & sampling_time_2 == '35' ~ '7_35')) %>%
  mutate(diff_trial = case_when(trial == 'CA' & trial_2 == 'CA' ~ 'CA',
                                trial == 'CB' & trial_2 == 'CB' ~ 'CB')) %>%
  mutate(diff_pen = case_when(pen == pen_2  ~ 'same')) %>%
  drop_na(diff, diff_trial, diff_pen)

# Phylogenetic
betadiv_p_dis <-
  beta_q1p %>% 
  rename(animal_code = 'first') %>%
  left_join(chicken_metadata, by = 'animal_code') %>%
  left_join(metadata_2, by = 'second') %>%
  mutate(diff = case_when(
    sampling_time == '7' & sampling_time_2 == '21' ~ '7_21',
    sampling_time == '21' & sampling_time_2 == '35' ~ '21_35',
    sampling_time == '7' & sampling_time_2 == '35' ~ '7_35')) %>%
  mutate(diff_trial = case_when(trial == 'CA' & trial_2 == 'CA' ~ 'CA',
                                trial == 'CB' & trial_2 == 'CB' ~ 'CB')) %>%
  mutate(diff_pen = case_when(pen == pen_2  ~ 'same')) %>%
  drop_na(diff, diff_trial, diff_pen)

# Functional KEGG
betadiv_f_dis_kegg <-
  beta_q1f_kegg %>%
  rename(animal_code = 'first') %>%
  left_join(chicken_metadata, by = 'animal_code') %>%
  left_join(metadata_2, by = 'second') %>%
  mutate(diff = case_when(
    sampling_time == '7' & sampling_time_2 == '21' ~ '7_21',
    sampling_time == '21' & sampling_time_2 == '35' ~ '21_35',
    sampling_time == '7' & sampling_time_2 == '35' ~ '7_35')) %>%
  mutate(diff_trial = case_when(trial == 'CA' & trial_2 == 'CA' ~ 'CA',
                                trial == 'CB' & trial_2 == 'CB' ~ 'CB')) %>%
  mutate(diff_pen = case_when(pen == pen_2  ~ 'same')) %>%
  drop_na(diff, diff_trial, diff_pen)

betadiv_n_dis %>%
  group_by(diff) %>%
  summarise(mean = mean(S), sd = sd(S))

# A tibble: 3 x 3
  diff   mean     sd
  <chr> <dbl>  <dbl>
1 21_35 0.458 0.100 
2 7_21  0.542 0.0907
3 7_35  0.597 0.0917

betadiv_p_dis %>%
  group_by(diff) %>%
  summarise(mean = mean(S), sd = sd(S))

# A tibble: 3 x 3
  diff   mean     sd
  <chr> <dbl>  <dbl>
1 21_35 0.168 0.0776
2 7_21  0.238 0.0883
3 7_35  0.297 0.103

betadiv_f_dis_kegg %>%
  group_by(diff) %>%
  summarise(mean = mean(S), sd = sd(S))

# A tibble: 3 x 3
  diff   mean    sd
  <chr> <dbl> <dbl>
1 21_35 0.205 0.246
2 7_21  0.252 0.291
3 7_35  0.227 0.262

6.2.2 Plot difference between sampling points for beta diversity

The resulting plot corresponds to Supplementary Figure S3.

# Neutral
p_n <- 
  betadiv_n_dis %>%
  filter(!diff == '7_35') %>%
  mutate(diff = factor(diff, levels = c('7_21', '21_35'))) %>%
  ggplot(aes(x = diff, y = S, fill = diff)) +
  geom_boxplot() +
  ylim(0,0.5) +
  theme_minimal() +
  theme(legend.position = 'none') +
  ggtitle('Neutral')

# Phylogenetic
p_p <- 
  betadiv_p_dis %>%
  filter(!diff == '7_35') %>%
  mutate(diff = factor(diff, levels = c('7_21', '21_35'))) %>%
  ggplot(aes(x = diff, y = S, fill = diff)) +
  geom_boxplot() +
  ylim(0,0.5) +
  theme_minimal() +
  theme(legend.position = 'none') +
  ggtitle('Phylogenetic')

# Functional KEGG
p_f_kegg <-
  betadiv_f_dis_kegg %>%
  filter(!diff == '7_35') %>%
  mutate(diff = factor(diff, levels = c('7_21', '21_35'))) %>%
  ggplot(aes(x = diff, y = S, fill = diff)) +
  geom_boxplot() +
  ylim(0,0.5) +
  theme_minimal() +
  theme(legend.position = 'none') +
  ggtitle('Functional')


grid.arrange(p_n, p_p, p_f_kegg, ncol = 3)

6.3 Effect of design in microbiome composition

6.3.1 Permanova values for different metrics of beta diversity

The resulting stats correspond to Supplementary Table S3.

# Dissimilarity tables
perm <- how(nperm = 999)
setBlocks(perm) <- with(chicken_metadata, pen)

# Neutral
b1n_dis_table <- 
  beta_q1n %>%
  bind_rows(beta_q1n %>% rename(first = second, second = first)) %>% 
  pivot_wider(names_from = second, values_from = S) %>%
  column_to_rownames('first') %>%
  mutate(across(everything(), ~replace(., is.na(.), 0))) %>% 
  as.dist()

adonis2(b1n_dis_table ~ sex * age + 
                        breed * age + 
                        treatment * age + 
                        trial * age,
        permutations = perm,
        data = chicken_metadata,
        by = "term")

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  with(chicken_metadata, pen) 
Permutation: free
Number of permutations: 999

adonis2(formula = b1n_dis_table ~ sex * age + breed * age + treatment * age + trial * age, data = chicken_metadata, permutations = perm, by = "term")
               Df SumOfSqs      R2       F Pr(>F)    
sex             1    0.497 0.00950  4.6173  0.001 ***
age             1    4.158 0.07950 38.6408  0.001 ***
breed           1    0.854 0.01632  7.9320  0.001 ***
treatment       2    0.824 0.01575  3.8282  0.001 ***
trial           1    2.919 0.05579 27.1201  0.001 ***
sex:age         1    0.382 0.00731  3.5534  0.001 ***
age:breed       1    0.271 0.00518  2.5184  0.009 ** 
age:treatment   2    0.562 0.01074  2.6100  0.001 ***
age:trial       1    1.379 0.02636 12.8118  0.001 ***
Residual      376   40.464 0.77355                   
Total         387   52.310 1.00000                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Phylogenetic
b1p_dis_table <- 
  beta_q1p %>%
  bind_rows(beta_q1n %>% rename(first = second, second = first)) %>% 
  pivot_wider(names_from = second, values_from = S) %>%
  column_to_rownames('first') %>%
  mutate(across(everything(), ~replace(., is.na(.), 0))) %>% 
  as.dist()

adonis2(b1p_dis_table ~ sex * age + 
                        breed * age + 
                        treatment * age + 
                        trial * age,
        permutations = perm,
        data = chicken_metadata,
        by = "term")

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  with(chicken_metadata, pen) 
Permutation: free
Number of permutations: 999

adonis2(formula = b1p_dis_table ~ sex * age + breed * age + treatment * age + trial * age, data = chicken_metadata, permutations = perm, by = "term")
               Df SumOfSqs      R2       F Pr(>F)    
sex             1    0.497 0.00950  4.6173  0.001 ***
age             1    4.158 0.07950 38.6408  0.001 ***
breed           1    0.854 0.01632  7.9320  0.001 ***
treatment       2    0.824 0.01575  3.8282  0.001 ***
trial           1    2.919 0.05579 27.1201  0.001 ***
sex:age         1    0.382 0.00731  3.5534  0.001 ***
age:breed       1    0.271 0.00518  2.5184  0.010 ** 
age:treatment   2    0.562 0.01074  2.6100  0.001 ***
age:trial       1    1.379 0.02636 12.8118  0.001 ***
Residual      376   40.464 0.77355                   
Total         387   52.310 1.00000                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Functional KEGG
b1f_dis_table_kegg <-
  beta_q1f_kegg %>%
  bind_rows(beta_q1n %>% rename(first = second, second = first)) %>%
  pivot_wider(names_from = second, values_from = S) %>%
  column_to_rownames('first') %>%
  mutate(across(everything(), ~replace(., is.na(.), 0))) %>%
  as.dist()

adonis2(b1f_dis_table_kegg ~ sex * age + 
                             breed * age + 
                             treatment * age + 
                             trial * age,
        permutations = perm,
        data = chicken_metadata,
        by = "term")

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  with(chicken_metadata, pen) 
Permutation: free
Number of permutations: 999

adonis2(formula = b1f_dis_table_kegg ~ sex * age + breed * age + treatment * age + trial * age, data = chicken_metadata, permutations = perm, by = "term")
               Df SumOfSqs      R2       F Pr(>F)    
sex             1    0.497 0.00950  4.6173  0.001 ***
age             1    4.158 0.07950 38.6408  0.001 ***
breed           1    0.854 0.01632  7.9320  0.001 ***
treatment       2    0.824 0.01575  3.8282  0.001 ***
trial           1    2.919 0.05579 27.1201  0.001 ***
sex:age         1    0.382 0.00731  3.5534  0.001 ***
age:breed       1    0.271 0.00518  2.5184  0.005 ** 
age:treatment   2    0.562 0.01074  2.6100  0.001 ***
age:trial       1    1.379 0.02636 12.8118  0.001 ***
Residual      376   40.464 0.77355                   
Total         387   52.310 1.00000                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.3.2 Community composition

perm <- how(nperm = 999)
setBlocks(perm) <- with(chicken_metadata, pen)

t_hel <- 
 hel %>% 
  t() %>% 
  as.data.frame()
  
adonis2(t_hel ~ breed * sampling_time +
                sex * sampling_time +
                treatment * sampling_time +
                trial * sampling_time,
        permutations = perm,
        data = chicken_metadata,
        by = "term")

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  with(chicken_metadata, pen) 
Permutation: free
Number of permutations: 999

adonis2(formula = t_hel ~ breed * sampling_time + sex * sampling_time + treatment * sampling_time + trial * sampling_time, data = chicken_metadata, permutations = perm, by = "term")
                         Df SumOfSqs      R2       F Pr(>F)    
breed                     1   0.3065 0.01190  6.7490  0.001 ***
sampling_time             2   6.2157 0.24127 68.4386  0.001 ***
sex                       1   0.0994 0.00386  2.1882  0.001 ***
treatment                 2   0.1430 0.00555  1.5749  0.001 ***
trial                     1   0.8039 0.03121 17.7039  0.001 ***
breed:sampling_time       2   0.1629 0.00632  1.7933  0.027 *  
sampling_time:sex         2   0.1366 0.00530  1.5036  0.079 .  
sampling_time:treatment   4   0.2162 0.00839  1.1901  0.184    
sampling_time:trial       2   0.8763 0.03402  9.6487  0.001 ***
Residual                370  16.8020 0.65219                   
Total                   387  25.7625 1.00000                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Best fitted model considered for the manuscript                                                                                                   
adonis2(t_hel ~ breed + sex + treatment + trial * sampling_time,
        permutations = perm,
        data = chicken_metadata,
        by = "term")

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  with(chicken_metadata, pen) 
Permutation: free
Number of permutations: 999

adonis2(formula = t_hel ~ breed + sex + treatment + trial * sampling_time, data = chicken_metadata, permutations = perm, by = "term")
                     Df SumOfSqs      R2       F Pr(>F)    
breed                 1   0.3065 0.01190  6.6903  0.001 ***
sex                   1   0.0980 0.00381  2.1404  0.001 ***
treatment             2   0.1552 0.00603  1.6944  0.001 ***
trial                 1   0.7519 0.02919 16.4145  0.001 ***
sampling_time         2   6.2568 0.24287 68.2924  0.001 ***
trial:sampling_time   2   0.8781 0.03408  9.5839  0.001 ***
Residual            378  17.3159 0.67214                   
Total               387  25.7625 1.00000                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.4 Community composition development

6.4.1 Distance-based RDA

The resulting plot corresponds to Figure 2b in the manuscript.

set.seed(4)
hel_bray <- vegdist(t(hel), method = 'bray')
hel_bray_sqrt <- sqrt(hel_bray)

pcoa <- cmdscale(hel_bray_sqrt, k = ncol(hel) - 1, eig = TRUE)
pcoa_scores <- pcoa$points

db_rda <- rda(pcoa_scores ~ sampling_time, data = chicken_metadata)

# Stats
anova(db_rda, by = 'term')

Permutation test for rda under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

Model: rda(formula = pcoa_scores ~ sampling_time, data = chicken_metadata)
               Df Variance      F Pr(>F)    
sampling_time   2 0.022776 27.979  0.001 ***
Residual      385 0.156701                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

RsquareAdj(db_rda)

$r.squared
[1] 0.1269005

$adj.r.squared
[1] 0.1223649

db_rda_scores <-
  data.frame(
    scores(db_rda, display = 'wa'),
    pen = chicken_metadata$pen,
    time = chicken_metadata$sampling_time
    )

db_rda_scores %>%
  ggplot(aes(x = RDA1, y = RDA2, colour = time)) +
  geom_point(size = 3) +
  scale_color_manual(values = c('#e6a024', '#cc6777', '#5bb4e5')) +
  theme_minimal() +
  theme(legend.position = 'none',
        axis.title.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 8),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank())

rm(list = ls())