生态学报  2022, Vol. 42 Issue (7): 2636-2647

文章信息

赵梓伊, 肖能文, 刘高慧, 李俊生
ZHAO Ziyi, XIAO Nengwen, LIU Gaohui, LI Junsheng
五种齿突蟾在横断山南潜在地理分布预测
Prediction of the potential geographical distribution of five species of Scutiger in the south of Hengduan Mountains Biodiversity Conservation Priority Zone
生态学报. 2022, 42(7): 2636-2647
Acta Ecologica Sinica. 2022, 42(7): 2636-2647
http://dx.doi.org/10.5846/stxb202103020572

文章历史

收稿日期: 2021-03-02
网络出版日期: 2021-12-15
五种齿突蟾在横断山南潜在地理分布预测
赵梓伊1,2 , 肖能文1 , 刘高慧1 , 李俊生1     
1. 中国环境科学研究院 国家环境保护区域生态过程与功能评估重点实验室, 北京 100012;
2. 兰州大学生命科学学院, 兰州 730000
摘要: 由于栖息地质量下降, 近年来齿突蟾属物种种群数量急剧减少, 明确齿突蟾属物种空间分布, 是监测、管理、保护齿突蟾属物种的基础。横断山区可能是齿突蟾属的起源中心和分化中心, 但齿突蟾属在横断山区的地理分布格局尚不明确。利用优化后Maxent模型, 首次预测西藏齿突蟾Scutiger boulengeri、刺胸齿突蟾Scutiger mammatus、胸腺齿突蟾Scutiger glandulatus、圆疣齿突蟾Scutiger tuberculatus、贡山齿突蟾Scutiger gongshanensis 5种高海拔齿突蟾属物种在横断山南生物多样性保护优先区域的潜在地理分布, 并分析其与环境因子的关系。结果显示, 5种齿突蟾属物种在横断山南的潜在地理分布格局存在差异, 西藏齿突蟾主要分布在横断山南的北部, 圆疣齿突蟾主要分布在横断山南东北部的四川省境内, 贡山齿突蟾主要分布在横断山南的西南部, 刺胸齿突蟾和胸腺齿突蟾的潜在分布格局较为相似, 在横断山南的中部、西北部地区都有较多分布, 但胸腺齿突蟾潜在分布区更为碎片化。另外, 横断山南北部地区的齿突蟾属丰富度明显高于南部地区。环境变量贡献率和刀切法结果显示温度因子和降水因子是决定横断山南齿突蟾属潜在分布的主要因素, 最冷季降水量对西藏齿突蟾、贡山齿突蟾、圆疣齿突蟾潜在分布有重要影响, 但它们对最冷季降水量的偏好存在差异。此外, 研究也显示, 通过评估潜在的Maxent参数组合, 选择最佳的Maxent模型是有效且必要的。
关键词: Maxent模型    潜在地理分布    横断山南    齿突蟾属    
Prediction of the potential geographical distribution of five species of Scutiger in the south of Hengduan Mountains Biodiversity Conservation Priority Zone
ZHAO Ziyi1,2 , XIAO Nengwen1 , LIU Gaohui1 , LI Junsheng1     
1. Key Laboratory of National Environmental Protection for Assessment of Region Ecological Processes and Functions, Chinese Research Academy of Environmental Sciences, Beijing 100012, China;
2. School of Life Sciences, Lanzhou University, Lanzhou 730000, China
Abstract: Due to the decreased habitat quality, the population size of Scutiger has declined dramatically in recent years. A clear understanding of the geographic distribution about Scutiger is the basis for monitoring, management, and conservation of Scutiger. The Hengduan Mountains region may be the origin and differentiation center of Scutiger. However, the geographical distribution pattern of Scutiger in the Hengduan Mountains region is still unclear so far. To explore this issue, we used the optimized Maxent model to predict the potential geographical distribution of five high-altitude species of Scutiger in this study, including Scutiger boulengeri, Scutiger mammatus, Scutiger glandulatus, Scutiger tuberculatus, and Scutiger gongshanensis, in the south of Hengduan Mountains Biodiversity Conservation Priority Zone. In addition, we further analyzed the relationship of potential geographical distribution with the environmental variables of five high-altitude species of Scutige. The results showed that the potential geographical distribution patterns of the five species of Scutiger were different in the south of the Hengduan Mountains Biodiversity Conservation Priority Zone. The potential distribution area of Scutiger boulengeri was mainly located in the northern part of the south of the Hengduan Mountains Biodiversity Conservation Priority Zone. Scutiger tuberculatus was mainly distributed among Sichuan Province in the northeast of the south of Hengduan Mountains Biodiversity Conservation Priority Zone. However, Scutiger gongshanensis was mainly distributed in the southwest of the south of Hengduan Mountains Biodiversity Conservation Priority Zone. The potential distribution pattern of Scutiger mammatus and Scutiger glandulatus was similar, mainly distributed in the central and northwest parts of the south of Hengduan Mountains Biodiversity Conservation Priority Zone, while the potential distribution of Scutiger glandulatus was more fragmented. What's more, the species richness of Scutiger in the northern part of the south of the Hengduan Mountains Biodiversity Conservation Priority Zone is significantly higher than that in the southern part. The relative contributions and jackknife test of environmental variables showed that temperature and precipitation were the main environmental factors, which determine the potential distribution of Scutiger in the south of the Hengduan Mountains Biodiversity Conservation Priority Zone. Precipitation of the coldest quarter had an important effect on the potential distribution of Scutiger boulengeri, Scutiger gongshanensis, and Scutiger tuberculatus, but their preferences for precipitation of the coldest quarter were different. In addition, our study also shows that it is effective and necessary to select the best Maxent model by evaluating the potential model combination of parameters for improving the Maxent model's predictive power.
Key Words: Maxent model    potential geographical distribution    the south of Hengduan Mountains Biodiversity Conservation Priority Zone    Scutiger    

由于栖息地丧失和破碎化、传染病以及全球气候变化的影响, 两栖动物的数量和种类正在快速的减少, 甚至部分物种可能灭绝[ 1- 4]。同时由于对两栖动物的地理分布缺乏了解, 两栖动物的保护行动受到阻碍, 因此, 明确两栖动物栖息地的空间分布格局是非常必要的。

