原核微生物群落随机性和确定性装配过程的计算方法   

软件

1. R version 3.4.3

实验步骤

1. 参照《Extraction and 16S rRNA Sequence Analysis of Microbiomes Associated with Rice Roots》(DOI:10.21769/BioProtoc.2884) (Edwards et al.，2018)中的方法对微生物群落进行16S rRNA序列分析获得样品的群落组成和系统发育关系树，对应的两个文件分别为OTU表格out_table.txt，和树文件rep_set.tre。
2. 系统发育群落组成的计算参数
          针对单个群落，使用净关系指数 (NRI) 和最近分类单元指数 (NTI) 来反映物种集聚的谱系结构。二者分别对应平均谱系距离 (MPD) 和最近谱系距离的 (MNTD) 的标准化测度。
          针对多个样本间的群落比较，则是基于MPD和MNTD使用βMPD和βMNTD计算群落间谱系beta多样性，而对于群落间聚集或分散受确定性或随机性因素影响的判断则是基于NRI和NTI使用βNRI和βNTI计算群落间系统发育beta多样性的显著性。
          以上参数均可在R软件包“picante”中完成计算 (Webb et al.，2002)。本文将分别基于上述参数介绍具体的含义与计算方法。
   
   2.1

   MNTD/MPD
          MNTD是mean-nearest-taxon-distance (最近种间平均进化距离) 的缩写，代表任意一个物种与群落中与自己亲缘关系最近的一个物种的平均进化距离；它们的值越小，物种间的亲缘关系越近。
   计算公式如下：
   
   fik：群落k中OTUi的相对丰度；
   nk：群落k中OTU的数量；
   ∆_ikjk：OTUi与同样在群落k中的所有其他OTUj之间的最小系统发育距离。
          MPD是mean phylogenetic distance (种间平均进化距离)的缩写，是群落中所有物种对间平均的进化距离。MPD分别计算样本内所有种两两之间 (两两之间成对) 的系统发育距离，再除以两两之间的对数。其计算可在R的picante包中使用mpd()函数完成。
   2.2

   NTI/NRI
          NRI和NTI这两个系数可以描述实际群落与零群落之间的平均谱系距离，它们都利用零群落中谱系距离的标准差进行了标准化。
          NTI是nearest-taxon-index (最近分类单元指数) 的缩写，是对样本中每个分类单元到最近分类单元的系统发育距离的标准化度量，量化终端聚类的程度，对应最近谱系距离的(MNTD)的标准化测度。为了评估系统发育群落结构的非随机程度，将OTU及其相对丰度在系统发育的各个尖端进行随机分配 (null.model=‘taxa.labels’in‘ses.mntd’)。NTI量化观察到的MNTD与零分布均值的标准偏差数(999个随机)。NTI是R picante包中“ses.mntd”输出的相反数。计算公式如下：
   
          NRI是net relatedness index (净关系指数)的缩写，是样本中分类单元相对于系统发育的平均成对系统发生距离的标准化度量，基于适当的物种库，量化一棵树上的分类单元的总体聚类程度，对应平均谱系距离(MPD)的标准化测度。NRI是R picante包中“ses.mpd”输出的相反数。计算公式如下：
   
          MPDobs (或MNTDobs) 为各样地观测值，MPDnull (或MNTDnull) 和sdMPDnull (或sdMNTDnull) 则分别对应由零模型得到的平均值和标准偏差。
          对于单个群落，NTI (或NRI) 大于+2表示共存的分类单元比偶然预期的关系更为紧密 (系统发生聚类)；NTI (或NRI) 小于-2表示共存的分类单元比偶然预期的亲缘关系更远 (系统发生过度分散)。在所有的群落中获得的平均NTI (或NRI) 与零值有显著差异，表示存在平均聚类 (>0) 或过度分散 (<0) (Kembel, 2009)。
   2.3

   βMNTD/βMPD
          系统发育随时间或空间的变化 (系统发育β多样性) 使用beta最近谱系距离βMNTD (或beta平均谱系距离βMPD) 和beta最近谱系指数βNTI (或beta平均谱系指数βNRI) 来量化，它们分别是MNTD (MPD) 和NTI (NRI) 的组间 (between-assemblage) 相似物 (Webb et al.，2009；Kembel et al.，2011)。
   βMNTD是发生在两个不同群落中的最近亲属间的丰度加权平均系统发育距离。这种度量量化了群落装配之间的系统发育距离；βMNTD被使用是因为它强调短的系统发育距离的系统发育周转率，其中大力支持系统发育信号的假设。其计算公式如下：
   
          min(∆_ikjm)：群落k中的OTUi与群落m中所有OTUj之间的最小系统发育距离(其他变量请参见上文)。
          βMPD是发生在两个不同群落中的任意两个物种间的丰度加劝平均系统发育距离。其计算公式如下：
   
          ∆_ikjm：群落k中的OTUi与群落m中所有OTUj之间的平均系统发育距离(其他变量请参见上文)。
   
   2.4

   βNTI/βNRI
          为了推断出对特定群落装配之间观察到的周转率起主要作用的生态因素，鉴于随机性生态过程的主导地位，使用零模型方法生成一个预期水平的βMNTD (βMPD)。为了量化观察到的βMNTD (βMPD) 值和零βMNTD (βMPD) 值分布之间偏差的大小和方向，使用了β-最近谱系指数βNTI (β-平均谱系指数βNRI)。通过将整个系统发育过程中的OTU随机化并重新计算βMNTD (βMPD) 999次，可以得到零βMNTD (βMPD) 分布。βNTI (βNRI) 是观察到的βMNTD (βMPD) 与零βMNTD (βMPD) 分布均值的标准偏差数。
   
          对于成对的群落比较，βNTI (βNRI) <-2或>+2分别表示小于或大于预期的系统发育周转。在所有成对比较中获得的平均βNTI (βNRI) 与零有显著差异，表示大于 (>0) 或小于 (<0) 预期平均周转 (Webb et al.，2009；Kembel et al.，2011)。

