Microbiome/关键类群的特殊代谢功能维持土壤微生物组稳定性

题目：Specialized metabolic functions of keystone taxa sustain soil microbiome stability

关键类群的特殊代谢功能维持土壤微生物组稳定性

期刊：Microbiome

三年均IF：10.402

上线时间：2021.2

通讯作者:张瑞福

摘要

土壤微生物组的生物多样性与生态系统稳定性之间的关系还不清楚。在这里 我们研究了土壤微生物组细菌系统发育多样性对功能性状与稳定性的影响。通过将连续稀释的土壤悬液接种到无菌土壤中，产生不同系统发育多样性的群落，并通过检测不同pH下的群落变化来评估微生物群的稳定性。通过DNA测序来检测分类学特征与功能性状。

我们发现高的细菌系统发育多样性群落更稳定，这表明高微生物多样性对扰动有更高的抗性。通过功能基因共发生网络与机器学习鉴定出特定的代谢功能，一些关键基因如氮代谢、磷酸盐与次磷酸盐代谢。分类注释发现关键功能是由特定的细菌分类群执行的，包括Nitrospira和Gemmatimonas等。

这项研究为我们理解土壤微生物生物多样性和生态系统稳定性之间的关系提供了新的见解，强调关键类群中的特殊代谢功能对土壤微生物稳定性至关重要的作用。

实验设计与微生物组稳定性

研究微生物群落稳定性大多通过分析推断（Microbiome/稻田土壤产甲烷菌的共存模式与（甲烷的产生和群落构建）密切相关，ISME/环境胁迫破坏微生物网络Environmental Microbiology/稀有类群维持作物真菌组的稳定性与生态系统功能）与实证证明（互作强度决定群落的生物多样性和稳定性），通过研究分析微生物共发生网络的特性以及计算恢复性与抵抗性指数。本文通过计算平均物种变异程度AVD来衡量微生物群落稳定性强弱。实验设计在上面一篇NC（NC：南农沈其荣、张瑞福团队揭示多样性激发的确定性细菌装配过程限制群落功能）不作赘述。在计算AVD之前将每个样本稀释至11020个reads，AVD计算公式：

ai表示一个otu的变异程度，xi表示一个样本otu稀释的丰度，_xi表示一个样本otu平均稀释的丰度，σi表示一个样本otu稀释的丰度的偏差。k是一个样本组的样本数，n是每个样本组中OTU的数量。

下面利用机器学习（Decision tree决策树、boosting 提升方法、bagging装袋算法、Nearest neighbor algorithm近邻算法、Support vector machine 支持向量机、Random forest 随机森林、Artificial neural network人工神经网）检验不同微生物功能类别在构建生态系统功能中的重要性。发现随机森林方法精度高，错误率低（一般群落分析多种机器学习方法中随机森林相对比较优良，例如ISME：南农沈其荣团队基于大数据准确预测土壤的枯萎病发生、Microbiome/稻种的驯化在生态进化上塑造水稻种子的细菌与真菌群落）。文章补充材料附上了各种方法的R代码，没有具体跑：

# 数据准备与处理

Fold=function(Z=10,Test,Dv,seed=1000){

n=nrow(Test)

d=1:n;dd=list()

e=levels(Test[,Dv])

T=length(e);set.seed(seed) # T is the dependent variable

for(i in 1:T){

d0=d[w[,Dv]==e[i]];j=length(d0)

ZT=rep(1:Z,ceiling(j/Z))[1:j]

id=cbind(sample(ZT,length(ZT)),d0);dd[[i]]=id}

mm=list();for(i in 1:Z){u=NULL;

for(j in 1:T)u=c(u,dd[[j]][dd[[j]][,1]==i,2])

mm[[i]]=u} # mm[[i]] is the ith subscript set: i=1, 2, …, Z

return(mm)} # Output 10 subscript sets

Test=read.csv(“Test.csv”)

