基于ArcGIS和R软件进行动物地理区划
Delineating Zoogeographical Regions Using ArcGIS and R Software   

研究背景

动物地理区反映了在历史发展过程中形成并在现代生态条件下聚集在一定地域的动物集合 (Morrone, 2015)。动物地理区的划分，可以了解不同动物区系的起源及其关系 (张荣祖，2011)，由此探讨物种的时空分布规律。因此它不仅为生物地理学、进化生物学和系统分类学等基础理论研究提供了清晰的空间框架 (Kreft and Jetz, 2010; 郜二虎等，2017)，而且对全球变化生物学的研究具有重要意义，如气候变化（Reygondeau et al. 2020），外来物种入侵（Capinha et al. 2015; Liu et al. 2020），疾病传播（Murray et al. 2015）等。
        阿尔弗雷德•拉塞尔•华莱士 (Alfred Russel Wallace) 在1876年提出了全球动物地理区划方案，并根据生物地理学上的发现，与达尔文共同创建了进化论。至今，该区划方案仍然是人们理解全球生物地理分布格局及分异过程的基础之一。然而，过去对动物地理区的划分主要依据博物学家或生物地理学家的经验，根据物种分布边缘、特有种分布或自然地理界线进行，定量分析手段不足，区划结果的可重复性较低。近年来，随着物种分布信息的完善、物种谱系关系的清晰、多元统计算法的发展，动物地理区划由定性研究转向定量研究，动物地理区划研究得到快速发展 (Kreft and Jetz, 2010; Holt et al., 2013; Vilhena and Antonelli, 2015)。目前，动物地理区划的研究范围已从区域尺度覆盖到全球尺度，研究对象由无脊椎动物到脊椎动物，从物种组成发展至谱系组成 (Kreft and Jetz, 2010; Rueda et al., 2010; Linder et al., 2012; Holt et al., 2013; He et al., 2017; 郜二虎等，2017)，这都展示出这一学科持续的生命力。
        ArcGIS和R是两个广泛应用于宏观生态学和生物地理学研究的软件。ArcGIS具有强大的空间数据采集、处理和分析能力，它在地理数据库构建、拓扑处理、地理坐标变换、数据可视化及制图等方面具有突出的优势。特别是内嵌在ArcMap中的ArcToolbox软件模块，具有许多高效实用的空间处理功能。R主要用于统计分析和制图，具有开源、免费、容易掌握和更新速度快等优势，而且收录了大量的软件包，为各专业领域提供了便利的交流平台。ArcGIS的矢量或栅格的空间数据可以利用R软件包” rgdal” (Bivand et al., 2021)读取和显示；R的统计分析结果可以输出shapefile矢量文件、csv属性文件或栅格文件并可直接在ArcGIS中读取。当今，运用ArcGIS和R进行动物地理区划已经形成一套较为成熟的技术方案 (Kreft and Jetz, 2010)，并被广泛应用于不同区域的动物地理区划研究中 (He et al., 2017; 郜二虎等，2017)。
本方案将详细介绍如何在ArcGIS中将物种分布范围 (面状数据) 转换为基于特定网格大小的物种分布数据。然后在R中调用相关软件包，计算网格间动物群落的β相异性，并通过排序和聚类分析进行动物地理区的划分。最后，根据区划结果进行空间可视化和制图。本方案经适当修改后可广泛应用于不同生物类群的地理区划中。

软件和相关软件包