物种分布模型可以有效地预测物种空间分布格局。近年来, 物种分布模型快速发展, 已被广泛应用于自然资源管理[ 5]、预测入侵物种和病虫害的潜在分布区[ 6]、预测珍稀濒危物种的适宜生境[ 7]以及评估全球气候变化对物种分布的影响[ 8- 10]等方面。最大熵模型(maximum entropy model, Maxent)是基于最大熵理论的物种分布模型, 它将已知物种分布点的环境变量特征作为约束条件, 寻找在此约束条件下熵最大的物种概率分布, 从而预测物种潜在分布范围[ 11]。与其他物种分布模型相比, Maxent模型只需要物种出现点数据, 在样本量较小的情况下, 也能得到准确的预测结果[ 11- 13]。近年来, Maxent模型已成为预测物种空间分布的首选[ 5, 9- 10, 14- 15], 在预测两栖物种潜在分布方面也有广泛应用。例如, Zank等利用Maxent模型预测新热带区24种黑昧蟾属蟾蜍(Melanophryniscus)的当前潜在分布区和未来潜在分布区[ 16]。黄勇杰等利用Maxent模型预测海南岛海南臭蛙(Odorrana hainanensis)的潜在分布以及影响海南臭蛙潜在分布的主要环境因子[ 17]。张凯龙等利用Maxent模型评估镇海林蛙(Rana zhenhaiensis)的潜在地理适宜性和影响其潜在分布的主要因素[ 18]

齿突蟾属(Scutiger)隶属两栖纲(Amphibia)无尾目(Anura)角蟾科(Megophryidae), 我国有已知物种18种[ 19], 主要分布在海拔3000-4500 m左右的高寒山区及高原地区[ 20]。齿突蟾属大多数物种为珍贵、稀有物种, 由于栖息地质量下降, 近年来种群数量急剧减少[ 21]

明确齿突蟾属物种空间分布, 是认识和保护齿突蟾属物种的基础, 齿突蟾属绝大多数种分布在我国西南山区, 其中横断山区可能是它们的起源中心和分化中心[ 20, 22], 但齿突蟾属在横断山区的具体分布情况仍不明确。

本研究以《中国生物多样性保护战略与行动计划》划定的横断山南段生物多样性保护优先区域(以下简称横断山南)为研究区, 横断山南分布有齿突蟾属物种8种[ 20, 23- 24], 包括西藏齿突蟾(Scutiger boulengeri)、刺胸齿突蟾(Scutiger mammatus)、胸腺齿突蟾(Scutiger glandulatus)、贡山齿突蟾(Scutiger gongshanensis)、九龙齿突蟾(Scutiger jiulongensis)、木里齿突蟾(Scutiger muliensis)、圆疣齿突蟾(Scutiger tuberculatus)、腾冲齿突蟾(Scutiger tengchongensis)。为了保证建模数据充足和模型预测的准确性, 本研究选取具有5个以上存在点数据的齿突蟾属物种, 分别是西藏齿突蟾、刺胸齿突蟾、贡山齿突蟾、圆疣齿突蟾以及胸腺齿突蟾, 通过优化Maxent模型首次预测五种齿突蟾属物种在横断山南的潜在地理分布, 分析影响五种齿突蟾属物种栖息地适宜性的主要环境因子, 为横断山南齿突蟾属物种种群动态监测、资源管理以及保护工作提供科学参考, 为相关生态学研究提供科学基础。

1 数据与方法 1.1 研究区概况

以《中国生物多样性保护优先区域范围》划定的横断山南段生物多样性保护优先区域为研究区, 横断山南总面积133656.15 km2, 包括西藏自治区的东部、四川省的西南部以及云南省西北部等区域。气候以南亚季风气候为主, 水系众多, 主要植被类型为灌丛植被, 地形起伏较大, 分布有较多高山深谷, 北部地区海拔较高, 南部海拔较低[ 25]

1.2 齿突蟾属在横断山南生物多样性保护优先区域分布数据

