杜仲(Eucommia ulmoides Oliv.)是我国重要的国家战略资源、特有的名贵药材和木本油料树种, 其产业应用涉及医药、食品、保健、饲料、橡胶、化妆品以及园林绿化等领域^[1]。我国对杜仲的利用已有两千多年的历史, 成书于汉代的《神农本草经》视其为“上品”^[2], 对我国不同历史时期杜仲产地分布的研究认为, 早期杜仲的地理格局奠定了当代杜仲资源分布和产业化的基础^[3]。在地质史上, 杜仲属植物曾广布于北半球, 由于第四纪大冰期的侵袭, 在欧洲、北美等其他地区相继消失, 只在中国存活至今^[4-6]。随着杜仲用途与开发价值不断被认知, 杜仲栽植范围日渐扩大, 目前我国杜仲栽植面积在35万hm², 在国内27个省(区、市)均有种植^[7], 并成功引种到美国、英国、法国、匈牙利、俄罗斯、日本、韩国等国家^[8]。

对目标地区树种的潜在适生性进行评估, 有利于良种资源的合理布局、科学引种以及加快生产效益。对杜仲的适生分布以及产区划分已开展了一定研究, 认为我国杜仲中心产区在黔北、黔西北、鄂北、鄂西北、湘北、湘西北、豫西南、川东、川北、滇东北以及陕南地区, 主要产区包括中亚热带和北亚热带及南温带南部的局部分散产区, 在此以外的地方不宜栽植杜仲^[8-9]。但早期杜仲适生区划分多通过文献归纳, 缺少量化依据。生态位模型是根据物种已知的分布数据和相关环境变量, 通过一定的算法构建模型用以判断物种的生态需求, 并将运算结果投射至特定时间和空间来预测物种的实际分布和潜在分布, 被广泛用于全球变化背景下物种分布、生物多样性保护、物种潜在分布区预测、保护区规划、外来物种入侵区域的预测等^[10-11]。

常用的生态位模型有规则集遗传算法模型GARP(genetic algorithm for rule-set prediction model)、最大熵模型MaxEnt(maximum entropy model)、生物气候模型BIOCLIM (biological climatic model)以及气候动态模拟软件CLIMEX(match climate and compare location model)等^{[10, 12]}。其中MaxEnt模型是目前使用最为广泛的生态位模型^[13], 通过计算物种和环境互作系统中具有最大熵时的状态参数, 确定物种和环境间的稳定关系, 并以此估计物种的适生分布^[14]。MaxEnt模型能有效处理各变量间的复杂交互关系, 模型预测表现优异^[15], 已有研究将MaxEnt模型及空间分析技术手段应用至贵州省杜仲的生态适宜性区划中, 取得了具有较高可信度的分析结果^[16], 此外在珙桐(Davidia involucrata)、毛红椿(Toona ciliata var. pubescens)、油松(Pinus tabuliformis)、构树(Broussonetia papyrifer)、山茱萸(Corni fructus)、蓝莓(Vaccinium corymbosum)、枸杞(Lycium barbarum)等树种的适生区划中也得到了满意结果^[17-23]。本研究在对杜仲分布数据全面收集和整理的基础上, 利用MaxEnt模型及Arcgis空间分析技术预测我国杜仲的潜在适生区, 解析影响杜仲分布的主导环境变量, 以期为合理布局杜仲资源培育及产业发展提供理论依据。

1 材料与方法 1.1 物种分布数据来源及处理