结果与分析

1. R中计算示例——针对一对样本的计算
          上述参数的计算可以在R中实现，现以两个来自海洋环境的原核微生物群落的数据为例子，计算上述参数。样品分别命名为S1和S2, OTU文件命名：out_table.txt；树文件命名：rep_set.tre；相应的脚本：betaNTI for one pair sample.R，可在该地址下载：https://doi.org/10.5281/zenodo.4506365。
   
   library(picante)
   library(ape)
   comm<-
   read.table("out_table.txt",sep='\t',row.names=1,header=T)
   #读入OTU文件
   #读入树文件
   phy<-read.tree("rep_set.tre")
   #使树文件和OTU表文件对齐
   prune_tree<-prune.sample(comm,phy)
   #计算每个OTU之间距离矩阵
   phydist<-cophenetic(prune_tree)
   #计算MNTD
   mntd<-
   ses.mntd(comm,phydist,null.model="taxa.labels",abundance.weighted=T, runs=999)
   
   该步骤计算结果如下：
   

   	ntaxa	mntd.obs	mntd.rand.mean	mntd.rand.sd	mntd.obs.rank	mntd.obs.z	mntd.obs.p	runs
   S1	498	0.13913	0.1712452	0.02944938	108	-1.090523	0.108	999
   S2	586	0.1432301	0.1591859	0.01302431	90	-1.225076	0.09	999

   其中，ntaxa Number of taxa in community，群落中OTU的数目；mntd.obs Observed MNTD in community，观测到的群落MNTD值；mntd.rand.mean Mean MNTD in null communities，随机群落的MNTD均值；mntd.rand.sd Standard deviation of MNTD in null communities，随机群落MNTD值的标准差；mntd.obs.rank Rank of observed MNTD vs. null communities，观察到的群落的MNTD值与随机群落MNTD值的排序；mntd.obs.z Standardized effect size of MNTD vs. null communities (=(mntd.obs-mntd.rand.mean)/mntd.rand.sd, equivalent to -NTI) 观察到的群落的MNTD与随机群落的MNTD值的标准化数值大小，相当于-NTI值；mntd.obs.p P-value (quantile) of observed MNTD vs. null communities (= mntd.obs.rank/runs + 1) 观察到的群落的MNTD与随机群落的MNTD值比较的p值；runs Number of randomizations，随机化的次数。
   
   #计算βMNTD值
   comdist.result<-comdistnt(comm,phydist)
   
   该步骤计算结果如下：
   

   	S1
   S2	0.06070395

   #以下为βMNTDnull的计算过程
   f<-function(){
   g<-randomizeMatrix(comm,null.model=c("frequency","richness","independentswap","trialswap"),iterations=1000)#针对两组样方实现随机群落构建
      fc<-comdistnt(g,phydist)#针对随机构建的两个样方进行群落间βMNTD的计算
   }
   mt<-apply(replicate(999,f()),1,mean)#上述f计算过程重复进行999次，得到βMNTD的平均值
   mt
   mt.sd<-apply(replicate(999,f(comm)),1,sd)#计算上述999次计算得到的βMNTD值的标准差
   mt.sd
   βNTI<-(comdist.result-mt)/mt.sd
          最终，得到的βNTI值为4.93，说明S1和S2两个样方的群落之间的变化主要由确定性过程影响。
