Chapter 5 Beta diversity

5.1 Metagenomics

load("resources/metagenomics/data_filtered.Rdata")
species_color <- c("#e5bd5b", "#6b7398", "#76b183")

Identifying outliers

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if( ~ !all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(
    x = NMDS1,
    y = NMDS2,
    color = Species
  )) +
  geom_point(size = 4) +
  scale_color_manual(values = species_color) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, 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") +
  labs(color = "Species")+geom_text_repel(aes(label = sample), size=3)

load("resources/metagenomics/beta_filtered.Rdata")
sample_metadata <- sample_metadata %>%
  filter(!sample %in% c("H06", "H20"))
genome_counts_filt <- genome_counts_filt %>%
  select(-H06, -H20) %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if( ~ !all(. == 0))

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if( ~ !all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if( ~ !all(. == 0)) %>%
  hillpair(., q = 1)

5.1.1 MAG level

Richness

betadisper(beta_q0n$S, sample_metadata$Species) %>% 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     2 0.007116 0.0035580 0.6172    999  0.545
Residuals 46 0.265159 0.0057643                     

adonis2(beta_q0n$S ~ Species, data = sample_metadata %>% arrange(match(sample, labels(beta_q0n$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.983486	0.09567599	2.433362	0.001
Residual	46	18.747800	0.90432401	NA	NA
Total	48	20.731286	1.00000000	NA	NA

pairwise.adonis(beta_q0n$S, sample_metadata$Species, perm = 999)

     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Eb vs Ha  1 0.9241415 2.303140 0.08435439   0.003      0.009   *
2 Eb vs Pk  1 0.7976225 1.909961 0.05985469   0.015      0.045   .
3 Ha vs Pk  1 1.1978229 2.967307 0.07424335   0.001      0.003   *

Neutral

betadisper(beta_q1n$S, sample_metadata$Species) %>% 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     2 0.05062 0.025309 1.7845    999  0.184
Residuals 46 0.65239 0.014183                     

adonis2(beta_q1n$S ~ Species, data = sample_metadata %>% arrange(match(sample, labels(beta_q1n$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.847557	0.08810564	2.22222	0.001
Residual	46	19.122233	0.91189436	NA	NA
Total	48	20.969790	1.00000000	NA	NA

pairwise.adonis(beta_q1n$S, sample_metadata$Species, perm = 999)

     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Eb vs Ha  1 0.6553806 1.641956 0.06163044   0.100      0.300    
2 Eb vs Pk  1 0.8860212 2.126283 0.06618514   0.010      0.030   .
3 Ha vs Pk  1 1.1405724 2.676921 0.06746797   0.002      0.006   *

5.1.1.1 NMDS

Richness

beta_q0n$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Species = factor(
      Species,
      levels = c("Eb", "Ha", "Pk"),
      labels = c("C. bottae", "H. ariel", "P. kuhlii")
      )
  ) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species)) +
  geom_point(size = 3) +
  scale_color_manual(values=species_color,  labels = c(expression(italic("C. bottae")),
             expression(italic("H. ariel")),
             expression(italic("P. kuhlii"))))+
  geom_segment(
    aes(
      x = x_cen,
      y = y_cen,
      xend = NMDS1,
      yend = NMDS2
    ),
    alpha = 0.9,
    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"
  ) +
  labs(color = "Species")

Neutral

beta_q1n$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(
    x = NMDS1,
    y = NMDS2,
    color = Species
  )) +
  geom_point(size = 3) +
  scale_color_manual(values=species_color,  labels = c(expression(italic("C. bottae")),
             expression(italic("H. ariel")),
             expression(italic("P. kuhlii"))))+
  geom_segment(
    aes(
      x = x_cen,
      y = y_cen,
      xend = NMDS1,
      yend = NMDS2
    ),
    alpha = 0.9,
    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"
  ) +
  labs(color = "Species")

5.1.2 At different taxonomic level

genome_metadata <- genome_metadata %>%
  mutate(
    class  = if_else(str_trim(class) == "", paste0(phylum, "_Unclassified"), class),
    order  = if_else(str_trim(order) == "",
                     if_else(str_detect(class, "Unclassified"), class, paste0(class, "_Unclassified")),
                     order),
    family = if_else(str_trim(family) == "",
                     if_else(str_detect(order, "Unclassified"), order, paste0(order, "_Unclassified")),
                     family),
    genus  = if_else(str_trim(genus) == "",
                     if_else(str_detect(family, "Unclassified"), family, paste0(family, "_Unclassified")),
                     genus)
  )

genome_tax <- genome_counts_filt %>% 
  left_join(genome_metadata[c(1,3,4,6)], by = "genome")

# Aggregate counts at the Family level
family_level_agg <- genome_tax %>%
  select(-genome) %>%
  group_by(phylum, class, family) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)), .groups = "drop")%>% 
  as.data.frame() %>% 
  select(-c(phylum,class))

# Aggregate counts at the Phylum level
phylum_level_agg <- genome_tax %>%
  select(-genome) %>%
  group_by(phylum) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)), .groups = "drop")%>% 
  as.data.frame() 

