谱系分析在群落生态学中的应用

研究背景

群落是如何构建的是群落生态学的研究核心。早期的研究通常聚焦于物种的分类单元，围绕以物种间竞争为核心理论的模型展开，例如戴蒙德法则 (Diamond’s rule) (Diamond, 1975) 和Fox优势状态模型 (Fox’s favored states model) (Fox, 1987) 等。近20年，随着分子生物学技术和系统发育研究方法的发展与应用，生态学家渐渐将系统发育关系与群落构建结合在一起，从进化的角度探究群落构建的过程 (Webb, 2000; Webb et al., 2002; Cavender-Bares et al., 2009)，其中Webb (2000) 根据物种的亲缘关系提出了净种间亲缘关系指数 (NRI) 和净最近种间亲缘关系指数 (NTI) 来检测群落构建的机制。
    Webb et al. (2002) 利用群落系统发育结构建立了一个启发式框架，用于揭示群落构建的主要驱动力。这个框架假设性状与环境的匹配 (即环境过滤，habitat filtering) 和竞争排斥 (competitive exclusion) 是构建群落结构的主要驱动力。当亲缘关系密切的物种聚集程度超过期望值 (系统发育聚集，phylogenetic clustering) 时，其根本原因可能是环境对具有相似适应性性状的物种的筛选 (性状保守性，trait conservatism)。反之，当群落中物种亲缘关系比预期小时 (系统发育过度离散，phylogenetic overdispersion)，可能是由竞争引起的保守性状 (conserved traits) 过度发散 (overdispersion) 导致的，也可能是环境对生态学上重要趋同性状 (convergent traits) 筛选的结果。但是，后续有研究发现竞争也会引起系统发育聚集 (Mayfield and Levine, 2010)，并取决于物种的生态位差异 (niche differences)、竞争能力差异 (competitive ability differences) 以及谱系信号 (phylogenetic signal) 的强弱 (Cavender-Bares et al., 2009; Mayfield and Levine, 2010; Davies, 2021)。此后，Cadotte et al. (2019) 总结了使用谱系和功能信息进行群落生态学分析的常见步骤：(a) 谱系信号检验；(b) 谱系和功能多样性相关性分析；(c) 零模型检验谱系和功能发散格局；(d) 谱系和功能多样性及发散格局与其它因子 (通常为环境因子) 的相关性 (Cadotte et al., 2019)。

一、软件和核心软件包