2. R中计算示例——针对多样本的计算
          以上是针对1对样本之间计算βNTI值的介绍，给出了每一步操作的目的，下面介绍针对多样本间βNTI值计算的方法。该方法参考Stegen et al. 2013，ISME (Stegen et al.，2013)，依然在R中进行操作，相应的脚本：betaNTI for multiple pairs of samples.R 可在前文地址下载，具体如下：
   
   ##切换到包含OTU表和系统发育树的文件目录
   #加载这个数据包
   library(picante)
   #读取OTU表，依然是样品名在行，OTU种类在列
   otu = read.csv("out_table.txt", header=T, row.names=1)
   #读取系统发育树
   phylo = read.tree("rep_set.tre")
   #确保系统发育树上的名称与otu表中的名称排序一致
   match.phylo.otu = match.phylo.data(phylo, t(otu))
   #计算betaMNTD，这里abundance.weighted=T表示考虑物种丰度，如果不设定该参数，则默认abundance.weighted=F，不考虑物种丰度
   beta.mntd.weighted = as.matrix(comdistnt(t(match.phylo.otu$data),cophenetic(match.phylo.otu$phy),abundance.weighted=T))
   #检查beta.mntd.weighted矩阵中的列名是否和系统发育树以及OTU表中一致，结果应该为T
   identical(colnames(match.phylo.otu$data),colnames(beta.mntd.weighted))
   #同上，只是一个检查，结果应该为T
   identical(colnames(match.phylo.otu$data),rownames(beta.mntd.weighted))
   ##计算betaMNTD随机值
   #随机化次数
   beta.reps = 999
   rand.weighted.bMNTD.comp=array(c(-999),dim=c(ncol(match.phylo.otu$data),ncol(match.phylo.otu$data),beta.reps));
   for (rep in 1: beta.reps) {
   #这一步执行的随机化是基于对系统发育树中OTU标签的随机重排进行的，如果沿用1中的方法对每一个群落进行随机构建后再计算betaMNTD,则可以选择使用rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(randomizeMatrix(comun, null.model = c("frequency", "richness","independentswap", "trialswap"), iterations = 1000),cophenetic(match.phylo.comun$phy),abundance.weighted=T,exclude.conspecifics = F))命令完成该步运算
   rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(match.phylo.otu$data),taxaShuffle(cophenetic(match.phylo.otu$phy)),abundance.weighted=T,exclude.conspecifics = F));
     print(c(date(),rep));
   }
   
   ##计算βNTI
   weighted.bNTI = matrix(c(NA),nrow=ncol(match.phylo.otu$data),ncol=ncol(match.phylo.otu$data));
   for (columns in 1:(ncol(match.phylo.otu$data)-1)) {
     for (rows in (columns+1): ncol(match.phylo.otu$data)) {
   #按照βNTI计算公式计算每一对样本间的βNTI值
     rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
     weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
   rm("rand.vals");
     };
   };
   rownames(weighted.bNTI) = colnames(match.phylo.otu$data);
   colnames(weighted.bNTI) = colnames(match.phylo.otu$data);
   weighted.bNTI;
   #得到全部样本对间的βNTI值并写入文件
   write.csv(weighted.bNTI,"weighted_bNTI.csv",quote=F);
   
   
