(a tight drop-in that restricts to Group9/Group10/Group11 and runs stats/plotting for Shannon only)

    ## ---- Subset to Group9/10/11 and keep factor levels tidy ---------------------
    library(dplyr)
    library(ggplot2)
    library(ggpubr)
    library(rstatix)
    library(knitr)
    library(kableExtra)

    shannon_g9_11 <- div.df2 %>%
    filter(Group %in% c("Group9", "Group10", "Group11")) %>%
    mutate(Group = factor(Group, levels = c("Group9","Group10","Group11"))) %>%
    droplevels()

    ## Quick table of the data used
    knitr::kable(shannon_g9_11[, c("Sample name","Group","Shannon")]) %>%
    kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"))

    ## ---- Summary stats (report in text/table) -----------------------------------
    sum_stats <- shannon_g9_11 %>%
    group_by(Group) %>%
    summarise(n = n(),
                mean = mean(Shannon, na.rm = TRUE),
                sd = sd(Shannon, na.rm = TRUE),
                median = median(Shannon, na.rm = TRUE),
                IQR = IQR(Shannon, na.rm = TRUE),
                .groups = "drop")

    write.csv(sum_stats, "Shannon_Group9_11_summary.csv", row.names = FALSE)
    knitr::kable(sum_stats, digits = 3) %>%
    kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"))

    ## ---- Overall test (3 groups) -------------------------------------------------
    # Nonparametric overall test (robust default)
    overall_kw <- compare_means(Shannon ~ Group, data = shannon_g9_11, method = "kruskal.test")
    knitr::kable(overall_kw, digits = 3) %>%
    kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"))
    write.csv(overall_kw, "Shannon_Group9_11_overall_Kruskal.csv", row.names = FALSE)

    ## Optional parametric overall (use if assumptions OK)
    overall_anova <- compare_means(Shannon ~ Group, data = shannon_g9_11, method = "anova")
    write.csv(overall_anova, "Shannon_Group9_11_overall_ANOVA.csv", row.names = FALSE)

    ## Assumption checks (optional; helps decide ANOVA vs KW)
    shapiro_res <- shannon_g9_11 %>% group_by(Group) %>% shapiro_test(Shannon)
    levene_res  <- shannon_g9_11 %>% levene_test(Shannon ~ Group)
    write.csv(shapiro_res, "Shannon_Group9_11_shapiro.csv", row.names = FALSE)
    write.csv(levene_res,  "Shannon_Group9_11_levene.csv",  row.names = FALSE)

    ## ---- Pairwise tests with BH correction --------------------------------------
    # Wilcoxon (nonparametric)
    pw_wilcox <- pairwise_wilcox_test(shannon_g9_11, Shannon ~ Group,
                                    p.adjust.method = "BH", exact = FALSE)
    write.csv(pw_wilcox, "Shannon_Group9_11_pairwise_wilcox_BH.csv", row.names = FALSE)

    # t-tests (parametric, optional)
    pw_t <- pairwise_t_test(shannon_g9_11, Shannon ~ Group, p.adjust.method = "BH")
    write.csv(pw_t, "Shannon_Group9_11_pairwise_t_BH.csv", row.names = FALSE)

    ## ---- Plot: box/jitter with overall & pairwise p-values ----------------------
    my_comparisons <- list(c("Group9","Group10"),
                        c("Group9","Group11"),
                        c("Group10","Group11"))

    p_shannon_g9_11 <- ggplot(shannon_g9_11, aes(x = Group, y = Shannon)) +
    geom_boxplot(outlier.shape = NA) +
    geom_jitter(width = 0.12, alpha = 0.7) +
    labs(y = "Shannon diversity", x = NULL,
        title = "Shannon diversity: Group9 vs Group10 vs Group11") +
    theme_bw() +
    stat_compare_means(method = "kruskal.test", label = "p.format", label.y.npc = "top") +  # overall
    stat_compare_means(comparisons = my_comparisons, method = "wilcox.test",
                        label = "p.signif", hide.ns = TRUE)

    ggsave("Shannon_Group9_10_11_boxplot.pdf", p_shannon_g9_11, width = 5.5, height = 4.2)

    #How to report: cite the Kruskal–Wallis p for the overall 3-group comparison and the BH-adjusted Wilcoxon p-values for pairwise contrasts (use the ANOVA/t-test outputs only if normality and homogeneity look acceptable from the Shapiro/Levene tables).