研究背景

       生态学、生物地理学、进化和保护生物学的研究中往往需要分析生物个体的地理空间分布信息，这需要基于准确的物种分布记录。随着大量的生物多样性数据库的公开，研究者们能够获得的数据资源变得越来越丰富。例如，全球生物多样性信息网络(https://www.gbif.org/，Global Biodiversity Information Facility, GBIF)是世界上最大的生物多样性资源库，截止到2021年3月已收录1,657,726,342条物种分布记录。该数据库中的物种分布信息来自于博物馆标本、野外观察、文献资料等不同渠道，已被广泛应用于气候变化影响、生物入侵、物种丰富度格局、物种分布预测、环境生态位和人类健康等领域的研究。然而，这些数据库中的数据来源纷杂且数量庞大，其中存在的物种分类鉴定错误和有误的分布信息可能导致后续分析得出错误结果(Maldonado et al., 2015; Zizka et al., 2020)。因此，物种分布数据下载后必须进行数据过滤，以便获得较为可靠的分布点用于后续分析。除了物种分布记录，气候数据也被广泛应用于宏观生物学的研究以探讨生物进化、生态和生物地理学上的问题。它们可用于预测生物入侵、物种对全球变暖的敏感性、生态位进化对物种多样性的影响、气候驱动的物种形成和灭绝过程等。关于历史、现在和未来环境下的温度和降水的详细信息现在也可以在线获得。在过去的研究中，人们常常是基于手动搜索、下载和处理这些数据，需要使用到多个软件，耗费大量的时间，非常不方便，这就迫切需要有功能强大的软件来高效地完成这些工作。       R是一个广泛应用于科学领域的免费软件，能在其中进行各类数据管理、可视化、统计分析和建模，有大量的软件包可供使用，是进行不同生物类群空间分布分析的理想选择。rgbif包于2012年发布，该软件包可以从GBIF数据集中搜索物种分布数据，更快地获取数据，并进行可重复的研究(Chamberlain and Boettiger, 2017)。此外，最近发布的CoordinateCleaner包还可以检测和删除错误的或不精确的地理分布信息(Chamberlain, 2016; Zizka et al.,2019)，能够识别超大型数据集中有问题的记录。学习如何使用这些软件包获取和分析与物种分布相关的气候条件并探讨它们与物种性状(如体型)的相互关系，能够提升检验相关生物学假设的能力，并减少读者在学习多个软件上所花费的时间，显著提升效率。       准确的物种分布记录是解析空间分布规律分析的前提条件。本数据分析流程将详细介绍如何使用R软件从GBIF数据库中搜索物种的地理分布点记录(具有经度和纬度信息且是来源于博物馆馆藏标本记录的数据)，识别和删除有问题的分布记录(地理坐标位置录入错误、人为饲养、重复记录等)。此外，还将介绍如何从WorldClim数据库(https://www.worldclim.org/)中提取对应的物种分布点的气候变量。最后，将展示如何生成可视化图形，并初步评估物种分布和气候之间的联系。

软件和相关软件包

分析案例

       本方案将使用九带犰狳(Dasypus novemcinctus Linnaeus, 1758)作为研究案例。这个物种原产于新热带地区，在美洲有广泛的分布(从阿根廷延伸到墨西哥)。但在1850年左右成功入侵到美国(Feijó et al., 2018)。本方案计划分析这种犰狳在其入侵范围(美国)内所面临的气候条件是否与其原生分布区域相似(Feijó et al. 2020)。

分析流程

本方案分为四个步骤(如图1所示)：

图1. 本方案的流程示意图       掌握本数据分析方案最好有一定的R语言基础知识(对于R初学者，可以从阅读手册''R入门''开始：https://cran.r-project.org/manuals.html)。通过本方案读者能够使用"rgbif"包从GBIF网站下载目标种相关数据，并学会如何检测和过滤在线数据库固有的问题分布记录。同时，读者将学会通过个体分布的经纬度信息从Worldclim全球数据库中提取气候数据。该流程的脚本经过修改后，也可用于提取其它以栅格文件格式存储的气候、植被、土壤等生态因子。最后，读者可以学习如何解析数据和制图，对气候数据的分布和模式有初步的了解。#打开Rstudio, 安装以下软件包(图2). > install.packages(c("rgbif", "ggplot2", "CoordinateCleaner", "dplyr", "maptools", "raster", "ggcorrplot")) > library(maptools) > library(sf) > library(raster) > library(ggplot2) > library(ggcorrplot) > library(CoordinateCleaner) > library(rgbif) > library(dplyr)

图2. Rstudio截图示意

第一步：使用''rgbif'' 包从GBIF中搜索与目标物种(Dasypus novemcinctus)相关的数据。       在出版物中使用从GBIF获得的数据，需要使用与每个数据集相关联的数字对象标识符 (Digital Object Identifier，DOI)，并在文章中正确地引用它(更多内容请参阅https://www.gbif.org/citation-guidelines)。rgibf包中的函数"occ_download()"可以检索DOI，并且更适合于申请大数据使用权限，更多信息请访问https://docs.ropensci.org/rgbif/articles/downloads.html。在这种情况下，读者需要注册GBIF的帐户。       使用函数"occ_search()"从GBIF中检索物种分布记录，并设定参数(hasCoordinate = TRUE)来搜索仅基于“Preserved_Specimen”(即馆藏标本)的具有经纬度信息的数据。"occ_search()"函数的默认值是最多返回500条记录。研究者可以根据需求进行更改，每次查询的最大限制为100,000条分布记录。在本方案中，将搜索并获取1000条目标物种的分布记录。这首先需要检查基于博物馆标本的目标物种在GBIF中有多少地理分布记录；接着是获得GBIF网站中与目标物种相关的唯一物种识别号(key number)，该识别号可用于rgibf包的多个功能中。

> key key [1] 2440779 # Dasypus novemcinctus 的物种识别号。 > occ_count(taxonKey=key, georeferenced=TRUE, basisOfRecord= "PRESERVED_SPECIMEN") [1] 1584 # 目前，GBIF数据库中有1584条基于博物馆的记录是具有经纬度信息的。需要注意的是，因为数据库一直在更新，这个数字可能会变动，因此使用本流程过程中可能会有部分数值有一些小范围的变动。

       读者可以通过定义特定的参数来定制查询，例如，仅从某个国家、或者基于经纬度限制的某个地区，或者特定的海拔(如1000m以上)。其他选项还包括：大陆、月份、年份等。在Rstudio中输入"?occ_search()"以了解可以选择的所有参数(图3)。

图3. 学习如何使用"occ_search"函数自定义GBIF搜索

       使用参数(hasCoordinate = TRUE)搜索具有经纬度信息的记录，并且只基于有标本的信息"preserved_sample"。> data str(data$data) # 这个命令将检查数据的整体结构(图4)。

图4. 从GBIF获取数据的原始结构

       str的输出显示数据包括 1000行(每一行代表一个单独的记录)和152列，每个记录包括的信息有国家名称、纬度和经度、海拔、采集时间、博物馆标本馆藏情况等)。根据科学问题，可以用多种方式解析这些数据。> names(data$data) #检查数据的列名。 > head(data$data[,c("countryCode","country")]) # 例如：检查54列和57列的列名(图5)。

       基于本方案的目的，只需要保留数据表中包含“物种名称”、“经度”、“纬度”、“国家名称”和“记录的不确定性(以米为单位)”(这将允许评估地理坐标的准确性)等数列。为了后续分析的方便，建议用简短好记的名称重命名各列，并使用"dplyr"R包将"species"和"country"作为因子。> datasel%                       dplyr::select(species, decimalLongitude, decimalLatitude,country,coordinateUncertaintyInMeters) %>%          rename(Species = species, Long = decimalLongitude, Lat = decimalLatitude, Uncertain=coordinateUncertaintyInMeters) %>%          mutate (Species = as.factor(Species), country= as.factor(country)) > datasel # 现在有了一个包含1000行、5列的data.frame(图 6)。

       读者已经可以从这个数据集中解析一些信息，可以看看每个国家按降序排列有多少条记录。 > datasel %>% group_by(country) %>% summarize(total= n()) %>% arrange(desc(total)) #图 7

       为了分析目标物种的原生地和入侵地的气候差异，需要创建一个新的列("范围'')，并将记录分为"原生"和"入侵"两类。同时，删除没有国家信息的记录(如上面只有一个记录显示为"NA"的列)。 > datasel% filter(!is.na(country)) %>%          mutate(Range = as.factor(if_else (country == "United States of America", "Invasive", "Native")))       现在，有了名为"datasel"的对象，其中包含从GBIF中提取的Dasypus novemcinctus分布记录。接下来将检查九带犰狳在原生地和入侵地范围内各有多少条分布记录。 > datasel %>% group_by(Range) %>% summarize(total= n()) %>%arrange(desc(total)) #图 8

图8. 原生和入侵范围的九带犰狳分布记录的数量

第二步：使用CoordinateCleaner包删除有问题的分布记录。       由于使用的数据来自GBIF，因此需要仔细检查并删除有问题的分布记录(Zizka et al，2020)。读者可以首先通过将分布记录的空间分布情况可视化地展示在地图上，直观地检查问题数据所在位置。分布记录点的可视化可以通过"maptools"包来实现。例如，如果该物种是一个国家特有的，就可以删除该国以外的分布记录。

> data(wrld_simpl) #加载世界地图。 > plot(wrld_simpl, col="tan2",bg="lightblue", border= NA) #图示地图。 > points(datasel$Long, datasel$Lat, bg="red", col='white', pch=21, cex=1.2) #将从GBIF下载的数据映射到地图上 (图 9)。

图9. GBIF检索到的九带犰狳的分布记录点

    从目标物种的物种分布记录图可以看出，大多数记录都符合已知的物种分布范围(从阿根廷到美国)，但海洋中出现了经纬度=0的分布点，说这明数据本身存在错误记录。使用以下代码检查有多少记录的经纬度等于零。 > datasel %>% filter(Long==0, Lat==0) # 统计经纬度等于零的记录 (图10) 。

图10. GBIF搜索到的21条记录(行)的经纬度为零

       正如读者所看到的，尽管country列显示的是"Mexico''，但有21个个体的经纬度信息被列为零。在R中有一些包可以帮助读者检查最常见的问题(例如scrubr, biogeo)。rgibf包还有一个函数"occ_issues()''，可以帮助检测和删除有问题的记录。这里将使用"CoordinateCleaner"包(Zizka et al. 2019)。在这个包中有几个可用的参数(参见Zizka et al. 2019中的表1)，可以通过在Rstudio中输入"?clean_coordinates()"来了解更多的参数选择。本教程将检测记录点坐标是否围绕首都、国家的中心，是否落入海洋，为零，或在饲养动物的博物馆(机构)周围。读者也可以使用世界自然保护联盟(International Union for Conservation of Nature，IUCN)的物种分布范围来剔除分布范围以外的所有记录。> problem_records summary(problem_records) # 46个有问题的分布点被发现(图 11).

图11. 由CoordinateCleaner包标识的问题记录摘要的屏幕截图

       图示出问题记录的分布地点。> datasel_problematic plot(wrld_simpl, col="tan2", bg="lightblue", border= NA) # 显示地图。> points(datasel_problematic$Long, datasel_problematic$Lat, bg="red", col='white', pch=21, cex=1) # 将有问题的点映射到地图上 (图 12)。

图12. GBIF检索到的九带犰狳有问题的分布记录点

       接下来将简单地排除所有有问题的记录。但在读者的研究中，建议仔细检查它们，特别是当研究物种的分布记录数量有限时。由CoordinateCleaner包标记的一些记录可能是正确的。此外，读者的数据集中可能还存在与分类错误有关的其它问题(例如，标本鉴定错误等)。数据过滤过程需要更仔细地查看整个数据集。使用R包的自动检测应该是检测和删除不准确记录的步骤之一，而不是唯一的步骤。       从数据集中排除所有有问题的记录。

       看看有多少个体有这样的信息，并把值的范围图像化。> summary(is.na(datasel_clean$Uncertain)) > hist(datasel_clean$Uncertain/1000, breaks = 100, main="Coordinate uncertainy in km", xlab= "Uncertainty in km") # 注意，GBIF提供的信息是以米为单位的，但是为了评估不确定性，使用千米为单位更加方便。 (图13)。

图13. GBIF检索到的分布记录的不确定程度(以km为单位)

       直方图显示，只有少数地方超过100公里(~1度)。在这里，将简单地删除不确定性超过100公里的记录。但如果有更精确的坐标信息便可以减少这种不确定性。 > datasel_cleanF% filter(is.na(Uncertain)| Uncertain/1000 str(datasel_cleanF) # 图 14。 图14. 过滤后的物种分布记录条数(688行6列)        经过过滤和去除重复步骤，最终的数据集共有688条。这些记录是基于分布在该物种原生和入侵范围的博物馆标本分布信息。以下代码可以图示这个最终数据集。> plot(wrld_simpl, col="tan2", bg="lightblue", border= "white", xlim=c(-132,-70), ylim=c(-58,69)) #这里可以通过设置经纬度的限制而只显示美洲的分布点。 > points(datasel_cleanF$Long, datasel_cleanF$Lat, bg="black", col='white', pch=21, cex=1) # 图15。 图15. 九带犰狳经过数据过滤后的分布记录

第三步:使用"Raster"包从WorldClim下载环境数据。        这一步是使用 "Raster"包从 WorldClim下载环境数据。网上有很多关于全球温度、降水、土壤、植被等信息的资源，收集了历史、当前和未来不同时期的气候数据。最常用的气候数据来源是WorldClim，这是一个全球气候数据库，包含19个温度和降水变量，以及太阳辐射、风速和水汽压。每一个变量图层都有四种空间分辨率：10 min、5 min、2.5 min和30 s。其中30 s代表大约1 km2网格，10 min代表大约340 km2网格。更多信息请点击http://www.worldclim.com/version2。本方案将使用R中的“Raster”包下载分辨率为5 min的19个气候图层。 > bioclim names(bioclim) # 图16

图16. 存储在"bioclim"对象中的不同气候变量层的名称

       接下来可以使用"plot"功能轻松地将这些图层中的信息可视化。需要注意的是，要从WorldClim数据中获得实际的气温值(°C)，需要将这些值除以10。> plot(bioclim$bio1/10, main="Annual Mean Temperature") # 气候数据库下载的年均温原始数据需要除以10，以便获得真实的数值, 图 17。 图17.全球年均温分布图(°C)> plot(bioclim$bio12, main="Annual Precipitation ")# 年均降雨， 图18。 图18. 全球年降水量分布图        还可以通过更改"zlim"来设定图中显示的值的范围。例如，如果只显示全球最冷季均温高于 0 °C的地方，可以将zlim设置为0。> plot(bioclim$bio11/10, zlim=c(0,28.9), main="Winter above 0°C") # 图19。 图19. 全球最冷季气温分布图(只显示值高于0 °C的地区)> plot(bioclim$bio11/10, zlim=c(-49.7,0)) # 显示温度低于零度的地方。图 20。 图20 . 全球最冷季气温分布图(只显示值低于0 °C的地区)

       现在，可以使用栅格包中的函数''extract()''从目标物种的每个分布记录中提取当地的环境数据。注意，这段代码可以应用到栅格格式的任何图层。如果读者已经有一个栅格图层保存在电脑中，使用函数"raster()"能够加载它到R中进行分析。可以在''?raster()"中读取更多参数信息。        首先，从GBIF数据集中只选择经度和纬度列。R有一个特定的顺序，经度优先。> longlat data_clim datasel_clim head(datasel_clim) #这个表格现在同时包含了物种记录和环境数据， 图 21。 图21. 数据摘要截图> write.csv(datasel_clim, "Dnomcinctus_GBIFrecords_climate.csv") # 保存表格，表格会被保存到当前工作目录，getwd()命令可以查看当前工作目录，而setwd()可以设定当前工作目录。 第四步: 解析数据并初步构建图表        一旦收集了每个地区的气候数据，就可以用很多方式来研究它。例如，可以利用这些信息来检验气候对体型的影响。在本教程中，感兴趣的是比较九带犰狳的原生分布区和入侵分布区的气候条件是否存在差异。        接下来可以从一些基本的问题开始，例如原生和入侵犰狳面临的平均、最低和最高年平均温度(AMT: bio1)有怎样的差异?(图 22)> datasel_clim % filter(!is.na(bio1)) #为了确保本流程中的分析顺利完成，需要删除标志为“NA”的数据，因为在GBIF数据中，有部分记录的经纬度信息登入有误。在读者自己的研究中，需要仔细核对这些数据。 > datasel_clim %>% group_by(Range) %>% summarize(meanAMT=mean(bio1/10),minAMT=min(bio1/10),maxAMT=max(bio1/10)) 图22. 九带犰狳入侵地和原生地分布区年均温、最低温度和最高温度的差异

       可以看到入侵区域的年均温度远低于原生区域，还可以通过更改上面代码中的"bio"对象来探索在其它气候变量上的差异。        接下来，使用"ggplot2"软件包中的一些图表，在原生和入侵范围数据中更多地分析犰狳所面临的气候条件。ggplot是一个广泛用于构建图的包。读者可以在ggplot书中阅读更多内容：https://ggplot2-book.org/index.html        从简单的箱线图和散点图开始(图 23和24)：# 使用箱线图显示年均温(bio1) > ggplot(datasel_clim, aes(Range, bio1/10)) + geom_boxplot(fill="tan2") + theme_classic() + labs(title= "Annual Mean Temperature",     x="", y= "Temperature in C") # 使用箱线图显示平均气温日差(bio2) > ggplot(datasel_clim, aes(Range, bio2/10)) + geom_boxplot(fill="tan2") + theme_classic() + labs(title= "Mean Diurnal Range",         x="", y= "Temperature in C") # 使用箱线图显示气温的季节性(bio 7). > ggplot(datasel_clim, aes(Range, bio7/10)) + geom_boxplot(fill="tan2") + theme_classic() + labs(title= "Temperature Seasonality",         x="", y= "Temperature in C") # 使用箱线图显示年降雨(bio12) > ggplot(datasel_clim, aes(Range, bio12)) + geom_boxplot(fill="tan2") + theme_classic() + labs(title= "Annual Precipitation",         x="", y= "Precipitation (mm)") 图23. 九带犰狳原生和入侵分布区温度和降水比较# 使用年均温(bio1)和年降水(bio12)制作散点图 > ggplot(datasel_clim, aes(x=bio1/10, y=bio12)) + geom_point(aes(col=Range), size=2) + theme_classic()+ labs(title= "Scatterplot",         y="Annual Precipitation (in mm)",         x="Annual Mean Temperature (C)") 图24. 九带犰狳原生和入侵分布区年均温和年降水量的散点图        一些气候数据是高度相关的，这可以很容易地用"ggcorrplot"软件包来可视化从WorldClim获得的19个气候变量之间的Pearson相关性(图 25)。> clima_cor ggcorrplot(clima_cor,type = "lower",lab_size = 3,method="square",                       hc.order = F, lab=T,                       title="Correlogram of WorldClim") 图25.九带犰狳原生和入侵分布点19个环境变量之间的皮尔逊相关性检验(暗红色代表高度正相关，深蓝色代表高度负相)。        读者可以看到很多变量都是高度相关的。最常用的方法是采用排序技术，如主成分分析 (PCA)。 # 首先选择包含气候数据的列 > climatedata summary(climatedata) # 图 26. 图26. 截图显示与物种发生相关的气候值的基本信息(最小值、中值、平均值、最大值)        注意，使用"summary"函数分布的数据，可以看到每个变量有不同的限制(最小和最大)。例如，bio1的最大值为287(代表28.7 C)， bio4的最大值为10237。这导致在主成分分析或其他统计检验时，与每个性状不同值范围相关的数学偏差。为了避免这个问题，可以对变量进行缩放，使其范围标准化，具有可比性(图 27)。> climate_scale summary(climate_scale) # 注意所有变量的均值现在都等于0 图27. 标准化后的与物种分布记录相关的气候值汇总(最小值、中值、平均值、最大值) > plot(climatedata[,1], climate_scale[,1])# 注意标准化后的数据与原始数据是完全相关的. > cor(climatedata[,1], climate_scale[,1]) # Pearson 相关性 = 1        现在可以使用"princomp"函数执行主成分分析(Principal Component Analysis， PCA，图28)。读者还可以跳过预先的缩放转换，并在运行princomp函数时添加"cor = TRUE”。 > pcaclimate summary(pcaclimate) # 显示每个维度上能够解释的变异，例如本案列中PCA1解释了~51.2%的变异，而PCA1和PCA2总共解释了71.6%的变异. 图28. 各主成分贡献值        PCA的缺点是每个主成分都与原始变量有不同程度的关联。以下柱状图(图29)显示每个变量与第一主成分(PCA1和PCA2)的关联强度。 > barplot(pcaclimate$loadings[,1]) # 条形图显示PCA1 > barplot(pcaclimate$loadings[,2]) # 条形图显示PCA2 图29. 19个气候变量与第一(上)和第二(下)主成分的相关性        通过绘制前两个主成分得分可以直观地看到九带犰狳在原生和入侵地区的气候条件分布概况(图30)。> plot(pcaclimate$scores[,1], pcaclimate$scores[,2], xlab="PCA1", ylab="PCA2") 图30. 第一和第二主成分散点图        但是，R作图有一定的局限性。为了制作更好的图表或在以后的测试中使用第一个PCA得分值作为新的分析变量，可以将得分值保存在一个新的列中，并将其合并到原始表中。> PCA1 PCA2 PCA3 datasel_climPCA str(datasel_climPCA) # 正如读者看到的这样，所有的信息现在都整合在一个表格里面，现在通过'write.csv' 函数将它保存下来。        为了构建一个漂亮的PCA图形，再次使用ggplot2包(图31)。 > ggplot(datasel_climPCA, aes(x=PCA1, y=PCA2, color= Range))+ geom_point(aes(color=Range), size=2)+stat_ellipse()+ labs(x= "PCA 1", y= "PCA 2")+ theme_classic) 图31. 原生和入侵九带犰狳分布地点总体气候差异散点图        最后，可以使用"ggsave"函数保存这个图。读者可以保存多种格式(例如:pdf, tif, jpg, png)。通过查看" ?ggsave() "可以获取更信息。> ggsave("PCA_armadillo_bioclim.pdf") # 保存为pdf格式 > ggsave("PCA_armadillo_bioclim.tiff", width = 15, height = 15, units = "cm", dpi= 300) # 保存为Tif格式       总的来说，这些初步的分析结果表明：原生地的九带犰狳与入侵地的九带犰狳栖息环境已存在明显的分化。在实际应用中，基于对数据的审核，以上的分析还需要进一步深化，如使用t检验进行两组间均值差异比较；使用方差分析(ANOVA)进行多组之间的比较等；使用ENMeval包进行生态位建模，对比研究不同目标类群潜在栖息地分布范围上存在的差异等。读者可以根据研究需要进行进一步的学习。

致谢

研究团队受到第二次青藏高原科学考察相关专题的资助(No. 2019QZKK0402, 2019QZKK0501)；Anderson Feijó博士受到中国科学院国际人才计划，国际博士后奖学金(2021PB0021，2018PB0040)的资助，特此致谢。

参考文献