气候是影响物种分布的重要驱动因素^[1]。联合国政府间气候变化专门委员会(IPCC)第六次评估报告指出, 在当前的温室气体排放情景下, 到2100年全球温升可达到2.2—3.5℃^[2]。由于全球气候变化的区域变异性以及物种本身的生态位特征, 不同的物种会表现出不同的响应模式。目前, 已有许多物种表现出向高纬度、高海拔地区扩散或转移的现象^[3—5]。但某些物种也会呈现出不同的变化趋势, 如我国新疆地区的鸟类向东扩张^[6], 以及部分物种的分布区萎缩^{[5, 7—8]}。此外, 由于宠物贸易、园林建设等人为引入的原因, 也导致了外来物种的定殖和扩张^[9—10]。

灰喜鹊(Cyanopica cyanus)隶属于雀形目(PASSERIFORMES)鸦科(Corvidae)灰喜鹊属, 主要分布于东亚, 其中我国北方是其主要分布区^[11]。主要以鳞翅目(LEPIDOPTERA)幼虫和直翅目(ORTHOPTERA)昆虫为食, 是重要的食虫鸟类, 在国内常被人为引种用以防治松毛虫(Dendrolimus)等森林害虫^{[9, 12—14]}, 对保持森林生态系统的稳定性具有重要作用。由于人为引种、笼养逃逸及放生等原因, 已有一定数量的灰喜鹊在其自然分布地以外的地方建立种群。1975年, 灰喜鹊作为逃逸物种在香港被记录^[15]; 广东广州、广西桂林和福建光泽分别在1985年、1986年和1992年引进灰喜鹊以防治森林虫害^{[12—13, 16]}; 1988年, 昆明动物园从北京引进灰喜鹊以供游客观赏, 后于1991年释放野外^[17]; 2011年, 王龙舞等人在海南海口发现本地笼养逃逸的灰喜鹊种群^[15]。目前, 灰喜鹊已在甘肃^[18]、云南^[17]、广西^[16]、广东^[13]、福建^[12]、海南^[15]、香港^[15]及台湾^[19]等省份定殖并呈现出扩散的趋势。

物种分布模型(Species Distribution Model, SDM)也常称为生态位模型(Ecological Niche Model, ENM), 是一种结合了物种分布数据和相关环境变量数据的数学模型^[20—21], 此类模型基于一定的算法, 在给定的环境变量构成的生态空间内, 根据分布数据及其所关联的环境变量参数来计算物种对生态位的需求, 并将此结果投射到给定的时间和空间范围内, 最后以概率的方式呈现物种对不同时空中生境的喜好程度^[21—23]。近年来, 物种分布模型被广泛地应用于保护生物学、入侵生物学及气候变化对物种分布影响的研究中^{[21—22, 24]}。最大熵(Maximum Entropy, MaxEnt)模型是目前应用最为广泛、性能可靠的物种分布模型之一^[24—27]。

本研究基于气候生态位理论, 结合灰喜鹊在我国境内的分布数据及气候变量, 通过构建MaxEnt模型, 旨在探讨以下科学问题：(1)灰喜鹊自然分布地种群和引入地种群对气候因子的响应情况; (2)在人为引入与气候变化背景下, 当前和未来时期灰喜鹊种群在我国境内的分布格局与变化趋势。这将有助于理解全球气候变化背景下物种的分布及变化规律, 并预测灰喜鹊对生态系统的潜在影响。

1 数据来源与研究方法 1.1 分布数据获取及处理

