3  多样性分析

3.1 Alpha 和 Beta 多样性

本章节复现论文 Fig. 2A-F，分析不同处理下土壤微生物群落的 alpha 多样性（Shannon 指数）和 beta 多样性（Bray-Curtis PCoA）。

核心发现：根相关微生物群落不是静态结构，而是随大豆发育阶段明显变化。

library(vegan)
library(Rmisc)
library(rstatix)
library(ggplot2)

treatment <- c('Ctrl', '-N', '-P', '-K')
stage <- c('D0', 'D1', 'D4', 'D7', 'D14', 'D28', 'D42', 'D60', 'D72')

metadata <- read.csv('metadata.csv', row.names = 1)
metadata$Treatment <- factor(metadata$Treatment, levels = treatment)
metadata$Stage <- factor(metadata$Stage, levels = stage)

qmp_data <- read.csv('qmp_data.csv', row.names = 1, check.names = FALSE)

# calculate Shannon diversity
shannon <- data.frame(Shannon = diversity(as.data.frame(t(qmp_data)), index = 'shannon', base = exp(1)))

3.1.1 Alpha 多样性 — 散装土壤 (BS)

Bulk Soil 的 Shannon 指数在各处理间差异较弱（论文标注为 n.s.），说明长期施肥处理对总土壤多样性并非简单线性改变。

group_bs <- metadata[metadata$Compartment == 'BS', ]
shannon_df <- merge(group_bs, shannon, by = 'row.names', all.x = TRUE)
shannon_df_summary <- summarySE(shannon_df, measurevar = 'Shannon', groupvars = 'Treatment')

for (i in 1:4) {
    shannon_df_summary$max <- max(shannon_df[shannon_df$Treatment == shannon_df_summary$Treatment[i], 'Shannon'])
}

dunn_res <- dunn_test(shannon_df, Shannon ~ Treatment, p.adjust.method = 'fdr')
dunn_res_df <- data.frame(dunn_res[-1])
for (n in 1:nrow(dunn_res_df)) {
    if (dunn_res_df$p.adj[n] <= 0.001)
        dunn_res_df$p.adj.signif[n] <- '***'
    else if (dunn_res_df$p.adj[n] <= 0.01)
        dunn_res_df$p.adj.signif[n] <- '**'
    else if (dunn_res_df$p.adj[n] <= 0.05)
        dunn_res_df$p.adj.signif[n] <- '*'
    else
        dunn_res_df$p.adj.signif[n] <- ''
}
Label <- c('', dunn_res_df$p.adj.signif[1:3])

ggplot(shannon_df, aes(x = Treatment, y = Shannon, color = Treatment)) +
    geom_boxplot(width = 0.68, outlier.shape = NA) +
    geom_point(size = 1) +
    geom_text(data = shannon_df_summary, aes(y = max * 1.05, label = Label),
              position = position_dodge(0.9), size = 3) +
    labs(x = '', y = 'Shannon index') +
    scale_y_continuous(limits = c(6, 7.5), labels = c(6.0, 6.5, 7.0, 7.5)) +
    theme_bw() +
    theme(legend.position = 'none')

散装土壤 Shannon 指数

3.1.2 Alpha 多样性 — 根际 (R)

根际 Shannon 指数随时间出现阶段性波动，说明植物发育过程会持续重塑根际细菌多样性。不同处理间的差异反映了不平衡施肥对根际微生物的影响。

group_r <- metadata[metadata$Compartment == 'R', ]
df_r <- merge(group_r, shannon, by = 'row.names', all.x = TRUE)
df_r_summary <- summarySE(df_r, measurevar = 'Shannon', groupvars = c('Day', 'Treatment'))

ggplot(df_r_summary, aes(x = Day, y = Shannon, fill = Treatment, color = Treatment, group = Treatment)) +
    geom_line(linewidth = 0.8) +
    geom_ribbon(aes(ymin = Shannon - sd, ymax = Shannon + sd, fill = Treatment), alpha = 0.1, colour = NA) +
    labs(x = 'Days post germination (d)', y = 'Shannon index') +
    scale_x_continuous(breaks = c(1, 4, 7, 14, 28, 42, 60, 72)) +
    theme_bw() +
    theme(panel.grid.minor = element_blank())

根际 Shannon 指数随时间变化

3.1.3 Alpha 多样性 — 根表 (RS)

根表（Root Endosphere）的 Shannon 指数波动幅度更大，反映根内微生物群落受到更强的宿主筛选和阶段性选择压力。

group_rs <- metadata[metadata$Compartment == 'RS', ]
df_rs <- merge(group_rs, shannon, by = 'row.names', all.x = TRUE)
df_rs_summary <- summarySE(df_rs, measurevar = 'Shannon', groupvars = c('Day', 'Treatment'))

ggplot(df_rs_summary, aes(x = Day, y = Shannon, fill = Treatment, color = Treatment, group = Treatment)) +
    geom_line(linewidth = 0.8) +
    geom_ribbon(aes(ymin = Shannon - sd, ymax = Shannon + sd, fill = Treatment), alpha = 0.1, colour = NA) +
    labs(x = 'Days post germination (d)', y = 'Shannon index') +
    scale_x_continuous(breaks = c(1, 4, 7, 14, 28, 42, 60, 72)) +
    theme_bw() +
    theme(panel.grid.minor = element_blank())

根表 Shannon 指数随时间变化

