气候变化是当前人类所面临的全球性环境问题[1]。联合国政府间气候变化专门委员会(Intergovernmental Panelon Climate Change, IPCC)评估报告显示：到21世纪末, 全球地表平均气温将升高0.3—4.8℃[2], 降水格局也会发生明显变化[2-3]。当气候改变时, 大多数物种会受到生理胁迫从而影响个体生长、种群动态，并最终导致物种分布区的消长、变迁[4]。最为显著的证据是, 在更新世末期, 冰期和间冰期交替, 气候波动剧烈, 引起地球上生物的大规模迁移, 甚至灭绝[5], 对现今的生态系统格局产生重大影响。因此, 研究气候变化下物种分布格局的变化, 不仅对物种的保护和资源利用至关重要, 而且对未来生态系统的可持续管理也具有重要意义。
物种分布模型主要利用物种的现代分布记录与环境数据, 依据特定的算法构建物种的生态位, 并投射到相应的气候情景中, 以概率的形式反映物种对生境的偏好程度, 从而估算出物种的潜在分布区域[6-7]。大量研究证明其对物种分布具有非常好的预测能力, 已广泛用于探究历史避难所、引种适生区以及入侵物种扩散地等研究中[6-9]。随着气候变化预估可靠性的大幅度提高, 将物种分布模型和未来气候数据相结合, 预测物种在未来全球气候变暖背景下的分布格局变化正成为生物地理学和全球变化生态学研究的热点。例如, 周佩和潘石玉等分别通过MaxEnt模型对费菜和何首乌在气候背景下的适生区分布格局研究发现研究物种的适生区范围均有缩小, 生境适宜度降低, 得出物种分布的稳定适生区, 这为物种的现阶段采用、保存及可持续开发利用提供了参考[10-11]。Rana通过对两种百合科濒危药用植物的分布研究发现两物种适生区向西北移动, 这对进一步开展两种药用植物的保存具有重要的指导意义[12]。但也有研究者指出, 物种的分布格局不仅受到气候因子的影响, 也受到地形、土壤、生物间相互作用等多因素的共同制约。目前, 在物种分布模型中, 大多只考虑气候因子[13-14]。因此, 此类模型预测的结果仅反映了物种存续与气候条件的匹配度, 而忽略了其他因子对物种分布的作用, 可能造成过度预测[8, 15-16], 所以, 在利用该模型进行物种分布预测时需考虑其他因素对预测结果的影响。
在自然界中, 物种常常只在特定的空间范围内分布, 即使部分区域内气候等生境条件适合, 物种也不能在这些区域存在。研究者将物种在空间内有限分布的现象称为空间约束并认为这在很大程度上与物种的扩散能力相关, 而物种扩散主要受到地理屏障阻隔、自身扩散特性和扩散时间的影响[15-19]。近年来, 研究者通过趋势面法、主轴邻距法等对物种分布区域的栅格中心点的经纬度坐标进行运算, 从而形成物种分布的空间约束因子变量[18]。研究表明, 在物种分布模型中加入空间约束因子变量后, 其模型预测效果优于单纯的气候变量模型[8, 18], 更为重要的是能够在很大程度上消除部分过度预测区域, 其结果符合物种迁移、分布的客观实际。近年来, 在预测未来物种分布时, 空间约束因子变量在国外得到了广泛的应用。例如, 一些欧美研究者以空间约束因子来反映植物的迁移作用, 并在预测未来气候变化下全球、各大洲以及生物多样性热点地区的植物分布范围、物种多样性和丰富度变化的研究中大量运用[20-22]。其结果表明由于物种的扩散速度较弱, 考虑空间约束因子后物种地理范围变化更加连续、集中, 有效去除了单纯气候因子导致的过渡预测的理论分布区, 这对深入认识当地植物对气候变化的响应程度以及预测生态系统服务功能变化具有重大意义[20-22]。此外, 对于移动、迁移速度较快的动物, 前人研究也表明其也受到空间约束因子的限制, 不能完全占据理论适生区[8, 12]。相比之下, 国内在预测物种未来潜在分布的报道中, 较少考虑空间约束作用的限制, 可能导致预测结果有所偏差[23]。
长苞铁杉(Tsuga longibracteata)是我国特有的古老孑遗植物[24-25], 属易危种, 主要集中分布于广西、湖南、贵州交界处以及福建南部等中亚热带地区, 历史分布研究较少, 孢粉数据未鉴定到种, 故对其历史分布情况尚不明晰。长苞铁杉在植物界地位独特, 对研究裸子植物的系统发育、古生态和古气候均具有重要价值[24], 同时也是珍贵的造林、用材树种。目前对长苞铁杉的研究主要集中于种苗更新[26-27]、年龄结构[24]、种间竞争[28]等方面。长苞铁杉对气候变化较为敏感, 已有研究证实气候变化会导致局部范围内的长苞铁杉分布区缩小, 居群濒临灭绝[24]。然而, 在未来气候变化背景下, 全国尺度下的长苞铁杉的空间分布变迁情况尚不清楚, 制约了对长苞铁杉的保护和可持续管理。
本文基于MaxEnt模型预测长苞铁杉在不同时期和气候变化情境下的潜在地理分布区, 探明其分布格局的变化趋势, 揭示引起潜在地理分布改变的关键因子和关键保育区域, 并比较将空间约束因子加入物种分布模型后对预测分布区的限制作用, 研究结果将为长苞铁杉的保护和可持续利用提供重要科学依据。
1 材料与方法
1.1 数据来源和预处理
通过全球生物多样性数据库(GBIF, https://www.gbif.org)、中国数字植物标本馆(http://www.cvh.org.cn)、国家标本平台(http://www.nsii.org.cn)、公开发表的中英文文献等收集长苞铁杉的现实分布点位数据(未收集到中国台湾的点位信息), 并在部分区域开展野外核查。对部分未标明地理坐标的标本数据, 按照采样地点描述, 通过谷歌地球软件进行分布点人工校正, 验证数据点位信息。随后, 将分布点在ARCGIS 10.3中剔除重复项(5 km×5 km网格中仅保留距中心点最近的分布点), 最后获得参与模型运算的采样点, 共计76个(图 1)。
图 1 长苞铁杉现代分布点位
Fig. 1 Distribution points of Tsuga longibracteata
图选项
本文的气候变量数据均来源于世界气候数据库(http://www.worldclim.org/), 分别选取现代(1960—1990年)和未来(2050年和2070年)RCP2.6和RCP8.5情境下的气候数据(目前此类模型模拟研究大多采样1960—1990年间的气候代表现代气候因子)[8, 10, 13]。其中, RCP2.6代表全球变暖幅度控制在不超过工业化以前的温度2℃的情景；RCP8.5代表温室气体较高排放的情景, 全球变暖趋势愈加显著, 温度增加超过2℃。每个时期的气候数据包含19个变量, 空间分辨率为30″。
考虑到环境变量的多重共线性会导致模型的过度拟合[19], 在构建模型之前, 参考Yinbo Zhang对中国气候变化下的植物分布研究, 对19个气候变量进行相关性分析, 在充分考虑气候因子生物学意义的基础上, 去除相关性大于0.75的气候因子[29], 仅保留不相关变量作为影响长苞铁杉分布的气候因子[19, 30], 降低生态空间维度, 提高模型转移能力[6]。经筛选得到9个气候变量(表 1)。
表 1 用于运算的环境变量
Table 1 Environmental variables in MaxEnt model
代码Code
环境变量Environmental variables
Bio-5
最暖月最高温Max. temperature of warmest month
Bio-6
最冷月最低温Min. temperature of coldest month
Bio-7
温度年较差Temperature annual range (Bio-5—Bio-6)
Bio-9
最干季均温Mean temperature of driest quarter
Bio-10
最暖季均温Mean temperature of warmest quarter
Bio-14
最干月降水量Precipitation of driest month
Bio-16
最湿季降水量Precipitation of wettest quarter
Bio-17
最干季降水量Precipitation of driest quarter
Bio-19
最冷季降水量Precipitation of coldest quarter
空间约束因子Spatial constraint
X、X2、X2Y、Y、Y2、XY2、X3、Y3、XY
表选项
趋势面分析(Trend surface analysis, TSA)是利用数学曲面模拟地理系统要素在空间上的分布及变化趋势的一种方法[31], 实质是通过回归分析原理, 运用最小二乘法拟合一个二维非线性函数, 模拟变量在空间上的分布规律, 从而找出研究区域内变量的空间分布格局的一种数学方法[32]。模拟某一属性要素在空间上的分布规律[33], 并展示该要素在地域空间上的变化趋势, 利用趋势面拟合函数所得的结果绘制出的等值图叫做趋势图[34]。其计算:
在空间分析中, 趋势面的函数方程主要有一次趋势面、二次趋势面、三次趋势面、高阶趋势面模型函数等[34]。本文参照Tovaranonte将经度和纬度作为空间约束因子, 运用ArcGIS 10.3将长苞铁杉现有主要分布区和周边区域划分为1 km×1 km的栅格[18], 提取网格中心点经纬度坐标, 进行三项式空间趋势面计算, 得出X、X2、X2Y、Y、Y2、XY2、X3、Y3、XY共9项[8], 重新赋值于每个栅格后, 得到研究区趋势面图层, 代表空间约束因子(表 1)。
1.2 模型运算与结果可视化
本文采用的物种分布模型为最大熵模型(Maximum Entropy, MaxEnt), 它能够在物种分布记录较少的情况下使用, 且其精度和稳定性常优于其他模型[6, 8]。在进行模型运算前, 将环境数据通过ArcGIS 10.3转换为ASCII, 同时将长苞铁杉分布点数据保存为CSV格式。分别将不同环境模式数据带入MaxEnt模型运算, 设置75%的分布点作为模型训练集, 25%作为检验集, 选择jackknife计算各气候因子贡献率, 重复运算30次, 最大迭代次数为5000次。通过计算受试者操作特性曲线(receiver operating characteristic, ROC)的下方面积(area under curve, AUC)反应模拟结果的准确性。AUC评估标准为: 0.5—0.6预测无效；0.6—0.7较差；0.7—0.8一般；0.8—0.9良好；0.9—1极准确。其中AUC＞0.85的结果就可采纳, AUC＞0.9表示模拟结果非常准确[35-36]。
在划分适生区时, 首先根据前人对濒危、珍稀物种的适生区划分标准[10]和长苞铁杉本身分布特点, 设定其阈值最低为0.1, 并在0.1—1之间划分为3份, 即0.1—0.3、0.3—0.6和0.6—1, 分别对应低、中、高适生区, 并计算其面积。然后根据不同模式的稳定适生区, 在ARCGIS 10.3中计算随时间变化存在的较为稳定的潜在适生区作为长苞铁杉应对未来气候变化的避难所。最后对结果进行可视化表达并比较分析两种模型(气候因子预测模型(C)和气候+空间约束因子预测模型(C+S))的结果差异。
2 结果分析
2.1 模型准确性分析
本文通过模型运算得出ROC曲线的AUC值来验证MaxEnt模型对长苞铁杉的地理分布与气候关系的模型适用性。根据9个气候因子构建的模型AUC值为0.911, 表明该模型预测精度达到了“极准确”水平, 可以较好地对该物种在气候变化背景下的潜在分布区进行预测；而在C+S模型中AUC值达到0.934, 这表明加入空间约束因子的模型模拟结果优于单纯的气候因子模型(表 2)。
表 2 模型AUC值
Table 2 The value of AUC
模型Model
AUC值The value of AUC
C模型Model C
0.911
C+S模型Model C+S
0.934
表选项
2.2 影响长苞铁杉分布的主导变量
因在本研究中C+S模型预测效果优于C模型, 故本文仅分析C+S模型中的各变量因子对长苞铁杉地理分布的贡献率。贡献率排前五位的因子分别为最干月降水量、Y、XY、X、温度年较差。S(空间约束因子)对其分布的总贡献率达到51%(表 3), 而该模型下气候因子中最干月降水量对其分布贡献率最大, 达到19.9%(表 3)。
表 3 气候+空间约束因子预测模型中前五位环境变量贡献率
Table 3 The contribution rate of the first five environmental variables in the prediction results using model of climate + spatial constraint factor model (C+S)
模型Model
C+S模型Model C+S
变量Variable
14
y
xy
x
7
贡献率Percent contribution
19.9
16.5
15.4
10.3
8.4
表选项
2.3 适生区等级划分
两种模型预测的长苞铁杉现代适生区范围较为一致, 分布中心均在贵州、湖南、广西等省区交界处, 并向东延伸至江苏、福建等省区, 其预测范围与现有分布点位所在区域较吻合。在两种预测模型下, 长苞铁杉在未来不同气候变化模式下的各等级适生区均低于其现代适生区(表 4), 且C+S模型低于C模型。此外, C+S模型预测的长苞铁杉适生区分布更为集中, 连续性较高；而C模型中, 部分适生区呈零星和间断分布, 地域跨幅更广。
表 4 不同气候模式下长苞铁杉潜在适生区面积情况/×103 km2
Table 4 The area of the potential adaption zone of Tsuga longibracteata in different climate models
模型Model
适生区等级Grade of suitable area
1990
RCP2.6
RCP8.5
2050
2070
2050
2070
C模型Model C
低适生区
307.02
222.35
288.73
192.04
133.24
中适生区
203.69
78.03
107.49
76.04
42.43
高适生区
53.52
15.13
9.49
5.44
2.82
总适生区
564.23
315.51
405.71
273.53
178.49
C+S模型Model C+S
低适生区
290.04
226.24
255.59
193.35
168.27
中适生区
137.3
102.97
144.09
80.41
19.99
高适生区
41.73
23.2
48.24
8.5
0.14
总适生区
469.07
352.41
447.92
282.26
188.41
表选项
两种模型预测结果均表现出长苞铁杉潜在适生区随时间推移而面积减小的趋势, 尤其是中高适生区面积较现代减少幅度较大。在C模型中两模式均表现为湖南、广西等地的适生区范围缩小, 但RCP8.5模式适生区面积变化和北移趋势均更显著, 且重庆、湖北交界处以及浙江沿海地区的适生区面积扩大(图 2)。在C+S模式中, 长苞铁杉中高等级适生区呈现出随时间变化而面积降低的趋势(表 4), 但其分布中心始终位于湖南西南部、贵州东南部、广西北部、广东西北部的交界处。RCP8.5时适生区有大幅度缩减, 部分省份生境丧失和破碎化, 且北移趋势显著(图 3)。
图 2 C模型预测长苞铁杉的潜在适生区分布
Fig. 2 The prediction of potential distribution of Tsuga longibracteata suitable area by model of climate
图选项
图 3 C+S模型预测长苞铁杉的潜在适生区分布
Fig. 3 The prediction of potential distribution of Tsuga longibracteata suitable area by model of climate + spatial constraint
图选项
2.4 避难所预测
避难所是指在物种面临恶劣环境条件仍能存活的区域[37-38], 本文将未来不同时期和不同RCP模式下的长苞铁杉稳定适生区作为其应对未来气候变化的避难所[39]。就C模型来看, 长苞铁杉稳定的适生区以重庆东南部、湖北西南部、湖南西北部、贵州北部、广西北部等交界处为分布中心, 台湾山脉和浙江沿海地带也有部分分布, 尤其是湘、渝、鄂边境地区(大巴山、巫山)面积最为集中。而C+S模型下, 其分布中心位于湘、黔、桂地区交界处(包括云贵高原东缘、武陵山和南岭部分地区), 与长苞铁杉现有的分布格局更加接近。与C模型相比, 总体地理位置更靠南, 面积较大且适生区连续性较好(图 4)。
图 4 不同模型条件下的稳定适生区
Fig. 4 Stable zones of Tsuga longibracteata under different model
图选项
3 讨论
3.1 影响长苞铁杉地理分布的主导气候因子
气候是制约生物的生长发育、物种分布以及生物多样性丰富程度的最主要自然因素之一[40]。水热条件共同决定植物地理分布格局[41], 但二者的作用针对不同的物种也有差异[42]。在C+S模型中, 最干月降水量是影响长苞铁杉地理分布的主导气候因子。在应用物种分布模型对孑遗植物裸果木、新疆野生果树、蒙古栎和红松等地理分布格局的研究中, 均发现水分对其种群的分布格局影响更大[43-46]。与之类似, Vedelsørense对棕榈树的研究结果也表明最干季降水量是影响其分布的主导因子[17]。这些已有报道与本文的最干月降水量作为影响物种分布的主导因子的结果一致。这个结论与铁杉喜湿的生物学特性有关, 并在一定程度上得到了控制实验的验证。如, 朱小龙等通过对长苞铁杉幼苗存活率的研究发现土壤含水量对其生长有显著影响, 最干月降水量变化在一定程度上会影响土壤最低含水量变化, 从而引起干旱胁迫。长苞铁杉幼苗叶片叶绿素含量和比例、渗透调节物质、抗氧化酶和叶片丙二醛含量等生理指标受此影响, 从而影响长苞铁杉生长[47]。因此, 本研究推测在最干月降水量比现代显著降低的区域不利于长苞铁杉的地理分布, 适度增加则有利于其生长和分布。由此, C模型中2070年最干月降水量较2050年有了较大的变化, 使得其南部部分适生区(广西、广东等地)消失；这也是C+S模型中2070年, RCP2.6时, 最干月降水量的适量增加使长苞铁杉在北部的渝、鄂、湘交界处的适生区面积略有扩大的原因。
3.2 长苞铁杉未来气候变化下分布格局变迁
气候变化对物种适宜生境面积的影响较为复杂, 根据物种自身对逆境抵抗力的不同, 呈现出扩张、收缩以及稳定不变的趋势[8, 10, 48-49]。在本研究中, 两种模型预测下的长苞铁杉潜在适生区均表现出面积减小, 且有不同程度的位置北移。在先前开展的利用物种分布模型模拟未来气候变化下喜马拉雅地区的百合科药用植物、落叶松和清香木的潜在分布格局的研究中, 也有近似的发现[12, 42, 50]。此外, 未来气候变化下长苞铁杉潜在适生区缩小也符合区域尺度上长苞铁杉地理分布变化的野外调查结果。邱迎君等人通过长期定点观察, 指出气候变化能导致长苞铁杉地理分布格局改变, 主要表现为居群范围缩小和分布位置向北和海拔更高的地区迁移[51], 从而导致长苞铁杉适生区范围锐减和生境破碎化[52]。未来气候变化将对长苞铁杉的适宜生境造成较大负面影响, 可能导致其受威胁程度加剧, 并危及未来该物种存续。
加入空间约束变量参与模型运算, 预测出来的潜在适生区范围缩减和北移的程度与未加入的存在明显差异。C模型预测的潜在分布区更分散, 北移趋势更突出(图 3)；而加入空间约束因素后, C+S模型的预测结果则表现为适生区集中, 北移趋势较小。物种的分布范围随外界环境改变逐渐变化, 物种受空间约束的限制, 不能在较短的时间内迁移到距离现代分布区域较远的地方, 因此其适生区分布应更加集中和连续[8]。本研究, C模型潜在分布区离散程度更高, 与之相反的是C+S模型下的连续性好, 这也在一定程度上印证了这一观点。同时, 该结果也意味着C模型的预测结果更趋近于长苞铁杉在理论条件下的最大气候适宜生境, 而C+S模型的预测结果可能更符合未来气候变化下其潜在适生区的实际变化状况。这一结论得到了相关研究的验证[8, 17]。
在C+S模型中, 空间约束变量因子对长苞铁杉未来分布的影响贡献率达到51%, 这表明空间约束作用在一定程度上对其分布起到主导作用。事实上, 气候、历史分布、地形、土壤、人类干扰和空间约束等因素对物种在不同空间尺度上的分布均有重要影响[16]。在大尺度上, 气候因子是物种潜在分布变化的主要驱动力[8, 53]；而空间约束等往往能在区域尺度上对物种分布产生不可忽视的影响。已有研究表明, 空间约束主要与物种迁移扩散相关, 物种扩散的时间长短、物种本身扩散能力和地理阻隔等因素都与空间约束相联系[8, 16]。植物迁移速率和适生区变化速率的差值往往是导致物种丧失原有生境的主要原因之一[54]。对长苞铁杉而言, 已有研究发现其种子萌发率较低, 且扩散距离小于50 m, 扩散速率较慢, 使得长苞铁杉的迁移能力较差, 不能快速进入和占据理论上的气候适宜生境[55-56]。长苞铁杉适生区范围内山脉纵横交错(雪峰山、武陵山、大娄山、南岭等), 河网密布(乌江、柳江、沅江等), 地理阻隔作用显著, 这也增加了其向北迁移扩散的难度, 进一步阻碍了其迁移到气候理想的适生区。这为本研究得出长苞铁杉潜在分布受到空间约束变量的明显制约, 气候变化下其适生区面积减少、分布区集中的结论提供了解释。
3.3 避难所与未来引种地研究
避难所在物种应对频繁的气候变化, 免遭劫难方面具有重要作用[5, 57]。前人研究指出现代分布格局与未来各个时间段的潜在适生区中的重叠部分, 是物种稳定分布的区域, 可以作为物种应对气候变化的避难所[39, 43]。C+S模型下稳定适生区位于湘、黔、桂交界处(包括云贵高原东缘山地、武陵山和南岭部分地区)以及湘、川、鄂结合部, 这与长苞铁杉现有的分布区较为一致。而在C模型中, 稳定适生区分布区大致位于湘、川、鄂交界处, 集中在大巴山、巫山等地, 范围更加靠北。高大山脉多作为植物应对气候变化的避难所的作用已被广大研究者所认识。已有研究发现在气候陡变时期, 山区由于地形复杂多样, 能保留或形成许多独特的稳定小气候区域, 有利于物种保存[38]。较高海拔落差的山地为物种迁移提供了独特的垂直生境带, 可以让物种随温度变化而上下位移, 一定程度上避免了高山、河流等导致的地理阻隔作用对物种水平迁移的限制[38, 58]。已有研究表明, 更新世末期的冰期—间冰期旋回过程中, 中国大陆经历了多次降温-升温气候变化过程, 造成了大量物种的灭绝[59]。在间冰期中, 温度的上升不利于耐寒的针叶物种的生存, 导致其分布范围大幅度缩小, 而山地系统尤其是西南横断山等高大山脉具备热带、亚热带至高山寒带各类生境带, 缓解了物种受温度变化的不利影响, 是针叶植物有效的避难所[39, 41, 60-61]。从地形条件上来看, 无论是长苞铁杉现有生境的武陵山、南岭, 还是C模型预测下的避难所大巴山、巫山等山区, 其山体高度大多超过了2000米, 在气温上升时能提供相对冷寒的生境, 因此推测这些区域可能成为长苞铁杉的避难所。
C+S模型考虑到了空间约束对物种分布的制约作用, 其预测结果相对更加符合实际。对避难所而言, C+S模型下的稳定适生区与原有长苞铁杉的分布区位置大致相同, 这表明其可能具有“原地避难”的特点, 湘、黔、桂结合部的山地应是未来的重点保育区, 应在这些区域内进一步加强保护力度, 通过建立自然保护区, 禁止滥砍滥伐, 优化工程布局等措施, 切实保护现有生境, 为气候变化下长苞铁杉的物种迁移创造良好的条件。C模型预测的大巴山、巫山等地区, 距离长苞铁杉现有核心分布区较远, 更倾向于是长苞铁杉的理论气候适宜区, 可将其作为长苞铁杉未来引种、迁地保存的区域。事实上, 对濒危物种而言, 就地保护和迁地保护同样重要[62]。建议在这些区域内, 预先开展引种实验和研究, 做好引种区规划, 为长苞铁杉今后的引种和异地保存奠定基础。
4 结论与展望
最干月降水量是影响长苞铁杉地理分布的主导气候因子, 空间限制对其分布也有重要作用。未来气候变化下, 长苞铁杉各等级适生区面积随时间推移呈下降趋势, 分布范围有不同程度的北移, 受胁程度进一步加剧。C+S模型的预测显示长苞铁杉未来的核心分布区仍位于现存的湘、黔、桂结合部, 更加符合其迁移、扩散特性, 应进一步加强对这些区域的保护工作。C模型预测渝、川、鄂结合部的大巴山等地区是未来气候变化下长苞铁杉的理论分布区域, 可作为长苞铁杉应对未来气候变化的引种地区, 应提早进行相关研究工作。
本研究在模型构建时加入了空间约束因子, 提高了模型适应性和精度。今后，应尝试加入更多的影响因素来构建模型, 以提高模型预测的精确度, 为研究气候变化条件下濒危、珍稀植物的物种分布区变化和物种保护提供科学依据。