Dv;Z=10;n=nrow(Test);mm=Fold(Z=10,Test,Dv,seed=7777) #Dv is the position of independent

variable

#1Decision tree

library(rpart)

E=rep(0,Z)

for(i in 1:Z){m=mm[[i]]; # mm[[i]] is the ith subscript set extracted by Fold() function: i=1, 2, …, Z

n1=length(m);a=rpart(avd~.,Test[-m,]) # avd is the independent variable

E[i]=sum(Test[m,Dv]!=predict(a,Test[m,])$avd)/n1}

mean(E)

#2boosting

library(adabag)

E=rep(0,Z)

for(i in 1:Z){m=mm[[i]]; # mm[[i]] is the ith subscript set extracted by Fold() function: i=1, 2, …, Z

n1=length(m)

a=boosting(avd~.,Test[-m,]) # avd is the independent variable

E[i]=sum(as.character(Test[m,Dv])!=predict(a,Test[m,])$avd)/n1}

mean(E)

#3bagging

library(adabag)

E=rep(0,Z)

for(i in 1:Z){m=mm[[i]]; # mm[[i]] is the ith subscript set extracted by Fold() function: i=1, 2, …, Z

n1=length(m)

a=bagging(avd~.,Test[-m,]) # avd is the independent variable

E[i]=sum(as.character(Test[m,Dv])!=predict(a,Test[m,])$avd)/n1}

mean(E)

#4Nearest neighbor algorithm

library(kknn)

E=rep(0,Z)

for(i in 1:Z){m=mm[[i]]; # mm[[i]] is the i th subscript set extracted by Fold() function: i=1, 2, …, Z

n1=length(m)

a=kknn(avd~.,k=6,train=Test[-m,],test=Test[m,]) # avd is the independent variable

E[i]=sum(Test[m,Dv])!=a$avd/n1}

mean(E)

#5Support vector machine

library(kernlab)

E=rep(0,Z)

for(i in 1:Z){m=mm[[i]]; # mm[[i]] is the i th subscript set extracted by Fold() function: i=1, 2, …, Z

n1=length(m)

a=ksvm(avd~.,data=Test[-m,]) # avd is the independent variable

E[i]=sum(Test[m,Dv]!=predict(a,Test[m,]))/n1}

mean(E)

#6Random forest

library(randomForest)

E=rep(0,Z)

for(i in 1:Z){m=mm[[i]]; # mm[[i]] is the i th subscript set extracted by Fold() function: i=1, 2, …, Z

n1=length(m)

a=randomForest(avd~.,data=Test[-m,]) # avd is the independent variable

E[i]=sum(Test[m,Dv]!=predict(a,Test[m,]))/n1}

mean(E)

#7Artificial neural network

library(nnet)

E=rep(0,Z)

for(i in 1:Z){m=mm[[i]]; # mm[[i]] is the i th subscript set extracted by Fold() function: i=1, 2, …, Z

n1=length(m)

mc=setdiff(1:n,m)

a=nnet(avd~.,data=Test,subset=mv,size=7,range=0.1,decay=0.01,maxit=200) # avd is the

independent variable

E(i)=sum(Test[m,Dv]!=predict(a,Test[m,])$avd)/n1}

mean(E)

这里再附上其他两种机器学习方法（参考文献：

Kim, H., Lee, K.K., Jeon, J., Harris, W. A., & Lee, Y. H. (2020). Domestication of Oryzaspecies eco-evolutionarily shapes bacterial and fungal communities in rice seed. Microbiome, 8(1), 20.doi:10.1186/s40168-020-00805-0）

#加载包：

library(phyloseq)

library(dplyr)

library(ggplot2)

# Data preparation

table <- otu_table(phy.clean.log)

table <- t(table)

k <- 129-75

result <- c(rep(0,75),rep(1,k))

table <- cbind(table, result)

dataset <- as.data.frame(table)