3. 多组样品计算结果示例
          针对多组样品多个样方的βNTI计算按照上述示例2分别计算两两样方之间的βNTI值，并按照设计的样方分组分析βNTI在不同分组中的变化范围。本例展示来自渤海湾6个不同站位表层海水样品在不同月份 (5、8、10月) 的βNTI统计值。18个样品的物种组成数据已上传至前文地址。按照示例2分别计算了unweighted βNTI和weighted βNTI，并将计算结果在Origin中进行box-plot制图，结果如下：
   
   
   图1. 6个站位在不同月份的β-NTI值分布
   
          如图1所示，我们分别计算了unweighted βNTI和weighted βNTI，unweighted βNTI是不考虑物种丰度的计算结果，weighted βNTI 是考虑物种丰度的计算结果。根据unweighted βNTI 的统计结果，显示5月份、8月份、10月份的海水样品βNTI值主要集中在-1~6之间，说明在各个月份内，不同的站位间的群落变化既有受随机性因素影响的情况，也有受决定性因素影响的情况；而两两月份间的样品分析显示βNTI值均大于2，说明月份间的微生物群落结构变化主要受决定性因素的影响。weighted βNTI 的统计结果与unweighted βNTI 的统计结果有所不同，显示月份内和月份间的样品对之间的βNTI值均有较大的波动，在-2~9之间，说明考虑物种丰度变化时，该区域的微生物群落结构受取样点和时间的影响均较大。
   

结果与分析