1. 文本中所有软件的运行及数据处理和分析均在macOS操作系统上进行 (在Windows及其它系统中部分代码格式可能需要自行调整)
2. R software v.4.0.3 (https://cran.r-project.org/src/base/R-4) 
3. RStudio v.1.4.1106 (https://www.rstudio.com/products/rstudio/)
4. ape v.5.4.1 (Paradis et al., 2004；https://cran.r-project.org/web/packages/ape/index.html)
5. evobiR v.1.1 (https://cran.r-project.org/web/packages/evobiR/index.html)
6. phytools v.0.7.70 (Revell, 2012；https://cran.r-project.org/web/packages/phytools/index.html)
7. FD v.1.0.12 (Laliberté et al., 2014；https://cran.r-project.org/web/packages/FD/ index.html)
8. picante v.1.8.2 (Kembel et al., 2010；https://cran.r-project.org/web/packages/pic ante/index.html)
9. reshape2 v.1.4.4 (Wickham, 2012; http://cran.ms.unimelb.edu.au/web/packages/ reshape2/)

二、研究方法

1. 数据准备
   本方案将以浙江千岛湖33个岛屿上的蚂蚁群落数据为例 (Zhao et al., 2020) (本文数据和代码请见附件)。
   1.1

   谱系树
   目前，一些类群已经有相对完整的谱系树 (即supertree)，研究者可以通过这些supertree获得仅包含自己研究物种的谱系树。例如植物 (V.phylomaker或phylomatic：https://phylodiversity.net/phylomatic/) (Jin and Qian, 2019)、鸟类 (birdtree：http://birdtree.org/) 和其他脊椎动物 (vertlife: https://vertlife.org/phyl osubsets/)。学者也可以直接通过分子测序数据构建谱系树。在得到谱系树后，可以通过ape包中的read.nexus或read.tree函数 (Paradis et al., 2004) 将谱系树读取进入R。前者读取Nexus格式的谱系树，后者读取Newick格式的谱系树。
   
   #读取谱系树
   library(ape) #加载ape包
   ant_tree <- read.nexus("islandtree.tre") #谱系树也可以是txt文件格式。
   
   1.2

   群落数据 community data
   矩阵 (matrix) 或数据框 (data.frame) 是群落生态学分析中常用的数据格式，其中行为地点或样点 (sites/samples)，列为物种 (taxa)。群落数据形式可以是0/1
   (presence/absence) 也可以是每个物种在不同样点里的多度 (abundance)。本文将以0/1数据为例：
   
   #读取群落数据
   #0/1群落数据
   ant_com <- read.delim("ant_com.txt", h = T, row.names = 1, sep = ",")
   
   1.3

   性状数据 trait data
   ant_trait <- read.delim("ant_trait.txt", h = T, row.names = 1, sep = ",") head(ant_trait)
   
   
   
   其中AL是体长 (Weber’s length)、ED_AL是相对眼宽 (relative eye width)、LL_AL是相对后足长 (relative leg length)、SL_AL是相对触角柄节长 (relative scape length)、ML_HW是相对上颚长 (relative mandible length) (Zhao et al., 2020)。
   
2. 数据分析
   2.1

   谱系信号 (phylogenetic signal)
   谱系信号表示近缘物种彼此相似的趋势，通常认为谱系关系近的物种，其功能性状越相似，即系统发育保守性。群落系统发育分析的前提假设就是物种的性状是系统发育保守的。Blombergs’ K (Blomberg et al., 2003) 和 Pagels’ λ (Pagel, 1999) 是两个定量区分性状是否具有系统发育保守性的指标。两者均是通过比较性状的观察值与布朗运动 (Brownian motion) 模型预测值的差异，以判断性状系统发育保守性的程度。两者皆是数值越大，性状系统发育保守性越强，其中λ的范围一般是0至1，K的值可以超过1。虽然λ和K的值可以反应性状系统发育保守性的程度，还需通过零模型检验其显著性。对于指标K，具体做法是将谱系树上物种对应的性状数据随机打乱多次 (通常为999次) 模拟没有谱系信号的情形 (Kraft and Ackerly, 2010)。而对于指标λ，则通过似然比检验 (likelihood ratio test) 比较观察值是否显著不同于λ = 0 (即没有谱系信号) 的零假设。上述计算可以通过phytools包中的phylosig函数实现 (Revell, 2012)。
   
   library(evobiR)
   library(phytools)
   #首先将性状数据按照谱系树上的名字顺序重新排列，详见?ReorderData()
   trait_sig <- ReorderData(ant_tree, ant_trait, taxa.names = "row names")
   
   #计算5个性状的Blomberg’s K值
   signal_AL_K <- phylosig(ant_tree, trait_sig$AL, method = "K", test = T, nsim = 999)
   signal_ED_AL_K <- phylosig(ant_tree, trait_sig$ED_AL, method = "K", test = T, nsim = 999)
   signal_LL_AL_K <- phylosig(ant_tree, trait_sig$LL_AL, method = "K", test = T, nsim = 999)
   signal_SL_AL_K <- phylosig(ant_tree, trait_sig$SL_AL, method = "K", test = T, nsim = 999)
   signal_ML_HW_K <- phylosig(ant_tree, trait_sig$ML_HW, method = "K", test = T, nsim = 999)
   
   #计算5个性状的Pagel’s λ值
   signal_AL_λ <- phylosig(ant_tree, trait_sig$AL, method = "lambda", test = T)
   signal_ED_AL_λ <- phylosig(ant_tree, trait_sig$ED_AL, method = "lambda", test = T)
   signal_LL_AL_λ <- phylosig(ant_tree, trait_sig$LL_AL, method = "lambda", test = T)
   signal_SL_AL_λ <- phylosig(ant_tree, trait_sig$SL_AL, method = "lambda", test = T)
   signal_ML_HW_λ <- phylosig(ant_tree, trait_sig$ML_HW, method = "lambda", test = T)
   
   我们以体长的谱系信号结果为例，虽然K值约为0.2，是小于1的，但是与没有谱系信号的随机结果存在显著差异 (P值小于0.05)。因此，我们认为体长这一性状是存在谱系信号的，但是其系统发育保守性程度并不高 (K = 0.2)。同样地，我们看到λ值约为0.87，接近1，说明通过λ指标计算得到的系统发育保守性程度较高，其P值小于0.05，说明体长具有显著的谱系信号。为什么λ值会高于K值，这是因为λ对于谱系信号的检测更加敏感 (详见Münkemüller et al., 2012)。
   
   signal_AL_K
   #Phylogenetic signal K: 0.202679
   #P-value (based on 999 randomizations): 0.004004
   
   signal_AL_λ
   #Phylogenetic signal lambda: 0.873473
   #logL(lambda): -87.8682
   #LR(lambda=0): 44.8614
   #P-value (based on LR test): 2.11482e-11
   
   2.2

   谱系多样性 (phylogenetic diversity) 
   最早的谱系多样性指标之一是Faith于1992年提出的phylogenetic diversity (PD) (Faith, 1992)，也称为Faith’s PD。其定义为一个局域群落内所有物种在谱系树上所占枝长的总和。可以通过picante包中的pd函数计算得到 (Kembel et al., 2010)。pd函数同时也可以计算出每个群落的物种数 (species richness, SR)。通过结果，我们可以发现，物种数高的群落，其Faith’s PD值也高。
   
   #Faith’s PD谱系多样性
   library(picante)
   PD <- pd(ant_com, ant_tree, include.root = TRUE)
   head(PD)
   
   # PD SR
   #S01 2.948611 66
   #S02 2.668410 59
   #S03 3.013944 68
   #S04 2.280161 54
   #S05 2.188034 48
   #S06 1.852288 40
   
   除了Faith’s PD，我们还可以从另一个角度去思考群落内物种的谱系多样性，即群落内任意两两物种的谱系关系有多近？即mean pairwise distance (MPD) 和mean nearest taxon distance (MNTD) 这两个谱系多样性指数 (Webb et al., 2002)。MPD表示群落内两两物种在系统发育树上距离的平均值，MNTD则表示群落内任意一个物种与群落中与自己亲缘关系最近的一个物种的平均谱系距离。picante包中的mpd和mntd函数可以计算上述两个谱系多样性指标。需要注意的是，与pd功能不同，mpd和mntd的输入数据格式是距离矩阵 (distance matrix) 格式而不是系统发育 (phylo) 格式，可以使用cophenetic命令将系统发育格式转化为距离矩阵格式。
   
   #MPD和MNTD谱系多样性
   MPD <- mpd(ant_com, cophenetic(ant_tree), abundance.weighted = FALSE)
   MNTD <- mntd(ant_com, cophenetic(ant_tree), abundance.weighted = FALSE)
   
   此外，为了与谱系多样性格局进行比较，我们可以将功能性状也转化为树状结构 (dendrogram)。本文先将功能性状转化为Gower相异距离矩阵 (Gower dissimilarity matrix) (Gower相异距离矩阵的优势在于可以处理分类变量的性状，并且可以为性状增加权重，但是本文并没有涉及上述两点，相关应用方法详见：Laliberté et al., 2014; Jarzyna et al., 2021)，这一步可以通过FD包中的gowdis函数实现 (Laliberté et al., 2014) (注意：在gowdis功能中，性状会被自动标准化转化，转化后所有的性状均值为0，标准差为1)；再使用UPGMA (unweighted pair group method with arithmetic mean) 方法将上述矩阵聚类，这一步通过hclust函数实现；最后使用ape包中的as.phylo函数将聚类的性状变为树状结构 (Paradis et al., 2004)，用于后续分析。
   
   library(FD)
   #Gower_trait = gowdis(ant_trait)
   #Hclust_trait = hclust(Gower_trait, "average")
   #trait_tree = as.phylo(Hclust_trait)
   
   由于使用树状结构可能影响功能信息的丢失，所以在进一步使用上述计算谱系多样性的方法计算功能多样性之前，我们先检测Gower相异距离矩阵和通过cophenetic命令转化功能树后得到的距离矩阵的均方离差 (mean squared deviation，即mSD) (Maire et al., 2015; Villéger et al., 2017; Jarzyna et al., 2021)。mSD的值等于0时，说明两个矩阵完美契合，值越大，说明两个矩阵的差异越大。
   
   #计算mSD
   dist.gow <- as.matrix(Gower_trait) # 将Gower相异距离矩阵转化为矩阵对象
   dist.dendr <- cophenetic(trait_tree) #通过cophenetic函数转化功能树后得到的距离矩阵转化为矩阵对象
   
   dist.gow[upper.tri(dist.gow)] <- NA #只保留下三角阵中的数据
   dist.dendr[upper.tri(dist.dendr)] <- NA #同上
   
   library(reshape2)
   
   dist.gow.m <- melt(dist.gow) #使用melt函数将上述下三角矩阵对象转化为数据框对象
   dist.dendr.m <- melt(dist.dendr) #同上
   
   dist.m <- cbind(dist.gow.m, dist.dendr.m [,3]) #按列合并上述数据框
   colnames(dist.m) <- c("sp1", "sp2", "gowerdist", "dendrdist") #重命名列名
   dist.m <- dist.m [complete.cases(dist.m),] #去掉NA值
   
   #根据Maire et al. (2015) 文中公式进行计算
   yst<- dist.m$dendrdist/max(dist.m$dendrdist) * max(dist.m$gowerdist)
   S <- 97 #物种数
   msd <- sum((dist.m$gowerdist-yst)^2)/((S*(S-1))/2)#0.042
   
   通过上述计算，得到mSD等于0.042，说明Gower相异距离矩阵和使用cophenetic转化功能树后得到的距离矩阵之间差异不大。当研究只考虑连续性状变量并且物种数和性状数量较小时，读者也可以选择使用原始的距离矩阵 (欧几里得距离和Gower距离等) 进行功能多样性的计算。注意：在性状 (例如大于10个) 和物种数很多的情况下，基于距离矩阵计算功能多样性是需要极大的计算力的，而使用树状结构计算功能多样性可以规避这一限制 (Maire et al., 2015; Zhao et al., 2021)。
   
   #功能多样性 (通过谱系多样性计算方法)
   FD_PD <- pd(ant_com,trait_tree)
   FD_MPD <- mpd(ant_com, cophenetic(trait_tree), abundance.weighted = FALSE)
   FD_MNTD <- mntd(ant_com, cophenetic(trait_tree), abundance.weighted = FALSE)
   
   2.3

   谱系和功能群落结构
   上述多样性的计算实际上是观察 (调查) 到的群落多样性格局，通过计算群落谱系和功能的标准化效应值 (Standardized effect size，SES)，可以得到谱系和功能的群落结构 (community structure) 的聚集或发散格局。标准效应值是通过观察 (调查) 值减去随机 (零模型) 得到的平均值的差，再除以随机得到的标椎差，公式如下：
   

   SES = (Obs - Mean_null)/SD_null

   其中Obs是实际观察到的谱系 (或功能) 多样性，Mean_null是随机多次 (通常为999次) 后得到的谱系 (或功能) 多样性的平均值，SD_null是随机多次后得到的谱系 (或功能) 多样性的标准差。可以通过picante包中的ses.mpd和ses.mntd函数实现。零模型有很多种选择，每种零模型都对应不同的前提假设，详见 (Gotelli, 2000; Swenson, 2014)，本研究选择使用“taxa.lables”零模型。该零模型的作用是将谱系树 (或功能树) 上物种的名字 (tip.label) 随机化，在保持群落数据的结构不变的情况下，进而模拟物种的谱系或功能随机分配的情况 (Swenson, 2014)，然后再使用物种名随机后的谱系树 (或功能树) 重新计算相应的多样性指数。
   SES值小于0表示群落内物种比预期的亲缘 (或功能) 关系更紧密，即群落聚集 (community clustering)；相反地，SES值大于0表示群落内的物种比预期的亲缘 (或功能) 关系更远，即群落发散 (community overdispersion)；如果SES等于0，则表示群落的谱系或功能多样性与零模型无差别，群落呈现随机格局。如果SES值小于–1.96，那么表示显著的群落聚集，如果SES值大于1.96，那么群落存在显著的发散格局 (Webb et al., 2002)。接下来，我们继续使用千岛湖33个岛屿上蚂蚁群落的数据进行演示：
   
   #谱系群落结构，即SES.MPD和SES.MNTD
   set.seed(10086) #在每次使用零模型前都设置一个seed (可以是任意数字)，保证之后重复运行的结果相同。
   result_mpd <- ses.mpd(ant_com,cophenetic(ant_tree), null.model = "taxa.labels", abundance.weighted = FALSE, runs = 100) #通常使用999次，这里为节省时间设置为100次
   set.seed(10086)
   result_mntd <- ses.mntd(ant_com,cophenetic(ant_tree), null.model = "taxa.labels", abundance.weighted = FALSE, runs = 100)
   head(result_mpd)
   
   
   
   结果中，ntaxa为每个岛屿的物种数，mod.obs为MPD的实际值，mpd.rand.mean为零模型得到的MPD的平均值，mpd.rand.sd为零模型得到的MPD的标准差，mpd.obs.rank为实际值在零模型结果中的分布位置，mpd.obs.z为SES值 (即SES.MPD)，mpd.obs.p为实际值在零模型结果中的分位数 (即P值，计算方式为mpd.obs.rank / runs + 1)，runs为零模型随机次数。
   
   #可以尝试使用上述的公式计算SES值，会得到相同的结果。
   (result_mpd$mpd.obs-result_mpd$mpd.rand.mean)/result_mpd$mpd.rand.sd
   
   同样的，我们将谱系树更换为功能树，可以计算功能群落结构:
   
   #功能群落结构，即SES.FD_MPD和SES.FD_MNTD
   set.seed(10086)
   result_fd_mpd <- ses.mpd(ant_com,cophenetic(ant_tree), null.model = "taxa.labels", abundance.weighted = FALSE, runs = 100)
   set.seed(10086)
   result_fd_mntd <- ses.mntd(ant_com,cophenetic(ant_tree), null.model = "taxa.labels", abundance.weighted = FALSE, runs = 100)
   
   2.4

   岛屿面积和隔离度对功能和谱系多样性及结构的影响
   那么千岛湖33个岛屿上蚂蚁群落的功能和谱系多样性 (以MPD和FD_MPD为例)，以及功能和谱系群落结构 (以SES.MPD和SES.FD_MPD为例) 随着岛屿面积和隔离度是怎样变化的呢？首先读取岛屿参数数据。
   
   #读取岛屿面积和隔离度数据
   area_and_isolation <- read.delim("area_and_isolation.txt", h=T, row.names=1, sep=",")
   
   #合并数据
   final_data <- cbind(as.matrix(result_mpd$mpd.obs), as.matrix(result_mpd$mpd.obs.z), as.matrix(result_fd_mpd$mpd.obs), as.matrix(result_fd_mpd$mpd.obs.z), area_and_isolation)
   colnames(final_data) <- c("MPD", "SES.MPD", "FD_MPD", "SES.FD_MPD", "LogArea", "Isolation")
   
   通过多元线性回归，判断岛屿面积和隔离度是否对多样性指标有影响。
   
   #多元线性回归
   summary(lm(MPD ~ LogArea + Isolation, final_data))
   
   
   
   summary(lm(SES.MPD ~ LogArea + Isolation, final_data))
   
   
   
   summary(lm(FD_MPD ~ LogArea + Isolation, final_data))
   
   
   
   summary(lm(SES.FD_MPD ~ LogArea + Isolation, final_data))
   
   
   
   #图形展示
   
   
   
   通过计算可以发现，岛屿面积而非隔离度，是影响千岛湖蚂蚁谱系和功能多样性的主要因子，并且随着岛屿面积的增加，谱系和功能多样性都升高，其中谱系多样性显著升高；此外，谱系和功能群落结构 (即SES.MPD和SES.FD_MPD) 也随着岛屿面积的增加而升高，说明小岛上的群落更聚集而大岛上的群落更发散。但是，值得注意的是，SES.MPD的值总体趋向于小于0，而SES.FD_MPD的值总体趋向于大于0，这说明谱系和功能可能反映出不同的生态过程：谱系群落结构反映出环境过滤过程，功能群落结构反映出竞争排斥过程，详见Zhao et al. (2020)。
   本文的研究基于alpha多样性，忽略了不同岛屿 (或群落) 间物种组成的差异。Beta多样性可以用来检验群落间组成的差异，相关研究方法可以参考赵维等 (2021) 的文章。
   

致谢

本文的研究案例来自论文 Zhao et al. (2020)。该项目受国家自然科学基金 (31500453、31872210、32101278、32101268) 和上海市科协启明星项目 (19QA1403300)资助。

参考文献

1. 赵维，王兴彪，侍浏洋，朱婉瑜，马磊，王敬敬，徐松，杨榕，张小霞，韩一凡，黄志勇. (2021). 原核微生物群落随机性和确定性装配过程的计算方法. Bio-101 e2003400. Doi: 10.21769/BioProtoc. 2003400.
2. Blomberg, S. P., Garland Jr, T. and Ives, A. R. (2003). Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution 57(4): 717-745.
3. Cavender-Bares, J., Kozak, K. H., Fine, P. V. and Kembel, S. W. (2009). The merging of community ecology and phylogenetic biology. Ecol Lett 12(7): 693-715.
4. Cadotte, M. W., Carboni, M., Si, X. and Tatsumi, S. (2019). Do traits and phylogeny support congruent community diversity patterns and assembly inferences? J Ecol 107(5): 2065-2077.
5. Diamond, J. M. (1975). Assembly of species communities. Ecology and Evolution of Communities, Belknap, Cambridge.
6. Davies, T. J. (2021). Ecophylogenetics redux. Ecol Lett 24(5): 1073-1088.
7. Fox, B. J. (1987). Species assembly and the evolution of community structure. Evol Ecol 1(3): 201-213.
8. Faith, D. P. (1992). Conservation evaluation and phylogenetic diversity. Biol Conserv 61(1): 1-10.
9. Gotelli, N. J. (2000). Null model analysis of species co-occurrence patterns. Ecol 81(9): 2606-2621.
10. Jarzyna, M.A., Quintero, I. and Jetz, W. (2021). Global functional and phylogenetic structure of avian assemblages across elevation and latitude. Ecol Lett 24(2): 196-207.
11. Jin, Y. and Qian, H. (2019). V.PhyloMaker: an R package that can generate very large phylogenies for vascular plants. Ecography 42(8): 1353-1359.
12. Kembel, S. W., Cowan, P. D., Helmus, M. R., Cornwell, W. K., Morlon, H., Ackerly, D. D., Blomberg, S. P. and Webb, C. O. (2010). Picante: R tools for integrating phylogenies and ecology. Bioinformatics 26(11): 1463-1464.
13. Kraft, N. J. B. and Ackerly, D. D. (2010). Functional trait and phylogenetic tests of community assembly across spatial scales in an Amazonian forest. Ecol Monogr 80(3): 401-422.
14. Laliberté, E., Legendre, P. and Shipley, B. (2014). FD: measuring functional diversity from multiple traits, and other tools for functional ecology. R package version 1.0-12. Vienna, Austria: R Foundation for Statistical Computing.
15. Maire, E., Grenouillet, G., Brosse, S. and Villéger, S. (2015). How many dimensions are needed to accurately assess functional diversity? A pragmatic approach for assessing the quality of functional spaces. Global Ecol Biogeogr 24(6): 728-740.
16. Mayfield, M. M. and Levine, J. M. (2010). Opposing effects of competitive exclusion on the phylogenetic structure of communities. Ecol Let 13(9): 1085-1093.
17. Münkemüller, T., Lavergne, S., Bzeznik, B., Dray, S., Jombart, T., Schiffers, K. and Thuiller, W. (2012). How to measure and test phylogenetic signal. Methods Ecol Evol 3(4): 743-756.
18. Pagel, M. (1999). Inferring the historical patterns of biological evolution. Nature 401(6756): 877-884.
19. Paradis, E., Claude, J. and Strimmer, K. (2004). APE: Analyses of phylogenetics and evolution in R language. Bioinformatics 20(2): 289-290.
20. Revell, L. J. (2012). Phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol 3(2): 217-223.
21. Swenson, N. G. (2014). Functional and Phylogenetic Ecology in R. Springer New York.
22. Villéger, S., Maire, E. and Leprieur, F. (2017). On the risks of using dendrograms to measure functional diversity and multidimensional spaces to measure phylogenetic diversity: a comment on Sobral et al. (2016). Ecol Let 20(4): 554-557.
23. Webb, C. O. (2000). Exploring the phylogenetic structure of ecological communities: an example for rain forest trees. Am Nat 156(2): 145-155.
24. Webb, C. O., Ackerly, D. D., McPeek, M. A. and Donoghue, M. J. (2002). Phylogenies and community ecology. Annu Rev Ecol Syst 33(1): 475-505.
25. Wickham, H. (2012). reshape2: Flexibly reshape data: a reboot of the reshape package. R package version 1.4-4. Vienna, Austria: R Foundation for Statistical Computing.
26. Zhao, Y., Dunn, R. R., Zhou, H. Si, X. and Ding, P. (2020). Island area, not isolation, drives taxonomic, phylogenetic and functional diversity of ants on land-bridge islands. J Biogeogr 47(8): 1627-1637.
27. Zhao, Y., Sanders, N.J., Liu, J., Jin, T., Zhou, H., Lu, R., Ding, P. and Si, X. (2021). β diversity among ant communities on fragmented habitat islands: the roles of species trait, phylogeny and abundance. Ecography 44(10): 1568-1578.