3.1.4 Beta 多样性 — 散装土壤 (BS)

PCoA 基于 Bray-Curtis 距离，展示样本间群落结构的差异。PERMANOVA 检验 Treatment 的解释度（R²）。

bray_dist <- vegdist(as.data.frame(t(qmp_data)), method = 'bray')
bray_dist <- as.data.frame(as.matrix(bray_dist))

tmp_bd <- bray_dist[row.names(group_bs), row.names(group_bs)]
tmp_dist <- as.dist(tmp_bd, diag = FALSE, upper = FALSE)

tmp_pcoa <- cmdscale(tmp_dist, k = 3, eig = TRUE)
tmp_eig <- round(tmp_pcoa$eig / sum(tmp_pcoa$eig) * 100, 2)

tmp_points <- as.data.frame(tmp_pcoa$points)
names(tmp_points) <- paste0('PCoA', 1:3)
points <- cbind(tmp_points, group_bs)

sig_res <- adonis2(tmp_dist ~ Treatment, data = group_bs)

ggplot(points, aes(x = PCoA1, y = PCoA2, shape = Treatment, color = Treatment)) +
    geom_point(size = 2, alpha = 0.75) +
    stat_ellipse(level = 0.8) +
    labs(subtitle = paste0('Treatment: R2 = ', round(sig_res$R2[1], 3)),
         x = paste0('PCoA1 (', tmp_eig[1], '%)'),
         y = paste0('PCoA2 (', tmp_eig[2], '%)')) +
    scale_shape_manual(values = c(1, 15, 16, 17)) +
    theme_bw()

散装土壤 PCoA 分析

3.1.5 Beta 多样性 — 根际 (R)

根际 PCoA 中不同 Stage 和 Treatment 的分离表明，发育阶段和施肥处理共同影响群落结构。R² 值反映了各因素对群落变异的解释程度。

tmp_bd <- bray_dist[row.names(group_r), row.names(group_r)]
tmp_dist <- as.dist(tmp_bd, diag = FALSE, upper = FALSE)

tmp_pcoa <- cmdscale(tmp_dist, k = 3, eig = TRUE)
tmp_eig <- round(tmp_pcoa$eig / sum(tmp_pcoa$eig) * 100, 2)

tmp_points <- as.data.frame(tmp_pcoa$points)
names(tmp_points) <- paste0('PCoA', 1:3)
points <- cbind(tmp_points, group_r)

sig_res <- adonis2(tmp_dist ~ Stage * Treatment, data = group_r)

ggplot(points, aes(x = PCoA1, y = PCoA2, shape = Treatment, color = Stage)) +
    geom_point(size = 2, alpha = 0.75) +
    labs(subtitle = paste0('Stage: R2 = ', round(sig_res$R2[1], 3), '   Treatment: R2 = ', round(sig_res$R2[2], 3)),
         x = paste0('PCoA1 (', tmp_eig[1], '%)'),
         y = paste0('PCoA2 (', tmp_eig[2], '%)')) +
    scale_shape_manual(values = c(1, 15, 16, 17)) +
    scale_color_manual(values = c('#00512c', '#00a04b', '#84d775', '#c2ebba', '#ffc680', '#ff8c00', '#f34a00', '#e60700')) +
    theme_bw()

根际 PCoA 分析

3.1.6 Beta 多样性 — 根表 (RS)

根表的 PCoA 呈现明显的发育轨迹，说明进入根内的菌群受到更强的宿主筛选和阶段性选择。

tmp_bd <- bray_dist[row.names(group_rs), row.names(group_rs)]
tmp_dist <- as.dist(tmp_bd, diag = FALSE, upper = FALSE)

tmp_pcoa <- cmdscale(tmp_dist, k = 3, eig = TRUE)
tmp_eig <- round(tmp_pcoa$eig / sum(tmp_pcoa$eig) * 100, 2)

tmp_points <- as.data.frame(tmp_pcoa$points)
names(tmp_points) <- paste0('PCoA', 1:3)
points <- cbind(tmp_points, group_rs)

sig_res <- adonis2(tmp_dist ~ Stage * Treatment, data = group_rs)

ggplot(points, aes(x = PCoA1, y = PCoA2, shape = Treatment, color = Stage)) +
    geom_point(size = 2, alpha = 0.75) +
    labs(subtitle = paste0('Stage: R2 = ', round(sig_res$R2[1], 3), '   Treatment: R2 = ', round(sig_res$R2[2], 3)),
         x = paste0('PCoA1 (', tmp_eig[1], '%)'),
         y = paste0('PCoA2 (', tmp_eig[2], '%)')) +
    scale_shape_manual(values = c(1, 15, 16, 17)) +
    scale_color_manual(values = c('#00512c', '#00a04b', '#84d775', '#c2ebba', '#ffc680', '#ff8c00', '#f34a00', '#e60700')) +
    theme_bw()

根表 PCoA 分析

3.1.7 结果解读

生态位	Alpha 多样性	Beta 多样性
BS（散装土壤）	处理间差异弱（n.s.），土壤背景群落稳定	Treatment R² 较低
R（根际）	随时间阶段性波动	Stage 和 Treatment R² 均较高
RS（根表）	波动幅度更大，宿主筛选作用强	明显发育轨迹

核心结论：根相关微生物群落随大豆发育阶段动态变化，不平衡施肥重塑了根际和根内微生物群落结构。