1. 软件Phylocom也对上述参数提供了标准化的解决方案，本文主要提供了基于R中picante数据包进行计算的方法 (Webb et al.，2002)。以上我们给出了MNTD、βMNTD、NTI、βNTI的计算方法，对应的MPD、βMPD、NRI、βNRI的计算则只需更改相应的函数即可获得。这两组参数在生物学意义上既相似又有不同，它们分别从不同角度描述了群落在系统发育树上的聚集模式，MNTD侧重指示分支末端的谱系关系，MPD代表了物种间亲缘关系的平均情况。基于两种谱系测度 (MPD和MNTD) 估计群落内部与群落之间成员亲缘关系的强度时，它们对解释变量的灵敏度可能不尽相同。在实际研究中，针对微生物群落的研究显示MNTD对肠道微生物群落的宿主选择性更敏感 (Youngblut et al.，2019)，而针对草本植物群落的研究则显示NRI和βMNTD与海拔梯度之间存在显著线性关系，而NTI和βMPD并未发现显著趋势 (赵鸣飞等，2017)。这很可能跟指标本身的算法特性有关。具体来说，对区域种库整体上具有谱系聚集趋势的情况，该两指标同时倾向于取得较高值；而谱系聚集发生于谱系树的多个浅层分支上时，仅NTI易于得到较高值 (Vamosi et al.，2009)。
2. 关于样品采样深度对该方法计算结果的影响，目前尚无系统性的报道，根据已有报道基于该方法进行的研究显示Stegen等 (Stegen et al.，2012) 研究随机和确定性过程对不同深度的沉积物样品的相对影响时，对样品的序列取样深度为500条；Wang等 (Wang et al.，2013) 对不同生态系统细菌群落系统发育多样性受随机性和决定性因素影响的研究中，每个样品的序列取样深度为1,000条；Yu等 (Yu et al.，2018) 研究空间尺度对种植小麦土壤中微生物群落随机性和决定性相对作用的影响时，对每个样品的序列取样深度为20,005条。可见，随着测序技术越来越成熟，用于进行群落系统发育分析的取样深度也越来越大。另外，关于物种数对系统发育多样性的影响，以MNTD和NTI为例，其他系统发育多样性指数类似：MNTDobs是观测到的MNTD，MNTDnull是随机打乱进化树之后得到的MNTD值，这个过程重复多次 (~999次)，可得到MNTDnull的分布。可近似转化为均值为0，标准差为1的标准正态分布。而这个分布其实是所有物种施加影响结果的集合。NTI和βNTI计算过程已经排除了物种数量对系统发育的影响。
3. 关于使用零模型 (null model) 方法解释原核微生物群落构建：零模型用于验证基于自然实验的生态学假说，它以计算机技术为基础，从已经获得的实验数据出发，构建受控的生态随机模型，最后用标准的统计检验方法验证理论假说。在解释原核微生物群落构建机制中，零模型是群落数据或者进化树数据进行随机化的方式。运行零模型主要是为了解，在随机状况下，物种系统发育关系的统计分布，从而推断真实的群落内物种的系统发育关系是否随机，以便了解环境过滤或者随机过程的相对重要性。所以在探讨这些生态学过程的相对重要性，计算NTI和βNTI的时候，必须要考虑零模型。目前，该方法使用基于randomizeMatrix函数提供的零模型构建主要包括物种内群落数据矩阵丰度 (维持物种出现频率) 的随机化 (frequency)、样本内群落数据矩阵丰度 (维持样本物种丰富度) 随机化 (richness)、利用独立交换算法 (Gotelli，2000) 对群落数据矩阵进行随机化 (保持物种发生频率和样本物种丰富度) (independentswap)、使用试交换算法 (Miklos and Podani，2004) 随机化群落数据矩阵 (保持物种发生频率和样本物种丰富度) (trialswap)。对于生态零模型，Gotelli指出“它是基于生态数据的随机化或者基于已知分布、设想分布的随机样本基础上的格局生成模型。与已知数据相比，一些数据的性质在零模型中保持不变，另外一些则允许随机变化，这样就产生了一个新的格局总体。”从中可以看出零模型的构建是从实际所得数据出发的，并非真正的实验。对于原核微生物群落而言，环境中的微生物种类只能尽最大的可能通过高通量、高分辨率的分析方法获得，但是仍然无法完全获知其中的全部微生物组成，因此，零模型的构建也只能依据已知样品已获得的物种信息进行预测，对于未能获得的物种信息则尚不能纳入到零模型的构建中。但是，针对样方之间的比较，我们更关注设计样方间的相对变化规律，因此，取样深度的一致性是分析对比的前提。
4. 注意，文中关于NTI/NRI、βNTI/βNRI的结果解释时其显著性的判断依据是其临界值是否大于|2|，实际上这是依据标准正态分布的95%置信区间得到的，一般认为NRI或者NTI大于1.96或者小于-1.96的结果，在95%的水平是显著的。这是因为计算NTI/NRI或者βNTI/βNRI都是基于Standard effective Size的变型，其本质都是基于多次随机化的群落或者进化树计算的，在随机样本量很大的情况下，MNTD/MPD或者βMNTD/βMPD的分布会接近正态分布(张金龙，2019)。所以文献中提及的这几个临界值接近于2，实际是接近1.96。
   

致谢

感谢国家重点研发计划重点专项项目 (2018YFA0901200)；天津市科技计划项目(18ZXSZSF00100,19YFZCSF00750)；国家自然科学基金项目 (41977200) 的支持。感谢Campbell O Webb，Steven W Kembel，James C Stegen等人先前的研究工作。感谢西北农林科技大学焦硕教授的热心指导！

参考文献