5.1.2.1 Family

beta_q0n_shot <- family_level_agg %>%
  select(-H09) %>% 
  column_to_rownames(., "family") %>% 
  hillpair(., q = 0)

adonis2(beta_q0n_shot$S ~ Species, data = sample_metadata %>% filter(sample %in% labels(beta_q0n_shot$S)) %>% arrange(match(sample, labels(beta_q0n_shot$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.912276	0.1024998	2.569633	0.001
Residual	45	16.744105	0.8975002	NA	NA
Total	47	18.656381	1.0000000	NA	NA

pairwise.adonis(beta_q0n_shot$S, sample_metadata %>% filter(sample %in% labels(beta_q0n_shot$S)) %>% arrange(match(sample, labels(beta_q0n_shot$S))) %>% pull(Species), perm = 999)

     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Eb vs Ha  1 0.9787096 2.719914 0.10179350   0.012      0.036   .
2 Eb vs Pk  1 0.7629461 2.012157 0.06285604   0.019      0.057    
3 Ha vs Pk  1 1.1033819 2.947326 0.07567466   0.001      0.003   *

beta_q1n_shot <- family_level_agg %>%
  select(-H19) %>% 
  column_to_rownames(., "family") %>% 
  hillpair(., q = 1)

adonis2(beta_q1n_shot$S ~ Species, data = sample_metadata %>% filter(sample %in% labels(beta_q1n_shot$S)) %>% arrange(match(sample, labels(beta_q1n_shot$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.578943	0.0843342	2.072284	0.01
Residual	45	17.143513	0.9156658	NA	NA
Total	47	18.722456	1.0000000	NA	NA

pairwise.adonis(beta_q1n_shot$S, sample_metadata %>% filter(sample %in% labels(beta_q1n_shot$S)) %>% arrange(match(sample, labels(beta_q1n_shot$S))) %>% pull(Species), perm = 999)

     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Eb vs Ha  1 0.5845919 1.605450 0.06269955   0.118      0.354    
2 Eb vs Pk  1 0.7388943 1.957133 0.06124244   0.052      0.156    
3 Ha vs Pk  1 0.9722440 2.461076 0.06398874   0.014      0.042   .

beta_q1n_shot$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Species = factor(
      Species,
      levels = c("Eb", "Pk", "Ha"),
      labels = c("C. bottae", "P. kuhlii", "H. ariel")
      )) %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species)) +
  geom_point(size = 3) +
  scale_color_manual(values=species_color,  labels = c(expression(italic("C. bottae")),
             expression(italic("H. ariel")),
             expression(italic("P. kuhlii")))) +  
  geom_segment(
    aes(
      x = x_cen,
      y = y_cen,
      xend = NMDS1,
      yend = NMDS2
    ),
    alpha = 0.9,
    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"
  ) +
  labs(color = "Species")

5.1.2.2 Phylum

beta_q0n_P_shot <- phylum_level_agg %>%
  select(-H09) %>% 
  column_to_rownames(., "phylum") %>% 
  hillpair(., q = 0)

adonis2(
  beta_q0n_P_shot$S ~ Species,
  by = "terms",
  data = sample_metadata %>%
    filter(sample %in% labels(beta_q0n_P_shot$S)) %>%
    arrange(match(sample, labels(beta_q0n_P_shot$S))),
  permutations = 999
) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	3.016295	0.3400081	11.59133	0.001
Residual	45	5.854949	0.6599919	NA	NA
Total	47	8.871244	1.0000000	NA	NA

pairwise.adonis(
  beta_q0n_P_shot$S,
  sample_metadata %>%
    filter(sample %in% labels(beta_q0n_P_shot$S)) %>%
    arrange(match(sample, labels(beta_q0n_P_shot$S))) %>%
    pull(Species),
  perm = 999
)

     pairs Df SumsOfSqs   F.Model        R2 p.value p.adjusted sig
1 Eb vs Ha  1 1.7879482 16.299998 0.4044665   0.001      0.003   *
2 Eb vs Pk  1 0.6958001  4.160066 0.1217815   0.006      0.018   .
3 Ha vs Pk  1 2.0007870 17.742599 0.3301403   0.001      0.003   *

beta_q1n_P_shot <- phylum_level_agg %>%
  select(-H19) %>% 
  column_to_rownames(., "phylum") %>% 
  hillpair(., q = 1)

adonis2(beta_q1n_P_shot$S ~ Species, data = sample_metadata %>% filter(sample %in% labels(beta_q1n_P_shot$S)) %>% arrange(match(sample, labels(beta_q1n_P_shot$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.144642	0.1882421	5.217624	0.011
Residual	45	4.936050	0.8117579	NA	NA
Total	47	6.080692	1.0000000	NA	NA

pairwise.adonis(beta_q1n_P_shot$S, sample_metadata %>% filter(sample %in% labels(beta_q1n_P_shot$S)) %>% arrange(match(sample, labels(beta_q1n_P_shot$S))) %>% pull(Species), perm = 999)

     pairs Df SumsOfSqs   F.Model        R2 p.value p.adjusted sig
1 Eb vs Ha  1 0.3345227  4.966266 0.1714500   0.014      0.042   .
2 Eb vs Pk  1 0.1592735  1.095082 0.0352172   0.371      1.000    
3 Ha vs Pk  1 1.0828550 10.015742 0.2176590   0.002      0.006   *

5.2 Amplicon

load("resources/amplicon/data_standard_filt.Rdata")

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if( ~ !all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if( ~ !all(. == 0)) %>%
  hillpair(., q = 1)

5.2.1 ASV level

Richness

betadisper(beta_q0n$S, sample_metadata %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample, labels(beta_q0n$S))) %>% pull(Species)) %>%
  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     2 0.004486 0.00224298 4.2258    999  0.021 *
Residuals 48 0.025478 0.00053078                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(beta_q0n$S ~ Species,  data = sample_metadata %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample, labels(beta_q0n$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.441876	0.06132966	1.568082	0.001
Residual	48	22.068375	0.93867034	NA	NA
Total	50	23.510250	1.00000000	NA	NA

pairwise.adonis(beta_q0n$S, sample_metadata %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample, labels(beta_q0n$S))) %>% pull(Species), perm = 999)

     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Eb vs Ha  1 0.8054324 1.770936 0.06155296   0.001      0.003   *
2 Eb vs Pk  1 0.6163397 1.329030 0.04242169   0.017      0.051    
3 Ha vs Pk  1 0.7428110 1.614406 0.03974960   0.001      0.003   *

Neutral

betadisper(beta_q1n$S, sample_metadata %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample, labels(beta_q0n$S))) %>% pull(Species)) %>%
  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     2 0.018083 0.0090413 5.5707    999  0.005 **
Residuals 48 0.077904 0.0016230                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(beta_q1n$S ~ Species, data = sample_metadata %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample, labels(beta_q0n$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.552653	0.06761273	1.740377	0.002
Residual	48	21.411254	0.93238727	NA	NA
Total	50	22.963907	1.00000000	NA	NA

pairwise.adonis(beta_q1n$S, sample_metadata %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample, labels(beta_q0n$S))) %>% pull(Species), perm = 999) %>% 
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Eb vs Ha	1	0.9795980	2.278986	0.07783691	0.002	0.006	*
Eb vs Pk	1	0.6422682	1.410368	0.04490136	0.069	0.207
Ha vs Pk	1	0.7371794	1.637698	0.04029996	0.002	0.006	*

5.2.1.1 NMDS

Richness

beta_q0n$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species)) +
 geom_point(size = 3) +
  scale_color_manual(values=species_color,  labels = c(expression(italic("C. bottae")),
             expression(italic("H. ariel")),
             expression(italic("P. kuhlii"))))+
  geom_segment(
    aes(
      x = x_cen,
      y = y_cen,
      xend = NMDS1,
      yend = NMDS2
    ),
    alpha = 0.9,
    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"
  ) +
  labs(color = "Species")

Neutral

beta_q1n$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species)) +
   geom_point(size = 3) +
  scale_color_manual(values=species_color,  labels = c(expression(italic("C. bottae")),
             expression(italic("H. ariel")),
             expression(italic("P. kuhlii"))))+
  geom_segment(
    aes(
      x = x_cen,
      y = y_cen,
      xend = NMDS1,
      yend = NMDS2
    ),
    alpha = 0.9,
    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"
  ) +
  labs(color = "Species")

#ggsave("Figures_ms/beta_neutral_metabarcoding.pdf", plot = p, width = 8, height = 4)

5.2.2 At different taxonomic level

sample_metadata_amplicon <- sample_metadata[c(1,9,10,13)] %>%
  mutate(sample = paste0(sample, "_ampli")) %>%
  mutate(method="Metabarcoding_standard")

genome_metadata_ampli <- genome_metadata %>%
  mutate(
    class  = if_else(is.na(class), paste0(phylum, "_Unclassified"), class),
    order  = if_else(is.na(order),
                     if_else(str_detect(class, "Unclassified"), class, paste0(class, "_Unclassified")),
                     order),
    family = if_else(is.na(family),
                     if_else(str_detect(order, "Unclassified"), order, paste0(order, "_Unclassified")),
                     family),
    genus  = if_else(is.na(genus),
                     if_else(str_detect(family, "Unclassified"), family, paste0(family, "_Unclassified")),
                     genus)
  )

genome_tax_ampli <- genome_counts_filt %>% 
  column_to_rownames(., "genome")%>%
  rename_with(~ paste0(., "_ampli")) %>%
  rownames_to_column(., "genome")%>%
  left_join(genome_metadata_ampli, by = "genome")

family_level_agg_ampli <- genome_tax_ampli %>%
  select(-genome) %>%
  group_by(phylum, class, family) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)), .groups = "drop")%>% 
    as.data.frame() %>% 
  select(-c(phylum,class))