杜仲在我国的地理分布信息通过查询中国数字植物标本馆(http://www.cvh.ac.cn/)、全球物种多样性信息库(GBIF, http://www.gbif.org/)、中国林业科学研究院经济林研究开发中心国家杜仲种质资源库调查档案以及整理相关文献获取, 分别得到1464条、222条、669条和1128条分布记录。通过排除经纬度缺失、无效以及重复登记信息, 并利用腾讯地图坐标拾取器校对(https://lbs.qq.com/tool/getpoint/)分布点数据, 获得727条有效地理分布记录。为降低人为调查偏向性造成分布记录的群集效应误差, 结合环境数据空间分辨率, 利用ArcMap软件设置10′×10′空间分辨率的格网, 在每网格中筛选距离中心最近的分布点, 最后得到有效分析样本地理分布记录601条。依照MaxEnt 3.4.1文件输入要求, 将数据整理为csv.格式备用。

1.2 环境变量数据收集及处理

从国家基础地理信息系统(http://nfgis.nsdi.gov.cn/)获取比例尺为1:4000000的中国行政区划图。环境变量共有54个, 包括19个气候变量、32个土壤变量和3个地形变量, 其中气候数据源自全球气候数据库网站(http://www.worldclim.org/, WorldClim 2.0)^[24], 空间分辨率为2.5′; 土壤数据和地形数据来自世界土壤数据库(http://www.fao.org/, HWSD 1.2)^[25], 空间分辨率为30″。根据MaxEnt 3.4.1模型分析对环境数据格式一致性的要求, 在Arcgis 10.2软件中将其空间分辨率统一重采样为2.5′。在物种分布模型的建模中, 环境变量的多重共线性将会使物种分布模型过度拟合, 导致预测结果精确性降低^[26]。利用Spearman秩相关分析分别考察气候变量、土壤变量中各变量之间的相关性, 相关系数矩阵见图 1—图 2, 保留相关系数小于0.75的变量, 对相关系数大于0.75的两个变量选择其中一个具有重要指示意义的环境变量。最终确定24个变量用于建模模拟, 如表 1所示。

1.3 模型预测及评估

将杜仲地理分布点和环境变量数据导入MaxEnt 3.4.1软件中建模运算, 随机选取75%分布数据为训练集, 剩余25%为测试集, 开启刀切法(jackknife)计算环境变量对杜仲分布的影响, 绘制各环境变量的响应曲线。设定模型运算500次, 设置10次重复, 其余选项采用默认设定^[27]。利用受试者工作特征曲线(receiver operating characteristic curve, ROC)对Maxent模型运算结果进行检验, ROC曲线以假阳性率(即1-特异度)为横坐标, 真阳性率(敏感度)为纵坐标, 曲线下的面积为AUC值。当AUC值介于0.70—0.80时, 认为预测结果较准确; 当AUC值介于0.80—0.90时, 认为预测结果很准确; 当AUC值大于0.90时, 认为预测结果极准确^[28-29]。

1.4 适生等级划分和主导环境因子分析

将Maxent分析得到的.asc文件在Arcgis 10.2中转化成Raster格式, 利用Arcgis 10.2空间分析工具重分类工具进行适生区划分。参照IPCC报告有关评估可能性的划分方法^[30], 并结合杜仲实际分布情况, 将杜仲潜在地理分布划分4个等级：0≤P≤0.1归为非适生区, 0.1＜P≤0.3归为低适生区, 0.3＜P≤0.6为中适生区, P≥0.6归为高适生区^[31]。将模型预测结果转化后导入ArcGIS 10.2绘制分布图, 并利用Raster Reclassify功能, 统计不同等级适生区分布面积。综合Jackknife检验、贡献率和置换重要值分析结果确定影响杜仲适生分布的主导环境变量, 根据杜仲潜在分布概率在不同主导环境变量中的响应曲线, 获得不同分布等级下主导环境因子的作用阈值。

2 结果与分析 2.1 Maxent模型的准确性评估

基于Maxent模型建模并重复运行10次后得出的AUC分析结果如图 3所示, 其中训练数据AUC值在0.898—0.903, 平均值为0.900;测试数据AUC值在0.861—0.907, 平均值为0.889。根据AUC值的评估标准^[32], 模型预测的准确性分别达到“极准确”与“很准确”的标准, 表明模型能很好的拟合杜仲物种分布数据, 其预测结果可信。

2.2 杜仲在我国的潜在适生区分布

基于MaxEnt模型预测的杜仲潜在适生区(包括高适生区、中适生区和低适生区)分布示意图如图 4所示, 潜在适生区面积总计为277.61万km², 约占国土总面积的28.92%。按照全国五级地貌等级分区体系^[33], 杜仲高适生区面积38.01万km², 集中分布在大巴山中低山谷地区、川东平行低山岭谷区、鄂西高原-大类山中低山丘陵谷地区、武陵山中低山谷地区、雪峰山中低山区、川南黔北滇东喀斯特高原中山区以及浙闽中低山丘陵谷地地区; 中适生区面积135.78万km², 集中分布在鲁东低山丘陵平原地区、武夷山中低山丘陵谷地区、江汉湖积冲积平原区、赣南低山丘陵盆地区、南阳盆地低山丘陵岗地平原区、汾渭洪积冲积平原台地区、盆西冲积平原区以及武陵山中低山谷地区; 低适生区面积103.82万km², 集中分布在燕山低山丘陵区、苏北湖积海积平原区、黄淮海冲积平原区、辽河下游冲海积平原区、鄱阳湖丘陵岗地冲积湖积平原区、幕阜山—罗霄山中低山丘陵盆地区、粤东低山丘陵平行岭谷区、桂西喀斯特低山丘陵盆地区、鄂尔多斯高原区及南横断山中高山区。

我国30个省级行政区中不同等级杜仲适生区分布面积预测结果如表 2所示, 以四川省杜仲潜在适生区面积最大, 达到23.092万km²; 四川省、云南省、湖南省、湖北省、河南省、贵州省、陕西省、广西壮族自治区、江西省、山东省、安徽省等地区杜仲潜在适生区分布较为集中, 约占全国总面积的68.32%。

2.3 影响我国杜仲分布的主导环境变量

利用Jackknife检验可以评估环境变量对预测结果的影响程度, 从而判断不同变量对杜仲潜在分布的重要性。Jackknife检验结果表明如图 5所示, 仅使用单独变量时, 正则化训练增益由高到低前10个环境变量依次为年平均温度(Bio1, 0.8534)、年均降水量(Bio12, 0.8238)、昼夜温差月均值(Bio2, 0.8025)、最暖季度降水量(Bio18, 0.7375)、温度季节性变化的标准差(Bio4, 0.5729)、降水量季节性变异系数(Bio15, 0.4182)、最暖月份最高温(Bio5, 0.3519)、最湿季度平均温度(Bio8, 0.3305)、海拔(Alt, 0.2974)、昼夜温差与年温差比值(Bio3, 0.2504);测试增益由高到低前10个环境变量依次为年平均温度(Bio1, 0.8648)、年均降水量(Bio12, 0.8354)、昼夜温差月均值(Bio2, 0.8105)、最暖季度降水量(Bio18, 0.7546)、温度季节性变化的标准差(Bio4, 0.5913)、降水量季节性变异系数(Bio15, 0.4321)、最暖月份最高温(Bio5, 0.3631)、最湿季度平均温度(Bio8, 0.3403)、海拔(Alt, 0.2994)、昼夜温差与年温差比值(Bio3, 0.2711)。贡献率和置换重要值也是评估各环境变量影响程度的重要指标, 由表 3可以看出, 年均降水量(Bio12, 53.8657)、年平均温度(Bio1, 24.1872)、温度季节性变化的标准差(Bio4, 11.1563)、海拔(Alt, 3.587)、坡度(Slp, 2.2436)5个环境变量的贡献率位居前列, 累积贡献率在95.0%以上; 置换重要值由高到低前5个环境变量依次为年平均温度(Bio1, 54.7322)、温度季节性变化的标准差(Bio4, 17.1309)、年均降水量(Bio12, 6.5668)、海拔(Alt, 6.5239)和昼夜温差月均值(Bio2, 6.0536)。

综合Jackknife检验、贡献率和置换重要值分析结果, 得出年平均温度(Bio1)、年均降水量(Bio12)、温度季节性变化的标准差(Bio4)、海拔(Alt)、坡度(Slp)是影响我国杜仲潜在适生区分布的主导因子。基于Maxent模型得出5个主导环境变量的响应区间如图 6所示, 以存在概率大于0.6作为杜仲最宜适生区的遴选条件, 5个主导环境变量的值域分别为：年平均温度介于11—16℃、年均降水量介于700—1450 mm、温度季节性变化的标准差在820以下、海拔在1400m以下、坡度在24°以下。

3 讨论

预估我国杜仲潜在的适生区分布, 对杜仲资源保护、引种栽培、精细化管理以及可持续利用具有重要意义。本研究利用601个杜仲地理分布点和22个环境变量数据, 通过MaxEnt模型对我国杜仲潜在适生区进行预测, 预测结果准确度在“很准确”的标准之上, 说明模型预测结果可信度高, 适用于我国杜仲适生区的估计与划分。与早期研究结果不同, 本研究认为杜仲高适生区分布还应包括鲁东低山丘陵平原地区、浙闽中低山丘陵谷地地区、武夷山中低山丘陵谷地区以及赣南低山丘陵盆地区, 这在实地考查以及一些文献报道中也初步予以证实^{[3, 34-37]}。根据MaxEnt模型及ArcGIS软件统计分析结果, 我国杜仲潜在适生区分布区面积达277.61万km², 其中高适生区面积达38.01万km², 这与目前国内杜仲栽植面积(约35万hm²)具有显著差距, 说明杜仲在我国各适生区引种栽植的潜力很大, 除四川、湖南、贵州、湖北、陕西等传统产区外, 今后在河南、山东、安徽、江西、浙江等省区也有巨大的发展空间。

综合Jackknife检验、贡献率和置换重要值分析结果, 年平均温度、年均降水量、温度季节性变化的标准差、海拔、坡度5个环境变量的重要性位居前列, 而12个土壤变量的贡献率之和仅在1%左右, 说明影响我国杜仲适生分布的主要环境因素为气候, 其次为地形, 而土壤因素的影响程度较小。温度和降水作为最主要的两个气候生态因子常被用来解释大尺度上植物的分布格局^[38]。年平均温度和年均降水量分别指征研究单元内重要的热量指标和水分指标; 温度季节性变化的标准差则主要反映研究单元内温度的季节分异, 用于指征气候的稳定性^[38-39]; 海拔包含了温度、水分、光照、土壤肥力等诸多环境因子的剧烈变化, 海拔梯度引起的环境异质性往往影响着植物的垂直分布格局^[40-41]; 坡度是水平方向上水分、土壤养分流的驱动因子, 对土壤厚度及其理化性质有显著影响, 进而影响着林木的定居和生长^[42-43]。本研究推算出杜仲最宜适生区5个主导环境变量的值域分别为：年平均温度介于11—16℃、年均降水量介于700—1450 mm、温度季节性变化的标准差在820以下、海拔在1400 m以下、坡度在24°以下, 这与部分杜仲生理生态研究所得出的结果相吻合^[44-46], 同时也进一步说明本研究所采用MaxEnt模型估算的合理性。