題目：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可能為預測微生物群落的變化提供進一步的見解，並增加操縱微生物群功能的可能性。並利用關鍵分類群對微生物群落穩定性的可能貢獻來改善生態系統服務。