phylum_level_agg_ampli <- genome_tax_ampli %>%
  select(-genome) %>%
  group_by(phylum) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)), .groups = "drop")%>% 
    as.data.frame()

5.2.2.1 Family level

beta_q0n_ampli <- family_level_agg_ampli %>%
  column_to_rownames(., "family") %>% 
  hillpair(., q = 0)

adonis2(beta_q0n_ampli$S ~ Species, data = all_samples %>% filter(sample %in% labels(beta_q0n_ampli$S)) %>% arrange(match(sample, labels(beta_q0n_ampli$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.495085	0.09224457	2.438839	0.001
Residual	48	14.712753	0.90775543	NA	NA
Total	50	16.207838	1.00000000	NA	NA

pairwise.adonis(beta_q0n_ampli$S, all_samples %>% filter(sample %in% labels(beta_q0n_ampli$S)) %>% arrange(match(sample, labels(beta_q0n_ampli$S))) %>% pull(Species),
  perm = 999) %>% 
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Eb vs Ha	1	1.0497877	3.756659	0.12214132	0.001	0.003	*
Eb vs Pk	1	0.5165122	1.558997	0.04939946	0.052	0.156
Ha vs Pk	1	0.7140756	2.332190	0.05642551	0.001	0.003	*

beta_q1n_ampli <- family_level_agg_ampli %>%
  column_to_rownames(., "family") %>% 
  hillpair(., q = 1)

adonis2(beta_q1n_ampli$S ~ Species,  data = all_samples %>% filter(sample %in% labels(beta_q1n_ampli$S)) %>% arrange(match(sample, labels(beta_q1n_ampli$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.508153	0.08460994	2.218332	0.004
Residual	48	16.316616	0.91539006	NA	NA
Total	50	17.824768	1.00000000	NA	NA

pairwise.adonis(beta_q1n_ampli$S, all_samples %>% filter(sample %in% labels(beta_q1n_ampli$S)) %>% arrange(match(sample, labels(beta_q1n_ampli$S))) %>% pull(Species),
  perm = 999) %>% 
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Eb vs Ha	1	1.2490797	4.112390	0.13217854	0.001	0.003	*
Eb vs Pk	1	0.5872306	1.628530	0.05148927	0.077	0.231
Ha vs Pk	1	0.5341731	1.530168	0.03775379	0.085	0.255

beta_q1n_ampli$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(all_samples, by = join_by(sample == sample)) %>%
  group_by(Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Species = factor(
      Species,
      levels = c("Eb", "Pk", "Ha"),
      labels = c("C. bottae", "P. kuhlii", "H. ariel")
      )) %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species)) +
  geom_point(size = 3) +
  scale_color_manual(values=species_color,  labels = c(expression(italic("C. bottae")),
             expression(italic("H. ariel")),
             expression(italic("P. kuhlii")))) +  
  geom_segment(
    aes(
      x = x_cen,
      y = y_cen,
      xend = NMDS1,
      yend = NMDS2
    ),
    alpha = 0.9,
    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"
  ) +
  labs(color = "Species")

5.2.2.2 Phylum level

beta_q0n_P_ampli <- phylum_level_agg_ampli %>%
  column_to_rownames(., "phylum") %>% 
  hillpair(., q = 0)

adonis2(
  beta_q0n_P_ampli$S ~ Species,
  by = "terms",
  data = all_samples %>%
    filter(sample %in% labels(beta_q0n_P_ampli$S)) %>%
    arrange(match(sample, labels(beta_q0n_P_ampli$S))),
  permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	0.8982446	0.1030059	2.756028	0.003
Residual	48	7.8220794	0.8969941	NA	NA
Total	50	8.7203240	1.0000000	NA	NA

pairwise.adonis(
  beta_q0n_P_ampli$S,
  all_samples %>%
    filter(sample %in% labels(beta_q0n_P_ampli$S)) %>%
    arrange(match(sample, labels(beta_q0n_P_ampli$S))) %>%
    pull(Species),
  perm = 999) %>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Eb vs Ha	1	0.6086184	4.449912	0.14149204	0.001	0.003	*
Eb vs Pk	1	0.4967691	2.856505	0.08693879	0.020	0.060
Ha vs Pk	1	0.2991202	1.732330	0.04252960	0.099	0.297

beta_q1n_P_ampli <- phylum_level_agg_ampli %>%
  column_to_rownames(., "phylum") %>% 
  hillpair(., q = 1)

adonis2(beta_q1n_P_ampli$S ~ Species, data = all_samples %>% filter(sample %in% labels(beta_q1n_P_ampli$S)) %>% arrange(match(sample, labels(beta_q1n_P_ampli$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	0.1959103	0.06600639	1.696107	0.19
Residual	48	2.7721407	0.93399361	NA	NA
Total	50	2.9680510	1.00000000	NA	NA

pairwise.adonis(beta_q1n_ampli$S, all_samples %>% filter(sample %in% labels(beta_q1n_P_ampli$S)) %>% arrange(match(sample, labels(beta_q1n_P_ampli$S))) %>% pull(Species), perm = 999)  %>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Eb vs Ha	1	1.2490797	4.112390	0.13217854	0.001	0.003	*
Eb vs Pk	1	0.5872306	1.628530	0.05148927	0.077	0.231
Ha vs Pk	1	0.5341731	1.530168	0.03775379	0.105	0.315

5.3 Comparison

Aggregate counts at family or phylum level in both datasets

load("resources/amplicon/data_standard_filt.Rdata")

sample_metadata_amplicon <- sample_metadata[c(1,9,10,13)] %>%
  mutate(sample = paste0(sample, "_ampli")) %>%
  mutate(method="Metabarcoding_standard")

genome_metadata_ampli <- genome_metadata %>%
  mutate(
    class  = if_else(is.na(class), paste0(phylum, "_Unclassified"), class),
    order  = if_else(is.na(order),
                     if_else(str_detect(class, "Unclassified"), class, paste0(class, "_Unclassified")),
                     order),
    family = if_else(is.na(family),
                     if_else(str_detect(order, "Unclassified"), order, paste0(order, "_Unclassified")),
                     family),
    genus  = if_else(is.na(genus),
                     if_else(str_detect(family, "Unclassified"), family, paste0(family, "_Unclassified")),
                     genus)
  )

genome_tax_ampli <- genome_counts_filt %>% 
  column_to_rownames(., "genome")%>%
  rename_with(~ paste0(., "_ampli")) %>%
  rownames_to_column(., "genome")%>%
  left_join(genome_metadata_ampli, by = "genome")

family_level_agg_ampli <- genome_tax_ampli %>%
  select(-genome) %>%
  group_by(phylum, class, family) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)), .groups = "drop")%>% 
    as.data.frame() %>% 
  select(-c(phylum,class))

phylum_level_agg_ampli <- genome_tax_ampli %>%
  select(-genome) %>%
  group_by(phylum) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)), .groups = "drop")%>% 
    as.data.frame()

load("resources/metagenomics/data_filtered.Rdata")

genome_metadata <- genome_metadata %>%
  mutate(
    class  = if_else(str_trim(class) == "", paste0(phylum, "_Unclassified"), class),
    order  = if_else(str_trim(order) == "",
                     if_else(str_detect(class, "Unclassified"), class, paste0(class, "_Unclassified")),
                     order),
    family = if_else(str_trim(family) == "",
                     if_else(str_detect(order, "Unclassified"), order, paste0(order, "_Unclassified")),
                     family),
    genus  = if_else(str_trim(genus) == "",
                     if_else(str_detect(family, "Unclassified"), family, paste0(family, "_Unclassified")),
                     genus)
  )

genome_tax <- genome_counts_filt %>% 
  left_join(genome_metadata[c(1,3,4,6)], by = "genome")

# Aggregate counts at the Family level
family_level_agg <- genome_tax %>%
  select(-genome) %>%
  group_by(phylum, class, family) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)), .groups = "drop")%>% 
  as.data.frame() %>% 
  select(-c(phylum,class))

# Aggregate counts at the Phylum level
phylum_level_agg <- genome_tax %>%
  select(-genome) %>%
  group_by(phylum) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)), .groups = "drop")%>% 
  as.data.frame() 