1. ArcGIS 10.2 (ESRI Inc., Redlands, CA, USA)
2. R software version 4.0.4 (https://www.r-project.org/)
3. betapart version 1.5.3 (Baselga et al., 2012; http://cran.r-project.org/package=betapart)
4. dendextend version 1.15.1 (Galili, 2015; http://cran.r-project.org/package=dendextend)
5. phyloregion version 1.0.6 (Daru et al., 2021; http://cran.r-project.org/package=phyloregion)
6. picante version 1.8.2 (Kembel, 2020; http://cran.r-project.org/package=picante)
7. recluster version 2.9 (Dapporto et al., 2020; http://cran.r-project.org/package=recluster)
8. rgdal version 1.5-23 (Bivand et al., 2021; http://cran.r-project.org/package=rgdal)
9. vegan version 2.5-7 (Oksanen et al., 2020; http://cran.r-project.org/package=vegan)
   

分析案例

本方案将以哺乳动物作为研究对象，对青藏高原及毗邻地区的动物地理区进行划分。如今，环绕着青藏高原分别有蒙古高原、华北、华南、东南亚、南亚、西亚和中亚等物种组成各异、演化历史独特的动物地理区，形成全球动物区系最为复杂的区域之一 (He et al., 2020)。

分析流程

本方案分为4个部分 (图1)：

1. 在ArcGIS中将物种分布范围转换成网格分布数据。
2. 在R软件里读取和处理物种分布和系统发育树数据。
3. 在R中计算相异性系数、进行排序分析及系统聚类。
4. 在R中进行排序和聚类结果的空间可视化和制图。
   
   
   图1. 本方案的流程框架图
   
   

一、在ArcGIS中将物种分布范围转换成网格分布数据

物种分布信息是划分动物地理区的基础。物种分布数据可从国际自然保护联盟(International Union for Conservation of Nature, IUCN)物种红色名录的空间分布数据(https://www.iucnredlist.org/resources/spatial-data-download) 获得。截止至2021年，全球红色名录共评估了134,400多个物种，其中超过82% (> 111,000种) 有空间分布数据。此外，鸟类的分布数据可从国际鸟类联盟 (BirdLife International) 获得，爬行类的分布可从全球爬行动物分布数据 (Roll et al., 2017; https://doi.org/10.5061/dryad.83s7k) 获得。这些数据都是以”shapefile”矢量格式提供的面状数据 (polygon)，包括图形数据 (.shp)、图形索引格式 (.shx)、坐标投影信息 (.prj)和属性表数据 (.dbf)等。这里，仅以从IUCN红色名录获取的哺乳动物分布范围 (面状数据) 为例 (图2)，展示如何在ArcGIS中将物种分布范围转换成网格分布数据。


图2. 哺乳动物分布范围 (面状数据) 示意图。其中红色线范围代表研究区，灰色线为每个物种的分布范围，其属性表信息包括物种名、现今是否分布、起源和分布的季节4个字段。

1. 生成研究区等面积投影网格
   运行ArcGIS的用户桌面软件ArcMap，打开ArcToolbox工具箱模块，在Data Management Tools – Feature Class目录下运行"Create Fishnet"。按照图3所示设置相关参数，可以生成研究区范围内等面积投影网格，假定网格大小设为220 km × 220 km。需要注意的是，预设的研究区（图2中红色线范围）应作为"模板数据集" (Template Extent)，并已具有等面积投影的坐标系。
   
   
   图3. 运用ArcMap的Create Fishnet工具生成网格的参数设置
   
   
2. 将物种分布范围转换成网格分布数据
   打开ArcToolbox工具箱模块，在Analysis Tools – Overlay目录下运行"Spatial Join"。按照图4所示设置相关参数，可以将物种分布范围转换成网格分布数据。结果如图5所示。
   
   
   图4. 运用ArcMap的Spatial Join工具将物种分布范围转换成网格分布数据
   
   通过网格和物种分布数据的"Spatial Join"空间联合分析，新生成的图层将会把每个物种分布范围涉及到的网格赋予该物种的信息。如图5所示，网格号为”1”的网格与物种”sp1”和”sp2”的分布范围相交，结果中网格号为”1”的网格会出现2条记录；网格号为网格号为”5”的网格与物种”sp1”、”sp2”和”sp3”的分布范围相交，结果中该网格会出现3条记录，以此类推。
   
   
   图5. 物种分布范围转换成网格分布前后的数据结构变化。（a）用“"Create Fishnet”生成的网格，其属性表包括每个网格中心点的经纬度；（b）每个物种的分布范围，不同颜色代表不同的物种；（c）网格和物种分布范围进行“spatial join”后的结果，其属性表反映每个网格包括的物种。
   
   值得注意的是，以上在ArcGIS中的操作，R软件也提供了可选的解决方案。例如， R包”letsR”提供了将IUCN面数据转换成"网格 × 物种" (site × species)的功能”lets.presab”。但是该功能在R中运行速度较慢，如果我们的研究涉及大空间尺度成千上万物种的空间数据时，ArcGIS处理效率要比R高得多。由于篇幅关系，有兴趣的读者可以下载”letsR”的指南进行学习。
   
   

二、在R软件里读取和处理物种分布和系统发育树数据

# 打开R软件，安装以下函数包。
install.packages(c("picante", "betapart", "phyloregion", "dendextend", "rgdal", "vegan"))
# 加载函数包
library(betapart)
library(dendextend)
library(phyloregion)
library(picante)
library(recluster)
library(rgdal)
library(vegan)
进行分析前，需要把已有的数据 (附录1) 存放到自定义的工作目录，然后在R中设置工作目录并读取数据。
setwd("e:/R/method") #设置工作目录
data <- read.csv("spe_site.csv", head = T) #在ArcMap中输出的Spatial Join 结果
str(data) #可以使用"str()"函数了解对象"data"的数据结构
head(data) #也可以使用"head()"函数展示"data"数据集的前6行数据 (图6)


图6. 查看对象"data"的数据结构及前6行数据

把物种的网格分布数据转为"网格 × 物种" (site × species) 的物种组成数据。物种组成数据以物种的有/无 (1/0) 表示，每一行表示网格，每一列表示物种 (图7)。
comm <- xtabs(present ~ grid + species, data) #转换为"网格×物种"的群落数据


图7. 查看"网格 × 物种"的群落数据的结构

三、相异性系数计算、排序及系统聚类

1. 相异性测度
   β相异性反映不同群落间物种组成的相似 (相异) 性程度 (Whittaker, 1960)。在生物地理学和生态学中，学者们提出了多种β相异性指数的测度方法，这些测度方法对物种组成变化测度的内涵有所差异 (Koleff et al., 2003)。其中，辛普森指数 (Simpson) 是总体β相异性 (Sorensen) 的替换组分 (turnover)，受物种丰富度影响较少，能更好地反映物种组成变化的现实 (Baselga et al., 2007)，因此被广泛应用于动物地理区划中 (Kreft and Jetz, 2010; He et al., 2017)。以下以辛普森指数为例，利用"betapart"函数包中的"beta.pair"计算网格之间物种β相异性的距离矩阵。我们可以看到结果的数据集中有3个相异性矩阵，分别是"beta.sim"，"beta.sne"和"beta.sor" (图8)。
   taxa_beta <- beta.pair(comm, index.family = "sorensen") # 计算Sorensen的β相异性
   
   
   图8. 查看beta.pair计算结果的数据结构
   
   物种β相异性的距离矩阵是基于物种有/无来测度的，这个算法假设每一个物种都是相互独立的，物种的有无只存在有/无 (1/0)的关系。但事实并非如此，例如人与黑猩猩的亲缘关系就比虎的亲缘关系要近。系统发育β相异性是指两个群落之间的系统发育距离，可以反映不同群落间的亲缘关系上的相异程度，克服不同类群间物种、属、科中分类等级差异的问题。因此，随着陆栖脊椎动物系统发育树的不断完善(可从http://vertlife.org/phylosubsets/下载不同类群的系统发育树)，基于谱系关系的β相异性可以把物种间系统发育关系考虑进去，近年来被广泛关注及应用 (Holt et al., 2013)。
   本方案下载了哺乳动物的系统发育树 (Upham et al., 2019)，该系统发育树包含了5,911种哺乳动物 (包括107种灭绝物种) 的进化关系。由于部分物种的DNA序列数据缺失，这些物种在系统发育树上是基于形态学的属或更高的分类单元内随机内插而成，因此产生了10,000棵不同组合的系统发育树。我们随机下载了其中的1,000棵系统发育树，并提前计算出最大可信进化分支树 (MCC, Maximum clade credibility tree)。在这里，我们将根据MCC中物种间的谱系关系计算不同群落之间的谱系β相异性。需要注意的是，"betapart"包中提供的谱系β相异性计算方法’phylo.beta.pair’，运算速度较慢，为减少计算时间，我们在这里运用220 km × 220 km网格的群落数据。
   tree <- read.tree("mammal_Upham_2019_mcc.txt") #读取哺乳动物系统发育树
   subphy <- match.phylo.comm(tree, comm)$phy #对应群落和系统发育树中的物种
   phylo_beta <- phylo.beta.pair(comm, subphy, index.family = "sorensen") #计算谱系β相异性
   # 至此得到了基于物种和谱系的β相异性矩阵，下面可以选择其中一个进行聚类和排序分析。
   # sim <- taxa_beta$beta.sim #基于物种β相异性的Simpson距离矩阵
   # sim <- phylo_beta$phylo.beta.sim #基于谱系β相异性的Simpson距离矩阵
   
2. 系统聚类 (hierarchical clustering)
   系统聚类结果可以非常清晰地反映动物地理区之间的分层关系，这对动物地理区划十分重要。R中的"hclust"函数中提供了8种不同的系统聚类方法，可通过查看"?hclust"获取详细信息。在实际运用中，我们可以分别用8种方法进行聚类分析，然后通过计算聚类结果和原始距离矩阵之间的共表型相关系数 (cophenetic correlation coefficients) 选择最佳的聚类方法(Kreft and Jetz, 2010)。下面，我们以应用最为广泛的非加权配对算术平均法 (UPGMA) 进行聚类分析，并确定最佳的聚类簇的数量。
   hc <- hclust(sim, method = "average") #执行UPGMA聚类分析
   k <- optimal_phyloregion(sim, method = "average", k = 15) #计算最佳聚类簇的数量
   optimal_k <- k$optimal$k #获取最佳聚类簇的数量
   
3. 排序 (ordination)
   排序是一种常用的降维技术，在动物地理区划中可以将网格间物种组成的β相异性直观地投射到二维空间进行可视化，并以此分析动物地理区之间的相互关系。排序分析在大尺度的动物地理区划中得到广泛运用，例如Smith (1983) 根据哺乳类分析了Wallace (1876) 划分的全球24个亚区之间的动物区系关系；Williams et al. (1999) 分析了非洲热带鸟类的生物地理过渡区等。其中，与主成分分析 (principal components analysis, PCA) 或主坐标分析 (principal coordinates analysis, PCoA) 等方法相比，非度量多维尺度法 (NMDS) 不需要假设数据呈现正态分布或线性关系。因此，本方案将以NMDS为例介绍排序在动物地理区划的运用。
   #执行NMDS排序分析
   nmds <- monoMDS(sim) # sim为基于物种或谱系β相异性的Simpson距离矩阵
   nmds #查看NMDS分析结果
   
   
   图9. 检查NMDS分析结果
   
   从NMDS排序分析的结果中，我们可以根据应力值 (Stress value) 评估排序结果和原始相异性矩阵之间的拟合度，范围从0到1，该值越低则表明NMDS排序的结果越好 (图9)。
   #在二维空间展示NMDS结果 (图10)
   points <- nmds$points #
   col <- recluster.col(points) #将二维矩阵投影到RGB空间中，四个角分别为红、绿、黄和蓝。
   color <- rgb(col[,3], col[,4], col[,5], max=255) #根据RGB空间，定义研究区每个网格的颜色
   dev.new(width=80, height=40) #出图框大小为宽80，高40
   par(mfrow = c(1,2)) #出图框为2幅图，图版显示1行2列
   gird_shp <- readOGR(dsn = "e:/R/method/shp", layer = "grid") #读取研究区网格的shp图层
   plot(gird_shp, col=color) #将研究区的网格用RGB空间中的颜色展示在地图上
   recluster.plot.col(col, text=F, pch=19, cex=2, xlab = "NMDS 1", ylab = "NMDS 2")
   
   
   图10. 基于网格内物种组成相异性的NMDS排序图及空间分布
   
   

四、排序和聚类结果的空间可视化和制图
在动物地理区划中，空间可视化和制图是十分重要的步骤。上面的分析中，我们基于网格间物种组成的β相异性矩阵，分别计算得到排序和聚类的结果。这里，我们将介绍把两个分析结果在地图上可视化的方法 (图11)。
# 设置出图框参数
dev.new(width=150, height=50) #出图框大小为宽150，高50
par(mfrow = c(1,3)) #出图框为3幅图，图版显示1行3列
# 根据最佳聚类簇的数量设置颜色，并给聚类树状图 (dendrogram) 定义颜色
col <- colorRampPalette(c('gray','blue','green','yellow','orange','red','black'))(optimal_k)
dend <- as.dendrogram(hc)
color_dend <- color_branches(dend,k=optimal_k, col=col)
plot(color_dend)
# 根据聚类树状图的颜色，定义每个网格的颜色
clus <- data.frame(hc$order,leaf_Colors(color_dend))
color <- clus[order(clus$hc.order),]
colnames(color) <- c("grid", "col")
gird_shp <- readOGR(dsn = "e:/R/method/shp", layer = "grid")
plot(gird_shp, col = as.character(color$col))
# 根据聚类树状图和网格的颜色，定义NMDS散点图的颜色
x <- nmds$points[,1]
y <- nmds$points[,2]
d <- data.frame(row.names(as.matrix(sim)), x, y)
plot(d$x, d$y, col = as.character(color$col), pch=19, cex=2, xlab = "NMDS 1", ylab = "NMDS 2")


图11. 基于网格内物种组成相异性的聚类和排序分析的空间可视化

通过本方案，读者将能够使用ArcGIS软件处理从公开的数据库中 (例如IUCN、BirdLife等) 获得的大尺度物种分布数据，学会如何把物种分布范围转换成网格分布数据。此外，读者将能够运用R软件包计算不同网格间的物种 (或谱系) β相异性，进行排序和聚类分析，将网格分布数据转换为动物地理区划，并把动物地理区划结果进行空间可视化和制图。本方案的代码经过适当修改后，可以用于其他尺度和其他生物类群进行群落划分研究。

致谢

本研究受国家自然科学基金 (31900324)、广东省基础与应用基础研究基金 (2020A1515011472) 以及全国第二次陆生野生动物资源调查相关专题的资助，特此致谢。

参考文献

1. 郜二虎，何杰坤，王志臣，徐扬，唐小平，江海声. (2017). 全国陆生野生动物调查单元区划方案. 生物多样性 25: 1321-1330.
2. 张荣祖. (2011). 中国动物地理. 科学出版社. 北京. ISBN: 9787030319203.
3. Baselga, A., Jimenez-Valverde, A. and Niccolini, G. (2007). A multiple-site similarity measure independent of richness. Biol Lett 3: 642-645.
4. Baselga, A. and Orme, C. D. L. (2012). betapart: an R package for the study of beta diversity. Methods Ecol Evol 3: 808-812.
5. Bivand, R., Keitt, T. and Rowlingson, B. (2021). rgdal: Bindings for the 'Geospatial' data abstraction library. R package version 1.5-23.
6. Capinha, C., Essl, F., Seebens, H., Moser, D. and Pereira, H. M. (2015). The dispersal of alien species redefines biogeography in the Anthropocene.Science, 348, 1248-1251.
7. Dapporto, L., Ramazzotti, M., Fattorini, S., Vila, R., Talavera G. and Dennis R. H. L. (2020). recluster: Ordination methods for the analysis of beta-diversity indices. R package version 2.9.
8. Daru, B. H., Karunarathne, P., Schliep, K. and Silvestro, D. (2020). phyloregion: R package for biogeographical regionalization and macroecology. Methods Ecol Evol 11: 1483-1491.
9. Galili, T. (2015). dendextend: an R package for visualizing, adjusting, and comparing trees of hierarchical clustering. Bioinformatics 31: 3718-3720.
10. He, J., Kreft, H., Gao, E., Wang, Z. and Jiang, H. (2017). Patterns and drivers of zoogeographical regions of terrestrial vertebrates in China. J Biogeogr 44: 1172-1184.
11. He, J., Lin, S., Li, J., Yu, J. and Jiang, H. (2020). Evolutionary history of zoogeographical regions surrounding the Tibetan Plateau. Commun Biol 3: 415.
12. Holt, B. G., Lessard, J. P., Borregaard, M. K., Fritz, S. A., Araújo, M. B., Dimitrov, D., Fabre, P. H., Graham, C. H., Graves, G. R., Jønsson, K. A., et al. (2013). An update of Wallace’s zoogeographic regions of the world. Science 339: 74-78.
13. 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: 1463-1464.
14. Koleff, P., Gaston, K. J. and Lennon, J. J. (2003). Measuring beta diversity for presence–absence data. J Anim Ecol 72: 367-382.
15. Kreft, H. and Jetz, W. (2010). A framework for delineating biogeographical regions based on species distributions. J Biogeogr 37: 2029-2053.
16. Linder, H. P., de Klerk, H. M., Born, J., Burgess, N. D., Fjeldså, J. and Rahbek, C. (2012). The partitioning of Africa: statistically defined biogeographical regions in sub‐Saharan Africa. J Biogeogr 39: 1189-1205.
17. Liu, X., Rohr, J.R., Li, X., Deng, T., Li, W. & Li, Y. (2020). Climate extremes, variability, and trade shape biogeographical patterns of alien species. Current Zoology, zoaa068.
18. Morrone, J. J. (2015). Biogeographical regionalisation of the world: a reappraisal. Aust Syst Bot 28: 81-90.
19. Murray, K.A., Preston, N., Allen, T., Zambranatorrelio, C., Hosseini, P.R. and Daszak, P. (2015). Global biogeography of human infectious diseasesk. Proc Natl Acad Sci U S A, 112 (41), 12746-12751.
20. Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., O'Hara, R. B., Simpson, G. L., Solymos, P., et al. (2020). vegan: Community Ecology Package. R package version 2.5-7.
21. Reygondeau, G., Cheung, W.W.L., Wabnitz, C.C.C., Lam, V.W.Y., Frölicher, T. and Maury, O. (2020). Climate Change-Induced Emergence of Novel Biogeochemical Provinces. Front Mar Sci, 7: 657.
22. Roll, U., Feldman, A., Novosolov, M., Allison, A., Bauer, A. M., Bernard, R., Bohm, M., Castro-Herrera, F., Chirio, L., Collen, B., et al. (2017). The global distribution of tetrapods reveals a need for targeted reptile conservation. Nat Ecol Evol 1: 1677-1682.
23. Rueda, M., Rodríguez, M. Á. and Hawkins, B. A. (2010). Towards a biogeographic regionalization of the European biota. J Biogeogr 37: 2067-2076.
24. Smith, C. H. (1983). A system of world mammal faunal regions. I. Logical and statistical derivation of the regions. J Biogeogr 10: 455-466.
25. Upham, N. S., Esselstyn, J. A. and Jetz, W. (2019). Inferring the mammal tree: Species-level sets of phylogenies for questions in ecology, evolution, and conservation. PLoS Biol 17: e3000494.
26. Vilhena, D. A. and Antonelli, A. (2015). A network approach for identifying and delimiting biogeographical regions. Nat Commun 6: 6848.
27. Wallace, A. R. (1876). The geographical distribution of animals. Harper & Brothers. New York.
28. Whittaker, R. H. (1960). Vegetation of the Siskiyou Mountains, Oregon and California. Ecological Monographs 30: 279-338.
29. Williams, P. H., Klerk, H. M. and Crowe, T. M. (1999). Interpreting biogeographical boundaries among Afrotropical birds: spatial patterns in richness gradients and species replacement. J Biogeogr 26: 459-474.