# Encoding the target feature as factor

dataset$result = factor(dataset$result, levels = c(0, 1))

dataset$result

# Splitting the dataset into the Training set and Test set

library(caTools)

set.seed(123)

split = sample.split(dataset$result, SplitRatio = 0.66)

training_set = subset(dataset, split == TRUE)

test_set = subset(dataset, split == FALSE)

#8朴素贝叶斯Naive bayes

library(e1071)

packageVersion('e1071') ## 鈥?1.7.0鈥?

classifier_nb = naiveBayes(x = training_set[-ncol(dataset)],

y = training_set$result)

nb_pred = predict(classifier_nb, newdata = test_set[-ncol(dataset)])

cm = table(test_set[, ncol(dataset)], nb_pred)

print(cm)

#9逻辑回归logistic regression

classifier_lr <- glm(formula = result ~.,

family=binomial,

data = training_set)

lr_pred <- predict(classifier_lr, type = 'response', newdata=test_set[-ncol(dataset)])

lr_pred <- ifelse(lr_pred > 0.5, 1, 0)

cm = table(test_set[, ncol(dataset)], lr_pred)

print(cm)

构建功能基因共发生网络代码如下

data <- read.table("gene_abundance.txt",header = TRUE,row.names = 1,sep = "\t")

# read in a gene abundance table

data <- t(data)

library(psych)

occor =corr.test(data,method = "spearman",adjust = " fdr") # calculate spearman’s

correlation

occor.r = occor$r # extract r-value of spearman’s correlation

occor.p = occor$p # extract p-value of spearman’s correlation

occor.r[occor.p>0.01|abs(occor.r)<0.75] = 0 # r-value and p-value filtering

write.table(occor.r,"pvalue.csv",sep="\t",quote=F,col.names=NA) # output results，构建的网络可以直接在Gephi或者在R内实现可视化详见（一文学会网络分析——Co-occurrence网络图在R中的实现，利用Gephi软件绘制网络图，Gephi可视化网络的一个简单示例操作视频）

群落构建过程指标βNTi代码在补充材料不放上面了。

主要结果

1. 稀释降低了土壤细菌群落的α多样性和稳定性

图1：重新组装的细菌群落的多样性和变异。A重组装群落系统发育多样性。B稀释与βNTI之间的关系。C物种丰富度与AVD之间的关系。D基于随机森林评估功能分类的重要性。

2. 稀释简化了功能基因共发生网络

表1：功能基因共发生网络拓扑性质特性与随机网络

图2：稀释梯度下的功能基因共发生网络

由表1与图2可知随着稀释微生物功能基因网络拓扑结构参数降低，网络更加松散，一些特殊功能会随着稀释消失。在稀释程度低的群落中特殊代谢功能占主导，而在稀释程度高的群落中广泛代谢功能占主导。

3. 特定功能的关键类群有利于微生物物群稳定性。

图3：两个关键功能的分类注释。A8个属于氮代谢与44个磷酸盐、次磷酸盐代谢在pH与稀释梯度上的相对丰度。B相应属的相对丰度。

由于“氮代谢”和“磷酸盐和次磷酸盐代谢”的特殊代谢功能是随机森林分类方法和功能基因共发生网络分析中与微生物群稳定性相关的关键节点的识别中最重要的功能，因此进行了分类注释以识别相关的分类群。

结论

结果表明，稀释显著降低了系统发育多样性，简化了功能基因共发生网络的模块化，降低了土壤微生物群落的稳定性。“氮代谢”和“磷酸盐和次磷酸盐代谢”的特殊代谢功能与土壤微生物群落的稳定性有关。这些关键功能和分类群Nitrospira 与Gemmatimonas可能为预测微生物群落的变化提供进一步的见解，并增加操纵微生物群功能的可能性。并利用关键分类群对微生物群落稳定性的可能贡献来改善生态系统服务。