merged_table <- full_join(family_level_agg_ampli,family_level_agg, 
          by = c("family"))%>%
  mutate(across(where(is.numeric), ~replace_na(., 0)))

merged_table_phylum <- full_join(phylum_level_agg_ampli,phylum_level_agg, 
          by = c("phylum"))%>%
  mutate(across(where(is.numeric), ~replace_na(., 0)))


sample_shot <- sample_metadata %>% 
  select(-c(Habitat,Longitude, Latitude, Region))%>%
  mutate(method="Metagenomics")

all_samples <- rbind(sample_shot, sample_metadata_amplicon)

beta_q0n <- merged_table %>%
  select(-H19, -H20) %>% 
  column_to_rownames(., "family") %>% 
  hillpair(., q = 0)

beta_q1n <- merged_table %>%
  select(-H19) %>% 
  column_to_rownames(., "family") %>% 
  hillpair(., q = 1)

beta_q0n_phylum <- merged_table_phylum %>%
  select(-H19, -H20) %>% 
  column_to_rownames(., "phylum") %>% 
  hillpair(., q = 0)

beta_q1n_phylum <- merged_table_phylum %>%
  select(-H19) %>% 
  column_to_rownames(., "phylum") %>% 
  hillpair(., q = 1)

5.3.1 Family