齿突蟾属分布数据主要有3种来源:全球生物多样性信息网络( https://www.gbif.org/)、公开发表的文献、出版物以及横断山南生物多样性保护优先区域的实地调查数据, 其中只选择有精确经纬度信息的分布数据, 并且去除错误、重复以及不在研究区内的分布信息。为了降低物种分布点的空间自相关[ 26], 对物种数据进行筛选, 利用ArcGIS 10.2.2中的SDMtoolbox插件使相距在1 km以内的物种分布数据只保留一个[ 27], 最终获得了36条西藏齿突蟾分布数据、23条胸腺齿突蟾分布数据、26条刺胸齿突蟾分布数据、10条圆疣齿突蟾分布数据以及5条贡山齿突蟾分布数据( 图 1)。

图 1 5种齿突蟾属物种的出现点 Fig. 1 Occurrence records of five species of Scutiger
图选项
1.3 环境变量

本研究总共选取了24个环境变量, 包含19个生物气候变量、3个地形变量(海拔、坡度、坡向)、2个栖息地变量(距河流距离、植被指数)。生物气候变量数据和海拔数据来源于WorldClim数据库( https://www.worldclim.org)中的1970-2000年监测环境数据, 空间分辨率为1 km, 并从海拔数据中提取出坡度、坡向。为便于后续建模分析, 将坡向数据减180°, 并取绝对值, 用以表示与正南阳坡的接近程度。河网数据来源于OpenStreetMap平台( https://www.openstreetmap.org/), 该数据能够全面反映区域的水系分布, 距河流距离通过欧式距离分析得到。植被指数来源于中国科学院资源环境科学与数据中心( https://www.resdc.cn/)的中国年度植被指数空间分布数据集, 经过处理得到2008-2018年平均归一化植被指数值, 分辨率为1 km。由于环境变量共线性会影响模型预测精度, 并夸大部分环境变量作用, 因此需要去除存在高度共线性的环境变量[ 28]。本研究利用方差膨胀因子(VIF)检测环境变量的共线性, 并去除方差膨胀因子大于10的环境变量, 最终剩余9个环境变量参与后续建模( 表 1)。

表 1 建模所用环境变量 Table 1 The environment variables for modeling
环境变量
Environmental variables
描述
Description
环境变量
Environmental variables
描述
Description
Bio3 等温性 Aspect 坡向
Bio5 最暖月最高温 Slope 坡度
Bio7 温度年较差 Distance to river 距河流距离
Bio13 最湿月降水量 NDVI 2008—2018年平均归一化植被指数值
Bio19 最冷季降水量
表选项
1.4 模型创建与评估 1.4.1 Maxent模型创建

由于物种采样强度在空间上往往不均衡, 存在采样偏差, 因此本研究利用SDMtoolbox插件创建高斯核密度偏差栅格文件[ 27], 使得在采样偏差大的区域选择更多的背景点。为避免训练数据集和验证数据集存在空间自相关, 导致模型预测性能被高估, 因此采用屏蔽地理结构分区法划分训练数据和测试数据[ 26], 将物种分布点在空间上等分为4块, 依次使用其中的三块构建Maxent模型, 剩余的一块验证模型, 共建立4个Maxent模型。由于圆疣齿突蟾和贡山齿突蟾的物种分布点数量较少, 不足以支持在空间上分区, 因此, 为了充分利用物种数据点, 使用留一交叉验证法对圆疣齿突蟾和贡山齿突蟾进行数据划分[ 29], 最终为圆疣齿突蟾建立了10个Maxent模型, 贡山齿突蟾建立了5个Maxent模型, 最终分析结果为各齿突蟾属物种所有模型的平均值。以上两种数据分区方法都是通过调用ENMeval包完成的[ 30]

已有研究表明, Maxent模型在默认设置下可能并不是最佳模型, 同时在样本量较小的情况下, Maxent模型更容易过拟合, 因此有必要进行模型参数调整[ 26, 31]。本研究通过调用SDMtune包[ 32], 将正则化系数设置为0.1-5, 间隔为0.1, 特征函数组合设置为L、LH、LQP、LQPH、LQPHT六种, 以AUCdiff(训练集AUC与测试集AUC的差值)和测试集AUC值作为指示, 利用遗传算法寻找Maxent模型的最佳正则化系数和特征函数组合设置, 并以最佳参数设置的Maxent模型输出结果作为最终结果。

利用Natural Breaks法将5种齿突蟾属物种潜在分布概率划分为不适宜区、低适宜区、中适宜区、高适宜区四类潜在分布区, 并且计算得到不同适宜等级栖息地面积和占比。其中, 对于西藏齿突蟾, 不适宜区(0-0.13), 低适宜区(0.13-0.34), 中适宜区(0.34-0.62), 高适宜区(0.62-1);对于刺胸齿突蟾, 不适宜区(0-0.13), 低适宜区(0.13-0.33), 中适宜区(0.33-0.59), 高适宜区(0.59-1);对于胸腺齿突蟾, 不适宜区(0-0.14), 低适宜区(0.14-0.35), 中适宜区(0.35-0.61), 高适宜区(0.61-1);对于圆疣齿突蟾, 不适宜区(0-0.08), 低适宜区(0.08-0.27), 中适宜区(0.27-0.59), 高适宜区(0.59-1);对于贡山齿突蟾, 不适宜区(0-0.08), 低适宜区(0.08-0.27), 中适宜区(0.27-0.55), 高适宜区(0.55-1)。

1.4.2 模型表现评估

本研究使用接受者操作特征曲线下面积(AUC)来评估模型表现, 接受者操作特征曲线(ROC)是以假阳性率为横坐标, 以真阳性率为纵坐标, 由绘制出的所有可能阈值点连接而成的曲线[ 11]。AUC值是ROC曲线与横坐标围成的面积, 它提供了一个不依赖于特定选择阈值的单一模型性能评估方法[ 11], AUC值的取值范围为0-1, AUC值越接近1, 模型的预测效果越好, 评估标准为:0.7-0.8可用, 0.8-0.9良好, 0.9-1.0优秀[ 16]

1.5 环境因子重要性分析

使用Maxent模型输出的环境变量贡献率和刀切法检验来评估环境因子的重要性, 其中刀切法检验是通过每次仅使用某一环境变量或排除某一环境变量而导致模型训练AUC值的变化来确定重要环境变量。通过创建环境变量响应曲线进一步分析重要环境因子对齿突蟾属物种存在概率的影响, 得到5种齿突蟾属物种潜在分布区的环境特征。

2 结果与分析 2.1 模型最佳参数设置与表现评估

根据AUCdiff和测试AUC值两个指标, 选择5种齿突蟾属物种最佳参数设置下的Maxent模型, 并进一步计算最佳Maxent模型的AUC值, 分析模型预测性能( 表 2)。与默认设置相比, 5种齿突蟾的最佳Maxent模型的测试AUC值均有所提高。最佳Maxent模型模拟西藏齿突蟾、刺胸齿突蟾、胸腺齿突蟾、圆疣齿突蟾、贡山齿突蟾潜在地理分布时的训练AUC值、测试AUC值均大于0.7, 表明Maxent模型对5种齿突蟾属物种分布的预测结果较准确可信。其中Maxent模型模拟贡山齿突蟾、圆疣齿突蟾潜在地理分布时的效果最好, 训练AUC值、测试AUC值均大于0.9, 预测结果精度极准确。

表 2 Maxent模型的最佳设置与预测性能 Table 2 The performance of Maxent model and best settings
物种
Species
最佳正则化系数
The best regularization multiplier
最佳特征函数组合
The best feature combination
最佳参数设置下训练集AUC值
The training AUC based on best settings
最佳参数设置下测试集AUC值
The test AUC based on best settings
最佳参数设置下AUCdiff
The AUCdiff based on best settings
默认设置下训练集AUC值
The training AUC based on default settings
默认设置下测试集AUC值
The test AUC based on default settings
西藏齿突蟾
Scutiger boulengeri
1.9 LH 0.885 0.786 0.098 0.906 0.774
刺胸齿突蟾
Scutiger mammatus
0.1 LQP 0.932 0.804 0.128 0.917 0.753
胸腺齿突蟾
Scutiger glandulatus
0.2 LQP 0.909 0.737 0.172 0.916 0.697
贡山齿突蟾
Scutiger gongshanensis
0.3 LQP 0.987 0.935 0.052 0.987 0.911
圆疣齿突蟾
Scutiger tuberculatus
0.6 LH 0.971 0.902 0.069 0.955 0.894
AUC:接受者操作特征曲线下面积Area Under Curve;L:线性linear;H:片段化hinge;Q:二次型quadratic;P:乘积型product
表选项
2.2 5种齿突蟾属物种的潜在分布区

5种齿突蟾属物种的潜在分布图表明( 图 2), 西藏齿突蟾的潜在分布区主要位于横断山南的北部地区, 高适生区主要分布在西藏芒康和四川甘孜州境内, 另外四川木里、云南德钦也有分布。刺胸齿突蟾的潜在分布区主要位于横断山南的中部、西北部地区, 高适生区主要分布在西藏芒康、云南德钦、香格里拉、四川得荣、乡城、巴塘等地区。胸腺齿突蟾在横断山南的中部、西北部和北部地区有广泛分布, 高适生区主要位于云南迪庆州以及四川甘孜州境内, 在西藏芒康、左贡、察隅也有零星分布。圆疣齿突蟾主要分布在横断山南东北部的四川省境内, 高适生区主要位于四川越西、美姑、甘洛、雷波、昭觉、冕宁、石棉等地区。贡山齿突蟾的潜在分布区位于横断山南的西南部, 高适生区集中分布在云南省怒江州的福贡、贡山、泸水、兰坪4个县。5种齿突蟾属物种的不同适宜等级栖息地面积和占比显示( 表 3):西藏齿突蟾在横断山南的适宜分布区范围最大, 约73485.48 km2, 占横断山南总面积的54.98%;贡山齿突蟾在横断山南的适宜分布区范围最小, 仅为24442.42 km2, 占横断山南总面积的18.29%。胸腺齿突蟾在横断山南的高适宜区范围最大, 为12286.70 km2, 其次为西藏齿突蟾。圆疣齿突蟾和贡山齿突蟾的高适宜区范围较小, 仅占横断山南总面积的2%左右。

图 2 5种齿突蟾属物种在横断山南的潜在地理分布 Fig. 2 Potential geographical distribution of five species of Scutiger in the south of Hengduan Mountains Biodiversity Conservation Priority Zone
图选项

表 3 5种齿突蟾属物种的不同适生等级栖息地面积 Table 3 Areas of the different classes of habitat distribution of the five species of Scutiger
西藏齿突蟾
Scutiger boulengeri
刺胸齿突蟾
Scutiger mammatus
胸腺齿突蟾
Scutiger glandulatus
圆疣齿突蟾
Scutiger tuberculatus
贡山齿突蟾
Scutiger gongshanensis
低适宜区面积
Area of low suitable habitats
39720.35 27931.81 33530.44 20038.47 14999.32
低适宜区占比
Percentage of low suitable habitats area in the south of Hengduan Mountains Biodiversity Conservation Priority Zone
0.2972 0.2090 0.2509 0.1499 0.1122
中适宜区面积
Area of moderate suitable habitats
21865.31 18176.04 20905.33 6520.75 6304.99
中适宜区占比
Percentage of moderate suitable habitats area in the south of Hengduan Mountains Biodiversity Conservation Priority Zone
0.1636 0.1360 0.1564 0.0488 0.0472
高适宜区面积
Area of high suitable habitats
11899.83 8111.38 12286.70 2734.58 3138.11
高适宜区占比
Percentage of high suitable habitats area in the south of Hengduan Mountains Biodiversity Conservation Priority Zone
0.0890 0.0607 0.0919 0.0205 0.0235
总的潜在分布区面积
Area of total suitable habitats
73485.48 54219.24 66722.46 29293.80 24442.42
总的潜在分布区占比
Percentage of total suitable habitats area in the south of Hengduan Mountains Biodiversity Conservation Priority Zone
0.5498 0.4057 0.4992 0.2192 0.1829
表选项
2.3 影响5种齿突蟾属物种潜在分布的主要环境因子

环境因子贡献率输出结果显示( 表 4), 最冷季降水量(Bio19)、温度年较差(Bio7)、距河流距离(Distance to river)是西藏齿突蟾模型模拟过程中贡献率排在前三位的环境因子, 对西藏齿突蟾潜在分布影响最大, 其中最冷季降水量的贡献率高达48.6%。最湿月降水量(Bio13)、最暖月最高温(Bio5)、距河流距离为刺胸齿突蟾模型模拟过程中贡献率排在前三位的环境因子, 贡献率皆超过10%, 对刺胸齿突蟾潜在分布影响较大。距河流距离、最冷季降水量、温度年较差是胸腺齿突蟾模型模拟过程中贡献率依次排在前三位的环境因子, 对胸腺齿突蟾潜在分布影响较大, 其中距河流距离对胸腺齿突蟾模型模拟的贡献率高达29.5%。归一化植被指数(NDVI)、等温性(Bio3)、最冷季降水量对圆疣齿突蟾模型模拟的贡献率排在前三位, 对圆疣齿突蟾的潜在分布影响较大, 其中归一化植被指数的贡献率高达41.17%。最冷季降水量和最暖月最高温对贡山齿突蟾模型模拟影响较大, 贡献率皆超过20%, 其中最冷季降水量对贡山齿突蟾潜在分布的贡献率高达64.66%。

表 4 5种齿突蟾属物种的环境因子贡献率/% Table 4 The environmental factors and their percentage contribution of the five species of Scutiger
环境变量
Environmental variables
西藏齿突蟾
Scutiger boulengeri
刺胸齿突蟾
Scutiger mammatus
胸腺齿突蟾
Scutiger glandulatus
圆疣齿突蟾
Scutiger tuberculatus
贡山齿突蟾
Scutiger gongshanensis
等温性Bio3 0.025 1.2 2.825 27.96 0.92
最暖月最高温Bio5 5.65 19 8.3 1.03 27.6
温度年较差Bio7 19.625 2.875 10.325 0 0.22
最湿月降水量Bio13 4.05 27.325 9.2 0.27 0.18
最冷季降水量Bio19 48.6 11.95 19.6 15.18 64.66
坡向Aspect 0.25 8.9 2.475 3.5 4.66
坡度Slope 5.875 8.05 10.15 9.12 0.36
距河流距离Distance to river 10.05 14.7 29.5 1.79 1.14
归一化植被指数NDVI 5.9 6 7.575 41.17 0.22
表选项

刀切法检验结果显示( 图 3), 在使用单一环境变量构建Maxent模型时, 温度年较差(Bio7)、最冷季降水量(Bio19)、最湿月降水量(Bio13)在西藏齿突蟾Maxent模型模拟时的训练AUC值最高, 是影响西藏齿突蟾潜在分布的关键因子。最湿月降水量(Bio13)、归一化植被指数(NDVI)、最暖月最高温(Bio5)在单独用于刺胸齿突蟾Maxent模型模拟时获得了最高的训练AUC值, 对刺胸齿突蟾潜在分布有重要影响。最暖月最高温(Bio5)、距河流距离(Distance to river)、最湿月降水量(Bio13)在单独用于胸腺齿突蟾模型构建时获得了最高的训练AUC值, 是影响胸腺齿突蟾潜在分布的重要因子。单独使用归一化植被指数(NDVI)、最冷季降水量(Bio19)、等温性(Bio3)构建圆疣齿突蟾Maxent模型时, 获得了最高的训练AUC值, 它们是影响圆疣齿突蟾分布格局的重要因子。最冷季降水量(Bio19)、温度年较差(Bio7)、最湿月降水量(Bio13)在单独用于贡山齿突蟾模型模拟时的训练AUC值最高, 对贡山齿突蟾潜在分布有重要影响。

图 3 5种齿突蟾属物种基于训练AUC值的刀切法检验结果 Fig. 3 Jackknife test based on train AUC of five species of Scutiger NDVI:归一化植被指数Normalized Difference Vegetation Index;Bio7:温度年较差Temperature annual range;Bio5:最暖月最高温Max temperature of the warmest month;Bio3:等温性Isothermality;Bio19:最冷季降水量Precipitation of the coldest quarter;Bio13:最湿月降水量Precipitation of the wettest month; AUC: 曲线下方面积Area under curve
图选项

综合环境变量贡献率和刀切法检验分析结果, 判断出影响5种齿突蟾属物种潜在分布的最主要环境因子:温度年较差(Bio7)、最冷季降水量(Bio19)是影响西藏齿突蟾潜在分布的最主要因子;最湿月降水量(Bio13)、最暖月最高温(Bio5)是影响刺胸齿突蟾潜在分布的最主要因子;距河流距离是影响胸腺齿突蟾潜在分布的最主要因子;归一化植被指数(NDVI)、最冷季降水量、等温性(Bio3)是影响圆疣齿突蟾潜在分布的最主要因子;最冷季降水量是影响贡山齿突蟾潜在分布的最主要因子。

2.4 5种齿突蟾属物种潜在分布区的环境特征

一般认为存在概率大于0.5时, 其对应的环境变量特征适合物种生存[ 33], 根据环境变量响应曲线, 西藏齿突蟾的生存受到温度年较差和最冷季降水量两个环境变量的限制, 适宜生存在温度年较差29.1-34.2℃、最冷季降水量6-12 mm的区域, 其中温度年较差越高, 最冷季降水量越低, 西藏齿突蟾的存在概率越高。刺胸齿突蟾适宜生存在最湿月降水量119-148 mm、最暖月最高温14.7-23.2℃的区域, 其中最湿月降水量值越低, 刺胸齿突蟾存在概率越高。胸腺齿突蟾的生存受到距河流距离的限制, 适宜生存在距河流距离小于3.7 km左右的区域。圆疣齿突蟾偏好生存在植被指数为0.71-0.8、最冷季降水量15.9-29.4 mm、等温性37.9-43.3的区域。贡山齿突蟾偏好生存在最冷季降水量为65-129 mm的区域, 其中, 最冷季降水量越丰富, 贡山齿突蟾的存在概率越高。

2.5 横断山南齿突蟾属物种丰富度分布格局

以训练数据特异性和敏感性加和最大判断出阈值[ 34], 通过阈值将Maxent模型输出的物种栖息地适宜性图转换为物种二元分布图, 并通过栅格计算器进一步统计出每个栅格内齿突蟾属物种总数, 获得横断山南齿突蟾属物种丰富度分布格局( 图 4)。横断山南齿突蟾属物种丰富度分布格局表明, 横断山南北部地区的齿突蟾属物种丰富度明显高于南部地区, 齿突蟾属物种丰富度较高的区域主要位于北部, 包括西藏芒康、四川稻城、理塘、巴塘、得荣等地。横断山南中部地区的齿突蟾属物种丰富度也相对较高, 主要是在云南香格里拉附近。

图 4 横断山南齿突蟾属丰富度分布格局 Fig. 4 The distribution pattern of potential species richness of Scutiger in the south of Hengduan Mountains Biodiversity Conservation Priority Zone
图选项
3 讨论 3.1 模型的合理性

由于Maxent模型具有对环境变量共线性的敏感性较低[ 12]、样本量较少时仍能保持稳健[ 13]、能够拟合复杂的变量关系[ 35]等优势, 近年来Maxent模型已被广泛应用于物种分布模拟研究中[ 5, 14]。然而, 有研究指出[ 36], Maxent模型中复杂的函数关系容易导致过拟合, 通常通过调整参数来优化Maxent模型[ 26, 31]。正则化系数和特征函数组合是优化Maxent模型的两个重要参数, ENMeval包[ 30]、Kuenm包[ 37]、SDMtune包[ 32]等相关开源软件可以有效地寻找Maxent模型的最佳参数设置。本研究采用SDMtune包寻找Maxent模型的最佳参数设置, 与ENMeval包、Kuenm包等相比, SDMtune包采用遗传算法寻找最佳参数, 极大的节省了计算机运行时间。在本研究中, 与默认参数设置相比, 5种齿突蟾最佳Maxent模型的测试集AUC均有所增加, 表明通过寻找最佳参数设置, Maxent模型的预测性能提升, 同时减弱了模型过拟合可能性。因此, 通过评估潜在的Maxent参数组合, 选择最佳的Maxent模型是有效且必要的。

模型预测性能评估是物种分布模型模拟过程中的重要步骤, 最理想的评估数据是与训练数据在时间和空间上完全独立的数据集[ 28], 然而受现实因素制约, 这种数据往往难以获取, 因此交叉验证成为一种较有效的替代工具。常见的交叉验证法由于没有保证训练数据与验证数据的空间独立性, 往往会高估模型的预测性能[ 38- 39], 因此本研究采用空间上独立的交叉验证法进行模型性能评估, 但这种方法也不能完全保证训练数据和验证数据在空间上独立, 因为相邻的两个数据分区仍可能存在空间自相关, 同时训练数据和验证数据可能具有相同的采样偏差, 这都会导致本研究中Maxent模型预测性能被高估。由于圆疣齿突蟾和贡山齿突蟾的样本量较小, 不足以进行空间分区, 因此选择留一交叉验证法进行模型性能评估[ 29], 保证物种数据被充分使用。尽管空间上独立的交叉验证法、留一交叉验证法能够较好的评估5种齿突蟾Maxent模型的预测性能, 但与自举法、重复分裂采样的交叉验证法相比, Maxent模型的重复次数过少, 模型的稳健性可能受到影响。

物种分布模型的预测性能与建模对象有关, 研究认为生态位狭窄的物种往往能更好的被预测[ 13], 在本研究中, Maxent模型对圆疣齿突蟾和贡山齿突蟾两个分布范围狭窄的物种的预测精度极准确, 训练AUC值、测试AUC值均大于0.9, 但这也可能是由于参与建模的物种分布点较少。样本量也会影响物种分布模型的预测性能[ 40], 样本量越大, 模型预测性能越好。在本研究中, 贡山齿突蟾只有5个物种分布点, 但这并不是因为在研究区内采样不充分, 没有完全捕获物种分布的实际生态位, 而是因为物种本身分布范围狭窄, 在这种情况下, 即使增加物种分布点, 也只能增加已有区域采样点的密度, 并不能提升模型预测性能[ 16]

3.2 5种齿突蟾属物种的潜在地理分布

横断山区可能是齿突蟾属物种的起源中心和分化中心[ 20, 22], 分布在横断山区的齿突蟾属物种丰富, 特有种多, 本文研究结果表明横断山南北部地区的齿突蟾属物种丰富度明显高于南部地区, 这可能是受海拔的限制, 横断山南北部海拔较高, 南部海拔较低, 而齿突蟾属偏好生存在海拔3000-4500 m左右的地区[ 20]

在本研究中, 5种齿突蟾属物种潜在地理分布格局存在差异, 西藏齿突蟾的潜在分布区主要位于横断山南的北部, 圆疣齿突蟾分布在横断山南东北部的四川省境内, 贡山齿突蟾的潜在分布区主要位于横断山南的西南部地区, 这可能是受到最冷季降水量的影响, 因为横断山南北部、东北部、西南部的最冷季降水量具有较明显的差异, 而在本研究中, 最冷季降水量是影响西藏齿突蟾、贡山齿突蟾、圆疣齿突蟾潜在分布的最主要环境因素。刺胸齿突蟾和胸腺齿突蟾的潜在分布格局较为相似, 在横断山南的中部、西北部地区都有较多分布, 但胸腺齿突蟾的潜在分布区更为碎片化, 可能是因为胸腺齿突蟾的潜在分布受距河流距离影响较大, 而横断山南水系发达, 河流众多, 距河流距离变化幅度较大。

5种齿突蟾属物种潜在分布区大小存在差异, 贡山齿突蟾和圆疣齿突蟾的潜在分布区最为狭窄, 这与它们的种群状态较一致, 贡山齿突蟾和圆疣齿突蟾都是中国特有种, 种群数量较少, 受威胁等级为VU。以往研究表明[ 20], 贡山齿突蟾实际分布区位于云南贡山、碧罗雪山, 圆疣齿突蟾实际分布区位于四川越西、昭觉、冕宁、西昌、盐边, Maxent模型预测结果表明, 与实际分布区相比, 贡山齿突蟾、圆疣齿突蟾的潜在分布区有扩大的趋势。特别是对于圆疣齿突蟾, 它的潜在分布区零星出现在云南香格里拉, 但目前该地区关于圆疣齿突蟾的记录较少, 仅在Subba等人的研究中出现[ 41]。鉴于云南香格里拉的归一化植被指数、最冷季降水量、等温性与圆疣齿突蟾主要栖息地的环境条件较为相似, 我们认为圆疣齿突蟾可能在云南香格里拉零星出现。

3.3 影响齿突蟾属物种潜在分布的主要环境因子及其特征

5种齿突蟾属物种潜在地理分布的差异反映出它们对不同栖息地环境的偏好。湿度是影响两栖动物空间分布的主要因素, 降水会改变两栖动物的繁殖场所, 影响两栖动物的生长发育速度、免疫功能以及存活率[ 42- 43]。在本研究中, 齿突蟾属物种潜在分布受降水因子影响较大, 但不同齿突蟾属物种对水分的偏好不同, 西藏齿突蟾适宜生存在最冷季降水量6-12 mm的区域, 圆疣齿突蟾适宜生存在最冷季降水量15.9-29.4 mm的区域, 贡山齿突蟾适宜生存在最冷季降水量为65-129 mm的区域, 这可能是导致西藏齿突蟾、贡山齿突蟾、圆疣齿突蟾潜在地理分布格局差异较大的主要原因。温度是影响两栖动物生理和行为的主要环境因子[ 44], 对两栖动物的繁殖和生境选择有重要影响。西藏齿突蟾适宜生存在温度年较差29.1-34.2℃的区域, 刺胸齿突蟾适宜生存在最暖月最高温14.7-23.2℃的区域, 圆疣齿突蟾适宜生存在等温性37.9-43.3的区域。栖息地特征也会影响物种潜在分布, 胸腺齿突蟾的生存受距河流距离的限制, 适宜生存在距河流距离小于3.7 km左右的区域, 这与胸腺齿突蟾本身的生物学特性一致, 它多栖息于中、小型山溪边[ 20], 偏好湿润环境。费梁等人[ 20]发现圆疣齿突蟾生活在山区林木茂密的流溪及溪岸两侧, 本研究发现圆疣齿突蟾偏好生存在植被指数为0.71-0.8的区域, 与之结果较为一致。

本文通过优化Maxent模型首次预测5种高海拔齿突蟾属物种在横断山南生物多样性保护优先区的潜在分布情况, 为横断山南齿突蟾属物种种群动态监测、资源管理、保护工作提供科学参考, 应重点关注齿突蟾属物种丰富度较高的区域以及受威胁齿突蟾属物种的高适宜区。另外, 本文仍存在不足, 如没有考虑部分齿突蟾属物种不同支系间的生态位可能存在差异[ 45], 而是认为种内生态位一致, 因此在后续研究中需要进一步考虑广布种不同支系间生态位的差异。

参考文献
[1]
Carvalho T, Becker C G, Toledo L F. Historical amphibian declines and extinctions in Brazil linked to chytridiomycosis. Proceedings of the Royal Society B: Biological Sciences, 2017, 284(1848): 20162254. DOI:10.1098/rspb.2016.2254
[2]
Hof C, Araújo M B, Jetz W, Rahbek C. Additive threats from pathogens, climate and land-use change for global amphibian diversity. Nature, 2011, 480(7378): 516-519. DOI:10.1038/nature10650
[3]
Becker C G, Fonseca C R, Haddad C F B, Batista R F, Prado P I. Habitat split and the global decline of amphibians. Science, 2007, 318(5857): 1775-1777. DOI:10.1126/science.1149374
[4]
Cohen J M, Civitello D J, Venesky M D, McMahon T A, Rohr J R. An interaction between climate change and infectious disease drove widespread amphibian declines. Global Change Biology, 2019, 25(3): 927-937. DOI:10.1111/gcb.14489
[5]
Holder A M, Markarian A, Doyle J M, Olson J R. Predicting geographic distributions of fishes in remote stream networks using maximum entropy modeling and landscape characterizations. Ecological Modelling, 2020, 433: 109231. DOI:10.1016/j.ecolmodel.2020.109231
[6]
Roy-Dufresne E, Saltré F, Cooke B D, Mellin C, Mutze G, Cox T, Fordham D A. Modeling the distribution of a wide-ranging invasive species using the sampling efforts of expert and citizen scientists. Ecology and Evolution, 2019, 9(19): 11053-11063. DOI:10.1002/ece3.5609
[7]
Li P X, Zhu W Q, Xie Z Y, Qiao K. Integration of multiple climate models to predict range shifts and identify management priorities of the endangered Taxus wallichiana in the Himalaya-Hengduan Mountain region. Journal of Forestry Research, 2020, 31(6): 2255-2272. DOI:10.1007/s11676-019-01009-5
[8]
Li R Q, Xu M, Wong M H G, Qiu S, Sheng Q K, Li X H, Song Z M. Climate change-induced decline in bamboo habitats and species diversity: implications for giant panda conservation. Diversity and Distributions, 2015, 21(4): 379-391. DOI:10.1111/ddi.12284
[9]
张华, 赵浩翔, 王浩. 基于Maxent模型的未来气候变化情景下胡杨在中国的潜在地理分布. 生态学报, 2020, 40(18): 6552-6563.
[10]
杨蕾, 杨立, 李婧昕, 张超, 霍兆敏, 栾晓峰. 东北地区5个物种潜在栖息地变化与优化保护规划. 生态学报, 2019, 39(3): 1082-1094.
[11]
Phillips S J, Anderson R P, Schapire R E. Maximum entropy modeling of species geographic distributions. Ecological Modelling, 2006, 190(3/4): 231-259.
[12]
Elith J, Phillips S J, Hastie T, Dudík M, Chee Y E, Yates C J. A statistical explanation of Maxent for ecologists. Diversity and Distributions, 2011, 17(1): 43-57. DOI:10.1111/j.1472-4642.2010.00725.x
[13]
Guisan A, Zimmermann N E, Elith J, Graham C H, Phillips S, Peterson A T. What matters for predicting the occurrences of trees: techniques, data, or species' characteristics?. Ecological Monographs, 2007, 77(4): 615-630. DOI:10.1890/06-1060.1
[14]
Alatawi A S, Gilbert F, Reader T. Modelling terrestrial reptile species richness, distributions and habitat suitability in Saudi Arabia. Journal of Arid Environments, 2020, 178: 104153. DOI:10.1016/j.jaridenv.2020.104153
[15]
Wiltshire K H, Tanner J E. Comparing maximum entropy modelling methods to inform aquaculture site selection for novel seaweed species. Ecological Modelling, 2020, 429: 109071. DOI:10.1016/j.ecolmodel.2020.109071
[16]
Zank C, Becker F G, Abadie M, Baldo D, Maneyro R, Borges-Martins M. Climate change and the distribution of neotropical red-bellied toads (Melanophryniscus, Anura, Amphibia): how to prioritize species and populations?. PLoS One, 2014, 9(4): e94625. DOI:10.1371/journal.pone.0094625
[17]
黄勇杰, 卢佳斌, 王锋堂, 林英华, 刘磊, 米红旭, 莫方群, 方精, 李佳灵. 基于Maxent模型预测海南岛海南臭蛙的潜在地理分布. 动物学杂志, 2017, 52(1): 30-41.
[18]
张凯龙, 杨坤, 沃钰斌, 童浩杰, 金园庭. 基于Maxent模型的镇海林蛙种群潜在地理适宜性评价. 生态学杂志, 2018, 37(1): 164-170.
[19]
王剀, 任金龙, 陈宏满, 吕植桐, 郭宪光, 蒋珂, 陈进民, 李家堂, 郭鹏, 王英永, 车静. 中国两栖、爬行动物更新名录. 生物多样性, 2020, 28(2): 189-218.
[20]
费梁, 叶昌媛, 江建平. 中国两栖动物及其分布彩色图鉴. 成都: 四川科学技术出版社, 2012.
[21]
万丽霞, 李宏伟, 孙立新, 窦静莉. 我国齿突蟾属(两栖纲: 无尾目)资源现状及保护对策. 动物学杂志, 2014, 49(3): 357-365.
[22]
魏刚, 李子忠, 江建平. 角蟾科Megophryidae系统研究进展. 四川动物, 2010, 29(5): 652-658.
[23]
中国科学院中国动物志编辑委员会. 中国动物志: 两栖纲(中卷)无尾目. 北京: 科学出版社, 2009.
[24]
中国科学院青藏高原综合科学考察队. 横断山区两栖爬行动物. 北京: 科学出版社, 1997.
[25]
李俊生, 靳勇超, 王伟, 赵志平, 吴晓莆. 中国陆域生物多样性保护优先区域. 北京: 科学出版社, 2016.
[26]
Radosavljevic A, Anderson R P. Making better Maxent models of species distributions: complexity, overfitting and evaluation. Journal of Biogeography, 2014, 41(4): 629-643. DOI:10.1111/jbi.12227
[27]
Brown J L. SDMtoolbox: a python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. Methods in Ecology and Evolution, 2014, 5(7): 694-700. DOI:10.1111/2041-210X.12200
[28]
Guisan A, Thuiller W, Zimmermann N E. Habitat Suitability and Distribution Models: With Applications in R. Cambridge: Cambridge University Press, 2017.
[29]
Shcheglovitova M, Anderson R P. Estimating optimal complexity for ecological niche models: a jackknife approach for species with small sample sizes. Ecological Modelling, 2013, 269: 9-17. DOI:10.1016/j.ecolmodel.2013.08.011
[30]
Muscarella R, Galante P J, Soley-Guardia M, Boria R A, Kass J M, Uriarte M, Anderson R P. ENMeval: an R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in Ecology and Evolution, 2014, 5(11): 1198-1205. DOI:10.1111/2041-210X.12261
[31]
Morales N S, Fernández I C, Baca-González V. MaxEnt's parameter configuration and small samples: are we paying attention to recommendations? A systematic review. PeerJ, 2017, 5: e3093. DOI:10.7717/peerj.3093
[32]
Vignali S, Barras A G, Arlettaz R, Braunisch V. SDMtune: an R package to tune and evaluate species distribution models. Ecology and Evolution, 2020, 10(20): 11488-11506. DOI:10.1002/ece3.6786
[33]
郭杰, 刘小平, 张琴, 张东方, 谢彩香, 刘霞. 基于Maxent模型的党参全球潜在分布区预测. 应用生态学报, 2017, 28(3): 992-1000.
[34]
Liu C R, Newell G, White M. On the selection of thresholds for predicting species occurrence with presence-only data. Ecology and Evolution, 2016, 6(1): 337-348. DOI:10.1002/ece3.1878
[35]
Elith J, Graham C H, Anderson R P, Dudík M, Ferrier S, Guisan A, Hijmans R J, Huettmann F, Leathwick J R, Lehmann A, Li J, Lohmann L G, Loiselle B A, Manion G, Moritz C, Nakamura M, Nakazawa Y, Overton J M M, Peterson A T, Phillips S J, Richardson K, Scachetti-Pereira R, Schapire R E, Soberón J, Williams S, Wisz M S, Zimmermann N E. Novel methods improve prediction of species' distributions from occurrence data. Ecography, 2006, 29(2): 129-151. DOI:10.1111/j.2006.0906-7590.04596.x
[36]
Peterson A T, Papeş M, Eaton M. Transferability and model evaluation in ecological niche modeling: a comparison of GARP and Maxent. Ecography, 2007, 30(4): 550-560. DOI:10.1111/j.0906-7590.2007.05102.x
[37]
Cobos M E, Peterson A T, Barve N, Osorio-Olvera L. Kuenm: an R package for detailed development of ecological niche models using Maxent. PeerJ, 2019, 7: e6281. DOI:10.7717/peerj.6281
[38]
Wenger S J, Olden J D. Assessing transferability of ecological models: an underappreciated aspect of statistical validation. Methods in Ecology and Evolution, 2012, 3(2): 260-267. DOI:10.1111/j.2041-210X.2011.00170.x
[39]
Čengić M, Rost J, Remenska D, Janse J H, Huijbregts M A J, Schipper A M. On the importance of predictor choice, modelling technique, and number of pseudo-absences for bioclimatic envelope model performance. Ecology and Evolution, 2020, 10(21): 12307-12317. DOI:10.1002/ece3.6859
[40]
Mi C R, Huettmann F, Guo Y M, Han X S, Wen L J. Why choose Random Forest to predict rare species distribution with few samples in large undersampled areas? Three Asian crane species models provide supporting evidence. PeerJ, 2017, 5: e2849. DOI:10.7717/peerj.2849
[41]
Subba B, Sen S, Ravikanth G, Nobis M P. Direct modelling of limited migration improves projected distributions of Himalayan amphibians under climate change. Biological Conservation, 2018, 227: 352-360. DOI:10.1016/j.biocon.2018.09.035
[42]
Brannelly L A, Ohmer M E B, Saenz V, Richards-Zawacki C L. Effects of hydroperiod on growth, development, survival and immune defences in a temperate amphibian. Functional Ecology, 2019, 33(10): 1952-1961. DOI:10.1111/1365-2435.13419
[43]
Chen Y H. Habitat suitability modeling of amphibian species in southern and central China: environmental correlates and potential richness mapping. Science China Life Sciences, 2013, 56(5): 476-484. DOI:10.1007/s11427-013-4475-3
[44]
Todd B D, Scott D E, Pechmann J H K, Gibbons J W. Climate change correlates with rapid delays and advancements in reproductive timing in an amphibian community. Proceedings of the Royal Society B: Biological Sciences, 2011, 278(1715): 2191-2197. DOI:10.1098/rspb.2010.1768
[45]
Lin X Q, Shih C, Hou Y M, Shu X X, Zhang M H, Hu J H, Jiang J P, Xie F. Climatic-niche evolution with key morphological innovations across clades within Scutiger boulengeri (Anura: Megophryidae). Ecology and Evolution, 2021, 11(15): 10353-10368. DOI:10.1002/ece3.7838

深圳SEO优化公司大庆品牌网站设计报价昌都seo网站优化报价自贡网站推广工具报价金华SEO按天扣费哪家好廊坊网站优化推广哪家好长治SEO按天计费多少钱乐山关键词排名推荐张家口网站搜索优化哪家好长沙百度标王哪家好鞍山网站制作价格厦门企业网站制作哪家好湛江百度竞价包年推广哪家好成都品牌网站设计哪家好无锡网络广告推广百色seo网站推广公司太原seo哪家好商丘网站优化按天扣费多少钱新乡seo排名价格厦门模板制作价格德州网站seo优化多少钱百度网站优化多少钱漯河网站优化报价庆阳企业网站设计公司锦州如何制作网站松原网站推广工具舟山SEO按天计费哪家好龙岩建设网站推荐茂名网站开发多少钱哈密建站多少钱伊犁seo网站优化歼20紧急升空逼退外机英媒称团队夜以继日筹划王妃复出草木蔓发 春山在望成都发生巨响 当地回应60岁老人炒菠菜未焯水致肾病恶化男子涉嫌走私被判11年却一天牢没坐劳斯莱斯右转逼停直行车网传落水者说“没让你救”系谣言广东通报13岁男孩性侵女童不予立案贵州小伙回应在美国卖三蹦子火了淀粉肠小王子日销售额涨超10倍有个姐真把千机伞做出来了近3万元金手镯仅含足金十克呼北高速交通事故已致14人死亡杨洋拄拐现身医院国产伟哥去年销售近13亿男子给前妻转账 现任妻子起诉要回新基金只募集到26元还是员工自购男孩疑遭霸凌 家长讨说法被踢出群充个话费竟沦为间接洗钱工具新的一天从800个哈欠开始单亲妈妈陷入热恋 14岁儿子报警#春分立蛋大挑战#中国投资客涌入日本东京买房两大学生合买彩票中奖一人不认账新加坡主帅:唯一目标击败中国队月嫂回应掌掴婴儿是在赶虫子19岁小伙救下5人后溺亡 多方发声清明节放假3天调休1天张家界的山上“长”满了韩国人?开封王婆为何火了主播靠辱骂母亲走红被批捕封号代拍被何赛飞拿着魔杖追着打阿根廷将发行1万与2万面值的纸币库克现身上海为江西彩礼“减负”的“试婚人”因自嘲式简历走红的教授更新简介殡仪馆花卉高于市场价3倍还重复用网友称在豆瓣酱里吃出老鼠头315晚会后胖东来又人满为患了网友建议重庆地铁不准乘客携带菜筐特朗普谈“凯特王妃P图照”罗斯否认插足凯特王妃婚姻青海通报栏杆断裂小学生跌落住进ICU恒大被罚41.75亿到底怎么缴湖南一县政协主席疑涉刑案被控制茶百道就改标签日期致歉王树国3次鞠躬告别西交大师生张立群任西安交通大学校长杨倩无缘巴黎奥运

深圳SEO优化公司 XML地图 TXT地图 虚拟主机 SEO 网站制作 网站优化