分布数据的来源包括全球生物多样性信息数据库(https://www.gbif.org/)、中国观鸟记录中心(http://www.birdreport.cn/)以及文献记录^{[12—13, 15]}。共获取到61026条分布数据, 分布数据收集时间截止2022年7月。

首先删除信息错误、缺失、模糊及重复的分布数据, 然后将剩余数据在ArcMap 10.6中与使用的环境变量叠加, 删除在环境变量范围外的分布数据。为了降低空间自相关对模型构建的干扰作用, 使用R 4.1.3中的Wallace包对分布数据进行筛选, 每个2.5 min×2.5 min的栅格里仅保留1个分布点^[28—29]。最终保留有效分布点2458个。

以《中国动物图谱》^[30]中灰喜鹊分布界线为参照, 将位于该区域内的分布点定义为自然分布点, 共计2203条。将位于自然分布地外且有明确的证据表明灰喜鹊为引入物种的地点定义为引入地, 最终确定新疆乌鲁木齐、云南昆明、贵州贵阳、广西桂林、广东广州、福建光泽、海南、香港和台湾等9个省市为引入地^{[12—13, 15—17, 19]}, 筛选出有效分布点129个(图 1)。

1.2 环境变量获取及筛选 1.2.1 环境变量获取

本研究共选取世界气候数据库(https://www.worldclim.org/)中的19种气候变量用于模型分析(表 1), 分辨率为2.5 minutes。未来时期的环境变量选择2041—2060年(2050s时期)和2061—2080年(2070s时期)两个阶段, 气候模式选择北京气候中心发布的BCC-CSM2-MR气候系统模式, 该模式对于中国的气候变化趋势具有较好的模拟效果^[31], 情景模式选择ssp245, 即中等辐射强迫与中等社会脆弱性的组合情景^[32]。

1.2.2 环境变量筛选

环境变量在物种分布模型中的过度拟合问题受到了广泛的关注, 为了避免模型的过度拟合, 剔除部分环境变量可以使模型的预测性能得到一定的提升^[33—34]。学者们普遍采用Pearson相关系数分析法对模型使用的环境变量进行多重线性相关性分析, 以环境变量间的相关性系数结合环境变量的贡献率作为筛选的依据^{[22, 35]}。使用ArcMap 10.6中的“Band Collection Statistics”工具对使用的19种环境变量进行相关性分析, 相关性系数r表示两两变量之间线性关系的密切程度, 本研究选择| r |>0.8作为环境变量间有强相关性的标准^[35—36]; 然后将19种环境变量导入到MaxEnt模型中试运行, 使用模型结果中环境变量的贡献值作为筛选的依据^{[27, 33, 37]}, 环境变量的贡献率表示环境变量的重要程度, 贡献率越高的环境变量对预测结果的贡献越大^[38]; 最后结合以上两种方法的结果来筛选环境变量, 当两个环境变量之间的相关性系数r>0.8时, 则保留贡献率更高的环境变量。最终优选出bio3、bio4、bio10、bio11、bio12和bio15等6种环境变量用于综合分布模型与自然分布地模型, bio2、bio3、bio6、bio7、bio8、bio13、bio14和bio15等8种环境变量用于引入地模型。

1.3 MaxEnt模型构建 1.3.1 模型参数优化

使用R 4.1.3中的Wallace包分别对自然分布地模型、引入地模型及综合分布模型进行计算, 优选出每个模型最佳的调控倍频(Regularization multiplier, RM)和特征组合(Feature combination, FC)。在Wallace中设置调控倍频的范围为0.5—4.0, 间隔为0.5;特征组合项选择：H、L、LQ、LQH、LQHP(H, hinge, 片段化; L, linear, 线性; Q, quadratic, 二次型; P, product, 乘积型)^{[29, 39—41]}。delta.AICc是一种模型评估指标, 在候选模型中delta.AICc值最低的模型被确定为最优模型^[42], 根据Wallace计算结果中最低的delta.AICc值选择3个模型最优的RM和FC：自然分布地模型使用RM=0.5, FC=H; 引入地模型使用RM=3, FC=LQHP; 综合分布模型使用RM=0.5, FC=LQHP。

MaxEnt模型的其他参数设置如下：选择25%的分布数据作为测试集, 75%的分布数据作为训练集; 设置最大背景点数为10000;模型重复次数为10次; 重复运行模式选择Crossvalidate, 即交叉验证法; 最大迭代次数为1000次; 选择10 percentile training presence阈值规则作为适生区与非适生区的划分依据; 其余模型参数为默认参数。

1.3.2 适生区划分

采用10 percentile training presence阈值作为划分适生区与非适生区的界限, 适宜值低于该阈值的地方为非适生区, 高于或等于该阈值的地方为适生区^{[35, 43]}, 其中自然分布地模型的阈值为0.2629, 引入地模型的阈值为0.1170, 综合分布模型的阈值为0.2579。然后将适生区按照适宜值由低到高三等分为低适生区、中适生区和高适生区^[38]。

1.3.3 模型预测精度评估

以AUC值作为模型预测精度评估的标准, 值越大表示模型预测结果越有规律, 环境变量与物种地理分布之间的相关性越大, 模型预测效果越好^{[26, 38, 44]}。当AUC值为0.5—0.6时, 代表模型预测失败; 当AUC值为0.6—0.7时, 代表模型预测效果较差; 当AUC值为0.7—0.8时, 代表模型预测效果一般; 当AUC值为0.8—0.9时, 代表模型预测效果良好; 当AUC值为0.9—1.0时, 代表模型预测效果优秀^[45—46]。

2 结果与分析 2.1 模型评估结果

自然分布地模型、引入地模型和综合分布模型的AUC值分别为0.851、0.962和0.826, 为良好或优秀水平, 表明模型预测结果较好。

2.2 环境变量对灰喜鹊分布的影响

对自然分布地种群分布影响较大的环境变量为bio12(年均降雨量)、bio4(温度季节性变化)和bio11(最冷季度均温)(图 2)。自然分布地种群较适宜(存在概率>0.5)的年均降雨量的范围为507.45—760.84 mm和868.06—1345.63 mm, 较适宜的温度季节性变化标准差的范围为853.59—1138.49, 较适宜的最冷季度均温的范围为-3.80—6.29℃, 自然分布地种群对于3种环境变量的响应曲线均呈现出先增加后减少的“单峰”曲线(图 3)。

对引入地种群分布影响较大的环境变量为bio6(最冷月份最低温)和bio14(最干月降雨量)(图 2)。引入地种群较适宜的最冷月份最低温的范围为7.99—17.30℃, 较适宜的最干月降雨量的范围为16.65—202.32 mm, 引入地种群的存在概率均随着最冷月份最低温和最干月降雨量的增加而升高(图 4)。

对灰喜鹊种群分布影响较大的环境变量为bio12(年均降雨量)、bio11(最冷季度均温)和bio4(温度季节性变化)(图 2), 表明灰喜鹊更偏好温暖潮湿的东部季风区。灰喜鹊种群较适宜的年均降雨量的范围为509.30—1347.55 mm和2310.55—4176.65 mm, 较适宜的最冷季度均温的范围为-4.00—6.41℃和15.34—19.62℃, 较适宜的温度季节性变化标准差的范围为849.26—1142.89, 由于综合分布模型的分布数据主要来自自然分布地, 所以贡献率较高的环境变量与自然分布地模型一致, 但也由于引入地分布数据的影响, 灰喜鹊种群对这3种环境变量的响应曲线均呈现出了先增加后减少、再增加再减少的“双峰”曲线(图 5)。

2.3 潜在适生区的分布与变化趋势 2.3.1 自然分布地模型

通过对比当前和未来时期的适生区面积和分布发现：自然分布地种群适生区的面积和适宜值都将随着气候的变化而大幅增加, 扩张的区域主要位于高纬度的东北地区和高海拔的青藏高原(图 6, 图 7)。当前时期的适生区面积共175.46万km², 占中国国土总面积的18.24%, 主要分布于自然分布地内的华北、华中和华东地区, 其中多数为中适生区, 高适生区主要分布于湖北、安徽、江苏及渤海周围的陆地区域, 自然分布地外, 仅青藏高原存在较多的适生区, 西藏的札达县还出现了小面积的高适生区; 2050s时期的适生区面积共252.24万km², 占中国国土总面积的26.23%, 相比当前时期增加了43.76%, 新增的适生区主要位于东北地区和青藏高原。其中高适生区在河北、山东、河南、安徽、江苏和湖北等地大幅增加, 表明灰喜鹊不仅有明显扩张的现象, 自然分布地种群的密度可能也会随适宜值的升高而增加; 2070s时期的适生区面积共283.66万km², 占中国国土总面积的29.49%, 相比2050s时期增加了12.46%, 适生区面积小幅增加, 分布上的变化较小, 仅新增了湖南和川东地区, 表现出了轻微的南扩现象。质心分析表明灰喜鹊适生区总体向高纬度方向转移, 相比当前时期, 2050s时期的质心朝北偏东23°方向转移了253.38 km, 2070s时期的质心朝北偏东25°方向转移了229.16 km(图 7)。

2.3.2 引入地模型

通过对比当前和未来时期的适生区面积和分布发现：引入地种群适生区的面积和适宜值均随着气候的变化而小幅增加, 分布上的变化较小(图 6, 图 8)。当前时期的适生区面积共103.34万km², 占中国国土总面积的10.74%, 主要分布于自然分布地外的新疆北部、四川东部、云南、广西、广东、福建东南部、海南和台湾, 此外在自然分布地内的重庆、湖北、江苏和上海也存在少量的适生区, 其中中、高适生区主要分布于广西、广东、海南和台湾等地, 这表明引入地种群很可能在这些沿海地区扩张和定殖; 2050s时期的适生区面积共125.82万km², 占中国国土总面积的13.08%, 相比当前时期增加了21.75%;2070s时期灰喜鹊的适生区面积共128.98万km², 占中国国土总面积的13.41%, 相比2050s时期增加了2.5%。质心分析表明灰喜鹊适生区总体向高纬度方向转移, 相比当前时期, 2050s时期的质心朝北偏东37°方向转移了88.85 km, 2070s时期的质心朝北偏东46°方向转移了67.69 km(图 8)。

2.3.3 综合分布模型

通过对比当前和未来时期的适生区面积和分布发现：灰喜鹊种群适生区的面积和适宜值都将随着气候的变化而大量增加, 扩张的区域主要位于高纬度的东北地区和新疆北部、高海拔的青藏高原以及南部的四川东部和两广地区(图 6, 图 9)。当前时期的适生区面积共226.82万km², 占中国国土总面积的23.58%, 主要分布于自然分布地内的华北、华中和华东地区, 与自然分布地模型类似, 自然分布地内的低、中、高适生区的分布格局也与自然分布地模型基本一致, 此外在自然分布地外的新疆北部、青藏高原、广西南部、广东南部、海南和台湾也存在一定面积的适生区, 但适宜程度普遍较低; 2050s时期的适生区面积共378.97万km², 占中国国土总面积的39.4%, 相比当前时期增加了67.08%, 适生区面积大幅增加, 新增的适生区主要位于东北地区、新疆北部、青藏高原、四川东部和两广地区, 新增的适生区中主要为中、低适生区, 当前时期的中适生区大部分转变为高适生区; 2070s时期的适生区面积共431.53万km², 占中国国土总面积的44.87%, 相比2050s时期增加了13.87%, 适生区面积和适宜值进一步增加, 但分布上的变化较小。质心分析表明灰喜鹊适生区总体向高纬度方向转移, 相比当前时期, 2050s时期的质心朝北偏西19°方向转移了164.07 km, 2070s时期的质心朝北偏西34°方向转移了153.44 km(图 9)。

3 讨论 3.1 影响灰喜鹊分布的主要环境变量

气候因子对物种分布的影响与物种的适应性以及食物的可利用性有关^{[8, 47—48]}。在自然分布地模型与综合分布模型中, 贡献率最高的环境变量均为年均降雨量。降水主要通过影响植被和食物, 从而间接影响鸟类的分布^[48—50]。在降水量充沛的地区, 鸟类的栖息地和昆虫等食物资源也相对充足^[50]。此外松毛虫作为灰喜鹊的重要食物来源^[12], 松毛虫的卵、幼虫和蛹需要较湿的环境来保证其正常的生长发育, 在干旱的环境中松毛虫的密度会明显下降^[51]。所以在年均降雨量的响应曲线中, 灰喜鹊的存在概率在较低的年均降雨量范围内会随着年均降雨量的升高而迅速增加。有研究表明降雨量与鸟类丰度呈正相关^{[50, 52]}, 本文结果与之一致。在引入地模型中, 年均降雨量因贡献率较低而被最干月降雨量替代, 其原因可能是引入地大多为降雨量充足的南方地区, 在降雨量有保障的前提下, 年均降雨量就不再是影响灰喜鹊分布的主要环境变量, 而最干月降雨量则依然是影响引入地种群分布的主要环境变量。除年均降雨量外, 最冷季度均温(或最冷月份最低温)在3种模型中也表现出了较高的贡献率。最冷季度均温过低会影响鸟类的身体机能^[47], 此外在最冷季度均温较高的地区, 鸟类可利用的食物资源也相对较充足^{[47, 50]}。所以在环境变量响应曲线中, 最冷季度均温(或最冷月份最低温)较高的地方, 灰喜鹊的存在概率也较高。

比较3种模型的环境变量响应曲线的差异, 发现引入地种群适宜的最冷季度均温(或最冷月份最低温)和年均降雨量比自然分布地种群更高, 这和引入地多为温暖湿润的南方地区有关。表明灰喜鹊从自然分布地引入到引入地后, 气候生态位产生了重要的变化, 证明灰喜鹊具有较强的适应能力, 这在胡箭关于灰喜鹊在昆明城区的生态适应的研究中也得以验证^[17]。

3.2 灰喜鹊的分布格局与变化趋势

本研究结果表明, 当前时期灰喜鹊在我国境内的适生区主要分布于其自然分布地内的华北、华中和华东地区。在未来气候变化的背景下, 3种模型的结果均显示灰喜鹊有显著的扩张趋势, 主要表现为向高纬度的东北地区和高海拔的青藏高原地区扩散。许多国内外研究也表明, 在全球气候变化条件下, 物种有向高纬度、高海拔地区迁移的趋势^[8], 该现象在鸟类中也普遍存在^[5]; Hitch和Leberg研究了北美地区56种鸟类的扩散趋势, 发现分布在低纬度地区的鸟类平均每年向北移动2.35 km^[53]; Parmesan在对鸟类等物种的研究中发现, 这些物种都表现出了向高海拔地区转移的现象^[54]; 国内的一些研究也证明了我国的鸟类有同样的扩散趋势：杜寅等分析了我国鸟类20年间的分布数据, 发现共有120种鸟类的分布范围已向北扩张^[48]; 原分布于印度及东南亚地区的钳嘴鹳(Anastomus oscitans), 近年来已向北扩散至我国^[55]; 李东明等在青海西宁发现了约60只白头鹎(Pycnonotus sinensis), 表明白头鹎不仅存在向北扩散的现象, 还存在向高海拔地区扩散的趋势^[56]; 鲜莉莉在青藏高原东北部对灰喜鹊开展的繁殖研究中表示, 灰喜鹊是近年来通过殖民的方式扩散至青藏高原东北部的^[18], 这也从另一方面验证了本研究结果的正确性。此外, 部分鸟类也会呈现出不同的变化趋势：马鸣研究发现由于人为因素、生境的可入侵性、全球气候变化及生物地理上的空缺等原因, 在我国新疆有至少42种鸟类存在向东扩散的现象^[6]; 多数鸟类的分布范围在逐渐缩小, 特别是分布在高纬度、高海拔地区的鸟类, 几乎没有可供其拓展的生存空间^{[5, 8]}。

我国气候的变化趋势与全球气候变化总体趋势基本一致^[48]。我国大部分地区的气温将呈现出增长的趋势, 且高纬度地区的气温增幅较大^{[8, 48—49]}。我国的降水量总体上同样处于增加的趋势^[8], 有研究发现青藏高原的降水量正在变得越来越多^[57—58], 北半球中高纬度的降雨量也处于增长的趋势^[5]。本研究结果发现灰喜鹊更偏好温暖潮湿的气候, 我国气候总体上呈暖湿化趋势, 且高纬度、高海拔地区的温度和降水增加较快, 这恰好与本研究结果中灰喜鹊种群的扩散方向一致。所以本研究认为, 我国气候的暖湿化趋势将导致灰喜鹊向高纬度和高海拔地区扩散以及其种群数量的增加。

由于人为引种、笼养逃逸及放生等原因, 目前灰喜鹊已在其自然分布地外建立了多处可自我维持的种群, 在自然分布地外的西南地区和华南地区也存在较大面积、较高适宜值的适生区。灰喜鹊之所以能在其自然分布地外的南方地区定殖, 很大的原因在于灰喜鹊偏好暖湿的气候, 而南方地区相比其自然分布地更加温暖潮湿。其次灰喜鹊较强的适应能力^[17], 也使其能在气候条件完全不同的地方存活, 引入地种群的气候生态位可能已产生了重要变化。虽然引入灰喜鹊对控制松毛虫等森林害虫具有积极的意义, 但其较强的适应力^[17]和竞争力^[18]也可能会对部分食虫鸟类, 尤其是同样以松毛虫为食的松鸦(Garrulus glandarius)、喜鹊(Pica pica)、大杜鹃(Cuculus canorus)及四声杜鹃(Cuculus micropterus)等鸟类造成影响^{[11, 51]}。