Richness

beta_q0n$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(all_samples, by = join_by(sample == sample)) %>%
  group_by(method, Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Species = factor(
      Species,
      levels = c("Eb", "Pk", "Ha"),
      labels = c("C. bottae", "P. kuhlii", "H. ariel")
      )
  ) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species, shape=method)) +
  geom_point(size = 4) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2 ), alpha = 0.9, show.legend = FALSE) +
  scale_color_manual(values = species_color,
                     labels = c(
                       expression(italic("C. bottae")),
                       expression(italic("H. ariel")),
                       expression(italic("P. kuhlii"))
                     ))+
  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"
  ) +
  labs(
    x = "Methodological approach",
    y = "Alpha diversity at phylum level",
    shape = "Methodological approach"
  )#+labs(color = "method") +geom_text_repel(aes(label = sample), size=3)

betadisper(beta_q0n$S, all_samples %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample, labels(beta_q0n$S))) %>% pull(Species)) %>% 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     2 0.02434 0.0121708 3.3657    999  0.034 *
Residuals 97 0.35077 0.0036162                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

all_samples <- all_samples %>% 
  mutate(id=sample) %>% 
  mutate(id = str_remove(id, "_ampli$"))

adonis2(
  beta_q0n$S ~ Species*method,
  data = all_samples %>%
    filter(sample %in% labels(beta_q0n$S)) %>%
    arrange(match(sample, labels(beta_q0n$S))),
  permutations = 999, 
  strata = all_samples %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample,labels(beta_q0n$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	8.166032	0.2049572	4.846525	0.001
Residual	94	31.676598	0.7950428	NA	NA
Total	99	39.842630	1.0000000	NA	NA

adonis2(
  beta_q0n$S ~ Species*method,
  by = "term",
  data = all_samples %>%
    filter(sample %in% labels(beta_q0n$S)) %>%
    arrange(match(sample, labels(beta_q0n$S))),
  permutations = 999, 
  strata = all_samples %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample,labels(beta_q0n$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	2.129351	0.05344404	3.159414	0.001
method	1	4.660736	0.11697862	13.830689	0.001
Species:method	2	1.375945	0.03453450	2.041552	0.009
Residual	94	31.676598	0.79504284	NA	NA
Total	99	39.842630	1.00000000	NA	NA

adonis2(
  beta_q0n$S ~ Species+method,
  by = "margin",
  data = all_samples %>%
    filter(sample %in% labels(beta_q0n$S)) %>%
    arrange(match(sample, labels(beta_q0n$S))),
  permutations = 999, 
  strata = all_samples %>% filter(sample %in% labels(beta_q0n$S)) %>% arrange(match(sample,labels(beta_q0n$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	2.122636	0.05327551	3.082563	0.001
method	1	4.660736	0.11697862	13.536951	0.001
Residual	96	33.052543	0.82957734	NA	NA
Total	99	39.842630	1.00000000	NA	NA

Neutral

beta_q1n$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(all_samples, by = join_by(sample == sample)) %>%
  group_by(method, Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Species = factor(
      Species,
      levels = c("Eb", "Pk", "Ha"),
      labels = c("C. bottae", "P. kuhlii", "H. ariel")
      )
  ) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species, shape=method)) +
  geom_point(size = 3) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2 ), alpha = 0.9, show.legend = FALSE) +
  scale_color_manual(values = species_color,
                     labels = c(
                       expression(italic("C. bottae")),
                       expression(italic("H. ariel")),
                       expression(italic("P. kuhlii"))
                     ))+
  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"
  ) +
  labs(shape = "Methodological approach")

betadisper(
  beta_q1n$S,
  all_samples %>%
    filter(sample %in% labels(beta_q1n$S)) %>%
    arrange(match(sample, labels(beta_q1n$S))) %>%
    pull(method)
) %>%
  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.03764 0.037645 6.3219    999  0.015 *
Residuals 99 0.58951 0.005955                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(
  beta_q1n$S ~ Species*method,
  data = all_samples %>% filter(sample %in% labels(beta_q1n$S)) %>% arrange(match(sample, labels(beta_q1n$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q1n$S)) %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	5.39140	0.1353951	2.975356	0.001
Residual	95	34.42835	0.8646049	NA	NA
Total	100	39.81975	1.0000000	NA	NA

adonis2(
  beta_q1n$S ~ Species*method,
  by = "terms",
  data = all_samples %>% filter(sample %in% labels(beta_q1n$S)) %>% arrange(match(sample, labels(beta_q1n$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q1n$S)) %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	2.4292074	0.06100509	3.351522	0.001
method	1	2.1955523	0.05513727	6.058306	0.001
Species:method	2	0.7666401	0.01925276	1.057716	0.017
Residual	95	34.4283476	0.86460487	NA	NA
Total	100	39.8197474	1.00000000	NA	NA

adonis2(
  beta_q1n$S ~ Species+method,
  by = "margin",
  data = all_samples %>% filter(sample %in% labels(beta_q1n$S)) %>% arrange(match(sample, labels(beta_q1n$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q1n$S)) %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	2.419210	0.06075403	3.333761	0.001
method	1	2.195552	0.05513727	6.051105	0.001
Residual	97	35.194988	0.88385763	NA	NA
Total	100	39.819747	1.00000000	NA	NA

5.3.2 Phylum

Richness

beta_q0n_phylum$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(all_samples, by = join_by(sample == sample)) %>%
  group_by(method, Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species, shape=method)) +
  geom_point(size = 4) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2 ), alpha = 0.9, 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"
  ) +labs(color = "method") + geom_text_repel(aes(label = sample), size=3)

betadisper(beta_q0n_phylum$S, all_samples %>% filter(sample %in% labels(beta_q0n_phylum$S)) %>% arrange(match(sample, labels(beta_q0n_phylum$S))) %>% pull(Species)) %>% 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     2 0.0752 0.037598 2.5521    999  0.082 .
Residuals 97 1.4290 0.014732                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(beta_q0n_phylum$S ~ Species*method, 
        data = all_samples %>% filter(sample %in% labels(beta_q0n_phylum$S)) %>% arrange(match(sample, labels(beta_q0n_phylum$S))), permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	8.982544	0.3963025	12.34142	0.001
Residual	94	13.683335	0.6036975	NA	NA
Total	99	22.665879	1.0000000	NA	NA

adonis2(
  beta_q0n_phylum$S ~ Species*method,
  data = all_samples %>%
    filter(sample %in% labels(beta_q0n_phylum$S)) %>%
    arrange(match(sample, labels(beta_q0n_phylum$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q0n_phylum$S)) %>% arrange(match(sample,labels(beta_q0n_phylum$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	8.982544	0.3963025	12.34142	0.001
Residual	94	13.683335	0.6036975	NA	NA
Total	99	22.665879	1.0000000	NA	NA

adonis2(
  beta_q0n_phylum$S ~ Species*method,
  by = "terms",
  data = all_samples %>%
    filter(sample %in% labels(beta_q0n_phylum$S)) %>%
    arrange(match(sample, labels(beta_q0n_phylum$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q0n_phylum$S)) %>% arrange(match(sample,labels(beta_q0n_phylum$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	2.343946	0.10341297	8.051068	0.001
method	1	5.002004	0.22068430	34.362115	0.001
Species:method	2	1.636595	0.07220521	5.621432	0.004
Residual	94	13.683335	0.60369752	NA	NA
Total	99	22.665879	1.00000000	NA	NA

adonis2(
  beta_q0n_phylum$S ~ Species+method,
  by = "margin",
  data = all_samples %>%
    filter(sample %in% labels(beta_q0n_phylum$S)) %>%
    arrange(match(sample, labels(beta_q0n_phylum$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q0n_phylum$S)) %>% arrange(match(sample,labels(beta_q0n_phylum$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	2.408699	0.1062698	7.546872	0.001
method	1	5.002004	0.2206843	31.344292	0.001
Residual	96	15.319929	0.6759027	NA	NA
Total	99	22.665879	1.0000000	NA	NA

Neutral

beta_q1n_phylum$S %>%
  vegan::metaMDS(.,
                 trymax = 500,
                 k = 2,
                 trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(all_samples, by = join_by(sample == sample)) %>%
  group_by(method, Species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Species = factor(
      Species,
      levels = c("Eb", "Pk", "Ha"),
      labels = c("C. bottae", "P. kuhlii", "H. ariel")
      )
  ) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2, color = Species, shape=method)) +
  geom_point(size = 3) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2 ), alpha = 0.9, show.legend = FALSE) +
  scale_color_manual(values = species_color,
                     labels = c(
                       expression(italic("C. bottae")),
                       expression(italic("H. ariel")),
                       expression(italic("P. kuhlii"))
                     ))+
  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"
  ) +
  labs(shape = "Methodological approach")

betadisper(
  beta_q1n_phylum$S,
  all_samples %>%
    filter(sample %in% labels(beta_q1n_phylum$S)) %>%
    arrange(match(sample, labels(beta_q1n_phylum$S))) %>%
    pull(Species)
) %>%
  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     2 0.2208 0.110381 2.3734    999   0.11
Residuals 98 4.5577 0.046507                     

adonis2(
  beta_q1n_phylum$S ~ Species*method,
  data = all_samples %>%
    filter(sample %in% labels(beta_q1n_phylum$S)) %>%
    arrange(match(sample, labels(beta_q1n_phylum$S))),
  permutations = 999) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	1.388676	0.1386858	3.059314	0.003
Residual	95	8.624432	0.8613142	NA	NA
Total	100	10.013109	1.0000000	NA	NA

adonis2(
  beta_q1n_phylum$S ~ Species*method,
  data = all_samples %>%
    filter(sample %in% labels(beta_q1n_phylum$S)) %>%
    arrange(match(sample, labels(beta_q1n_phylum$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q1n_phylum$S)) %>% arrange(match(sample,labels(beta_q1n_phylum$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	1.388676	0.1386858	3.059314	0.035
Residual	95	8.624432	0.8613142	NA	NA
Total	100	10.013109	1.0000000	NA	NA

adonis2(
  beta_q1n_phylum$S ~ Species*method,
  by = "terms",
  data = all_samples %>%
    filter(sample %in% labels(beta_q1n_phylum$S)) %>%
    arrange(match(sample, labels(beta_q1n_phylum$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q1n_phylum$S)) %>% arrange(match(sample,labels(beta_q1n_phylum$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	1.0276371	0.10262917	5.659823	0.040
method	1	0.1674935	0.01672742	1.844977	0.065
Species:method	2	0.1935457	0.01932923	1.065974	0.145
Residual	95	8.6244323	0.86131417	NA	NA
Total	100	10.0131086	1.00000000	NA	NA

adonis2(
  beta_q1n_phylum$S ~ Species+method,
  by = "margin",
  data = all_samples %>%
    filter(sample %in% labels(beta_q1n_phylum$S)) %>%
    arrange(match(sample, labels(beta_q1n_phylum$S))),
  permutations = 999,
  strata = all_samples %>% filter(sample %in% labels(beta_q1n_phylum$S)) %>% arrange(match(sample,labels(beta_q1n_phylum$S))) %>% pull(id)) %>%
  broom::tidy() %>%
  tt()

term	df	SumOfSqs	R2	statistic	p.value
Species	2	1.0378110	0.10364524	5.708093	0.035
method	1	0.1674935	0.01672742	1.842471	0.070
Residual	97	8.8179780	0.88064340	NA	NA
Total	100	10.0131086	1.00000000	NA	NA