1. Chase, J. M. and Myers, J. A. (2011). Disentangling the importance of ecological niches from stochastic processes across scales. Philos Trans R Soc Lond B Biol Sci 366(1576): 2351-2363.
2. Dini-Andreote, F., Stegen, J. C., van Elsas, J. D. and Salles, J. F. (2015). Disentangling mechanisms that mediate the balance between stochastic and deterministic processes in microbial succession. Proc Natl Acad Sci U S A 112(11): E1326-1332.
3. Edwards, J., Santos-Medellín, C. and Sundaresan, V. (2018). Extraction and 16S rRNA Sequence Analysis of Microbiomes Associated with Rice Roots. Bio-protocol 8(12): e2884.
4. Gotelli, N. J. (2000) Null Model Analysis of Species Co-Occurrence Patterns. Ecology 81: 2606-2621.
5. Kembel, S. W. (2009). Disentangling niche and neutral influences on community assembly: assessing the performance of community phylogenetic structure tests. Ecol Lett 12(9): 949-960.
6. Kembel, S. W., Eisen, J. A., Pollard, K. S. and Green, J. L. (2011). The phylogenetic diversity of metagenomes. PLoS One 6(8): e23214.
7. Miklos, I. and Podani, J. (2004). Randomization of presence-absence matrices: Comments and new algorithms. Ecology 86-92.
8. Stegen, J. C., Lin, X., Konopka, A. E. and Fredrickson, J. K. (2012). Stochastic and deterministic assembly processes in subsurface microbial communities. ISME J 6(9): 1653-1664.
9. Stegen, J. C., Lin, X., Fredrickson, J. K., Chen, X., Kennedy, D. W., Murray, C. J., Rockhold, M. L. and Konopka, A. (2013). Quantifying community assembly processes and identifying features that impose them. ISME J 7(11): 2069-2079.
10. Vamosi, S. M., Heard, S. B., Vamosi, J. C. and Webb, C. O. (2009). Emerging patterns in the comparative analysis of phylogenetic community structure. Mol Ecol 18(4): 572-592. 
11. Vellend, M. (2010). Conceptual synthesis in community ecology. Q Rev Biol 85(2): 183-206.
12. Wang, J., Shen, J., Wu, Y., Tu, C., Soininen, J., Stegen, J. C., He, J., Liu, X., Zhang, L. and Zhang, E. (2013). Phylogenetic beta diversity in bacterial assemblages across ecosystems: deterministic versus stochastic processes. Isme J 7: 1310-1321.
13. Webb, C. O., Ackerly, D. D. and Kembel, S. W. (2008). Phylocom: software for the analysis of phylogenetic community structure and trait evolution. Bioinformatics 24(18): 2098-2100.
14. Webb, C. O., Ackerly, D. D., McPeek, M. A. and Donoghue, M. J. (2002). Phylogenies and community ecology. Annu Rev Ecol Syst 33: 475-505.
15. Youngblut, N. D., Reischer, G. H., Walters, W., Stalder, G., Ley, R. E. and Farnleitner, A. H. (2019). Host diet and evolutionary history explain different aspects of gut microbiome diversity among vertebrate clades. Nat Commun 10(1):2200.
16. Shi, Y., Li, Y., Xiang, X., Sun, R., Yang, T., He, D., Zhang, K., Ni, Y., Zhu, Y. G., Adams, J. M. and Chu, H. (2018). Spatial scale affects the relative role of stochasticity versus determinism in soil bacterial communities in wheat fields across the North China Plain. Microbiome 6(1): 27.
17. Zhou, J., Deng, Y., Zhang, P., Xue, K., Liang, Y., Van Nostrand, J. D., Yang, Y., He, Z., Wu, L., Stahl, D. A., Hazen, T. C., Tiedje, J. M. and Arkin, A. P. (2014). Stochasticity, succession, and environmental perturbations in a fluidic ecosystem. Proc Natl Acad Sci U S A 111(9): E836-845.
18. 赵鸣飞, 薛峰, 王宇航,等.(2017). 山西芦芽山针叶林草本层群落谱系结构与多样性的海拔格局. 植物生态学报, 41(7): 707-715. 
19. 张金龙. (2019). 群落系统发育分析常见问题合集. 科学网, 2019-1-6.