珍稀濒危植物伯乐树的潜在适生区预测

龚 维，夏 青，陈红锋，俞新华，伍 菲

(1 华南农业大学 生命科学学院，广东 广州510642;2 中国科学院 华南植物园/中国科学院 植物资源保护和可持续利用重点实验室，广东 广州510650;3 华南农业大学公共基础课实验教学中心，广东 广州510642;4 21 世纪空间技术应用股份有限公司，北京100096)

物种地理分布格局对于研究生物多样性保护具有重要意义，也是制定物种保护策略的重要依据．亚热带常绿阔叶林是我国最重要的森林群落类型之一，蕴藏着丰富的植物类群，其地理分布格局对于探讨我国常绿阔叶林植物响应全球气候变化具有重要意义．物种分布模型常用于物种潜在分布区预测，基于气候匹配法的物种分布模型已得到了不断的完善和发展．目前，最大熵算法模型(Maximum entropy algorithm modeling program，Maxent)在珍稀濒危动植物的潜在分布区预测中得到了广泛应用［1-6］．

伯乐树Bretschneidera sinensis 是伯乐树科唯一的一个种，也是我国的特有种，以我国为分布中心［7-9］．该种在我国主要分布于长江流域各省，多生于亚热带常绿阔叶林山地［9-10］．长期以来，伯乐树生境破坏严重，现存母树资源极为稀少．加上结实率低，天然更新困难，该种已处于珍稀濒危状态，被评估为濒危等级［11］，是我国一级重点保护野生植物，世界自然保护联盟(the International Union for Conservation of Nature，IUCN)红皮书已将其列为濒危种［12］．作为第三纪古热带植物区系的濒危孑遗种［13］，探讨其地理分布格局对气候变化的响应有利于开展该种的保护生物学研究．本研究对伯乐树进行生态位模型(Ecological niche modelling)构建和分析，根据该种已知地理分布点数据，基于气候相似性原理，利用Maxent 模型，结合地理分析软件，预测伯乐树当代及21 世纪70年代的潜在适生区．同时，结合该种野外调查资料、地理分布区海拔和植被数据，分析其地理分布格局及生态和气候影响因子，以期为将来进一步开展该种保护生物学研究、野外资源调查等积累数据、提供依据．

1 材料与方法

1．1 数据收集

伯乐树的物种分布数据主要来源于中国数字标本馆的记录(CVH)，包括中国科学院成都生物研究所、中国科学院西双版纳热带植物园、中国科学院广西植物研究所、中国科学院华南植物园、中国科学院昆明植物研究所、中国科学院庐山植物园以及西北农林科技大学等．同时，结合《中国植物红皮书:稀有濒危植物(第一册)》［13］、《中国植物志》［9］及各地方植物志．另外，查阅伯乐树野外调查数据、研究成果和文献资料［14-17］．最终，获得了全面的伯乐树地理分布数据，共13 个省份(地区)、114 个分布点．对于部分采集记录或相关文献资料中缺乏具体地理坐标的分布记录，根据百度地图，并借助Google Earth 软件，查找出相应的经纬度．

在伯乐树地理分布数据收集的基础上，利用DIVA-GIS 7.4.0 软件和ARCGIS 9.1 软件分别以海拔图层和植被覆盖图层作为底图绘制该种的地理分布图．海拔图层来自于DIVA-GIS 软件自带的地理信息数据库，植被覆盖图层来自于网站www．iscgm．org．所采用气候数据共包括19 个气候因子(表1)，来源于世界气候数据网站(http:∥www．worldclim．org)．当代气候数据为Worldclim 数据库1950—2000年的来自世界各地气象站的气候信息，21 世纪70年代气候数据为Worldclim 数据库校准后的GCMs 气候投影数据［18］．应用2.5 arcminutes 数据库，每个栅格单元大致相当于22 km2(表1)．分析用的底层地图是1∶400 万中国行政区划图，从国家基础地理信息系统网站(http:∥nfgis．nsdi．gov．cn/nfgis/chinese/c_xz．htm)下载．

1．2 数据处理

根据所获得的伯乐树经纬度数据，借助ARCGIS 9.1 软件，生成该种在中国境内的地理分布图，每个分布点至少有1 个标本采集或者文献的确切记录．分别以海拔图层和植被覆盖图层为底图，整合伯乐树地理分布点，利用ARCGIS 9.1 软件绘制该种的地理分布图．

将Worldclim 下载的气候数据通过DIVA-GIS 7.4.0 软件转换成ASCII格式，应用Maxent模型对伯乐树的潜在适生区进行预测．将地理分布点的经纬度数据和环境数据导入Maxent v3.3.3 软件，随机选取75%的数据作为训练集(Sample file)，用于建模，剩余25% 的分布点作为测试集(Test sample file)，用于验证．选择启发式估测和刀切法(Jackknife)，设置10 次重复，产生10 个预测随机模型，其它参数均为软件默认值．采取常用的ROC 曲线(Reciever operating characteristic curve)分析法预测精度，以AUC(Area under roc curve)值作为模型预测准确性的衡量指标，选择具有最高AUC 值的图层进行分布区预测分析．AUC 值越大，模型的预测能力越好．相应的标准是:AUC＜0.5，模型的预测能力低于随机模型的预测能力;AUC=0.5～0.7，模型的诊断价值低;AUC=0.7～0.9，模型的诊断价值中等;AUC=0.9～1.0，模型的诊断价值较高［1，19-20］．

模拟结果的输出选择ASCII 格式，在ARCGIS 9.1 软件中加载Maxent 运行结果，将结果转换成RASTER 格式．根据重分类程序将生境适应性按照自然划分法(Natural break)分成6 个等级(Ⅰ～Ⅵ)，其中第Ⅵ级为最适分布区．对各等级的图层进行统计分析，得到所研究的13 个省份或地区在各适生等级下的适生面积．另外，根据模拟结果中各气候因子的相对贡献率，筛选出影响伯乐树地理分布的主导气候因子．

2 结果与分析

2．1 伯乐树在中国的地理分布情况及潜在适生区预测

根据现有的伯乐树分布点数据，江西省(18个)、湖南省(18 个)、广东省(20 个)和贵州省(16个)的分布点最多．另外，福建省和广西壮族自治区分别有11 和10 个分布点，云南省和浙江省分别有6和7 个分布点，分布点最少的是四川省(4 个)、重庆市(1 个)、湖北省(1 个)、海南省(1 个)和台湾(1个)．其中，海南省的鹦哥岭是伯乐树在热带地区的唯一分布点．以海拔图层为底图绘制的伯乐树地理分布图(图1)表明，该种主要沿南岭山脉、罗霄山、武夷山分布，最常见于亚热带地区海拔500～2 000 m的中山地带，在海拔低于500 m 或近2 000 m的地区也有少量分布．以植被图层为底图绘制伯乐树地理分布图(图2)，结合野外调查数据，结果显示该种主要分布于常绿阔叶林、常绿针叶林和落叶阔叶林中，常散生在湿润沟谷、溪旁坡地的林丛和灌丛．部分群体出现在耕地边缘及其他自然植被的镶嵌区．根据Maxent 的输出结果，基于1950—2000年及21 世纪70年代的气候数据，训练集数据和测试集数据的AUC 值分别为0.994 和0.993．由图3 显示，基于1950—2000年气候数据，该种潜在适生区位于亚热带地区，20.0°～32.5°N、105.0°～122.5°E．其最适分布区位于长江以南山地，如广西猫儿山、大瑶山，广东北部山地，以及江西九连山等，沿南岭山脉以及罗霄山和幕阜山呈东西走向分布，在中国的总面积约为337 551 km2(表2)．基于21 世纪70年代气候数据，该种潜在适生区北移至35°N，在中国的最适分布区总面积减少为228 109 km2，川渝交界处新增为其最适分布区(表2 和图3)．

基于1950—2000年气候数据，该种在广西壮族自治区的最适分布区面积最大(74 592 km2)，其次是江西省(66 100 km2)，但是该种在广西壮族自治区现有地理分布点较少(10 个)．该种在广东省、贵州省、湖南省和福建省的最适分布区面积相当(25 336～37 161 km2)，在四川省、重庆市、浙江省、湖北省和云南省的最适分布区面积均低于10 000 km2，在海南省以及台湾均无最适分布区．基于21 世纪70年代气候数据，该种在江西省的最适分布区面积最大(78 680 km2)，其次是广西壮族自治区(66 488 km2)．该种在贵州省和福建省的最适分布区面积明显减小，低于10 000 km2，在浙江省、湖北省和云南省减少更为显著，均低于1 000 km2，但是，该种在川渝地区最适分布区面积均明显增加．

图1 伯乐树地理分布点海拔图Fig．1 Altitudes in the geographical distribution of Bretschneidera sinensis

图2 伯乐树地理分布的植被覆盖情况Fig．2 A landcover type in the geographical distribution of Bretschneidera sinensis

图3 基于Maxent 模型与19 个气候因子预测的伯乐树1950—2000年及21 世纪70年代的潜在适生区Fig．3 A prediction of the potential distribution range of Bretschneidera sinensis in 1950—2000 and the 2070s based on 19 climatic factors using Maxent model

2．2 主导气候因子及其对伯乐树地理分布区的影响

从表1 可见，基于1950—2000年气候数据，最暖季度平均雨量对伯乐树潜在适生区影响最大，贡献率为45.5%．此外，温度条件也对该种潜在适生区产生了一定影响，如昼夜温差月均值的影响居第2 位，其贡献率为10.1%．相对贡献率接近或大于5%的气候因子为温度季度化、最冷月份最低气温、昼夜温差与年温差比值、最干季度平均温度和最湿季节平均温度．相对贡献率大于10%的气候因子(主导因子)对伯乐树分布概率的影响见图4．当最暖季度平均雨量低于500 mm 时，伯乐树的分布概率随着降雨量的上升而迅速增加;当降雨量超过500 mm 时，其分布概率仍随着降雨量的增加而缓慢上升;当降雨量在2 750 mm 时，其分布概率达到平衡(图4A)．因此，最暖季度平均雨量＞2 750 mm 后将不利于伯乐树生长．图4B 显示，昼夜温差月均值在7.8～9.8 ℃范围内，伯乐树的分布概率达到最高，低于7.8 ℃或高于9.8 ℃均不利于伯乐树生长．

表2 伯乐树在中国13 个省(地区)6 个等级的分布区面积Tab．2 Distribution areas of Bretschneidera sinensis in 13 provinces(regions)of China based on 6 classes km2

图4 主导气候因子对伯乐树分布概率影响的曲线Fig．4 Response curves of distribution probability of Bretschneidera sinensis to two important climatic factors of the Maxent model

3 讨论与结论

3．1 Maxent 模型的适用性

基于本研究结果，利用气候因子和最大熵模型构建的伯乐树地理分布与气候关系模型的AUC 值高达0.993(P=0.01)，表明预测结果非常好，能够较为准确地反映该种的潜在适生区．越来越多的研究证明，Maxent 模型简单易操作，与物种的拟合度较高，能准确反映物种的潜在适生区．该模型运算已知物种的地理分布点和环境因子数据，已知物种分布点数据越多、越准确，对模型的构建越有利，精确度也越高［21-22］．本研究通过各个渠道尽可能多地获得该种的地理分布信息，包括CVH 标本数据库、各地植物志、已报道文献等，并加以实地考察，最大限度地保证数据来源全面、可靠，这是本研究具有较高AUC 值的前提之一．此外，本研究利用ARCGIS 软件进行数据处理和分析，使物种分布数据和气候数据在每栅格单元上相对应，进一步提高了数据的准确性，降低了误差．

3．2 环境因子和气候因子对伯乐树潜在适生区的影响

根据本研究结果，伯乐树最适宜分布在中海拔(500～2 000 m)山地的林中．部分群体位于耕地边、林缘，这些地区更容易遭受人为干扰，生境破坏严重，生态条件恶劣，致使伯乐树在这些地区的群体小，个体少，从而濒临危险，部分群体甚至灭绝．从各气候因子对伯乐树地理分布的影响来看，最暖季度降雨量的影响最为明显，充沛的降水和潮湿的环境有利于其生长．根据野外调查资料，该种大多生长于湿润沟谷、溪旁坡地，尤其是华南地区的一些群体，多生长于相对湿度在80%以上的环境．同时，该树种幼树耐荫，环境中的高大乔木可为其遮挡阳光．而在耕地或林缘的群体，缺乏遮光条件，极不利于伯乐树幼苗成长．另外，一些高海拔地区往往高寒、干燥，且昼夜温差较大，也不利于伯乐树的生长．因此，亚热带中海拔山地气候温暖湿润，为伯乐树生长提供了有利条件．

3．3 伯乐树潜在适生区特征及其对气候变化的响应

本研究显示，伯乐树潜在适生区位于我国亚热带常绿阔叶林、海拔500～2 000 m 的中山地带，分布区位于南岭地区以及幕阜山、罗霄山和大娄山地区．伯乐树为第三纪古热带植物区系的孑遗种，同其他古老孑遗植物类似，其地理分布区不可避免地会受到第四纪冰期气候波动的影响．通常，冰期气候剧烈波动的情况下，一些地形复杂多样的山区仍然保持相对稳定的生态环境和地质条件，从而成为一些物种冰期时的适合生境［23］．我国的西南山地和南岭山地，地形复杂多样，没有直接受到第四纪冰川的影响，富集植物物种多样性，是许多第三纪孑遗植物的冰期避难所，分布着许多第四纪气候性地理残遗种．位于我国中部地区的幕阜山和罗霄山在一定程度上阻挡了第四纪冰川的影响，也为一些物种的生存提供了适合的生境条件，独特的地形地貌和地质条件使得这些地区成为了伯乐树的潜在适生区．对其他常绿阔叶林植物进行最大熵模型的潜在适生区预测，同样表明这些地区为冰期避难所，如杜鹃、鄂报春等［24-25］．此外，有研究指出，随着CO2浓度增加，全球气候变暖，北半球植被带有向北方高纬度迁移的趋势［26-27］．本研究获得了相似的预测结果．基于21世纪70年代的气候数据，伯乐树未来的潜在适生区将由南向北迁移和扩张．同时，该种最适分布区面积将减少，说明全球气候变暖一定程度上将导致该种生境适应性降低．

3．4 伯乐树保护策略的建议

近年来，生态位模型在生物多样性保护研究中的作用和意义日益突出．在自然保护区的选择和设计过程中，构建相应的生态位模型常常具有重要的价值和意义［28-30］．根据本研究结果，我们对伯乐树群体就地保护和迁地保护工作提出几点建议．南岭地区以及幕阜山、罗霄山和大娄山地区具有较高的生物多样性，也是伯乐树在1950—2000年期间及21 世纪70年代的最适分布区，建议对这些地区的伯乐树群体进行就地保护．盲目开荒、乱砍滥伐和过度放牧是造成伯乐树原有生境破坏的重要原因．针对伯乐树最适分布区群体，建议选择这些地区建立自然保护小区，进行规范管理，避免人为干扰和破坏．在自然保护小区内，收集伯乐树种子，培育异龄幼苗．通过人工授粉，加强个体及群体间的交流，增加个体数量和群体规模．同时，对最适分布区群体的自然演替过程进行长期跟踪和监测，以保证最大限度地保护现有群体．生态位模型可以帮助我们判断未知群体的实际分布地［28-30］．广西壮族自治区虽为最适分布区，但其现有地理分布点却很少．我们建议核查该地区是否存在一些记录空缺的分布点，并对该地区进行全面的野外调查，发现未知群体，及时对这些群体进行生态学研究和种质资源保护．此外，迁地保护可以使物种占据新的分布点，是特有珍稀濒危植物保护策略中一个十分重要的环节．位于伯乐树最适分布区以外的群体，尤其是一些边缘群体，个体数量少，群体小，存在生境片断化和地理隔离的风险，建议采取迁地保护策略．根据该种生物学特性和生态学特性以及适宜的气候因子数据，结合已有的迁地保护成功实例，选择适宜地区建立规范的迁地保护基地，如南岭地区．在迁地保护区域内，建立种质资源库和种苗基地，进行长期的管理和监测，保证湿度和温度，同时进行野外回归实验，最大限度地保护伯乐树的种质资源．

致谢:感谢中国科学院华南植物园罗中莱博士和黄惠润博士在数据分析和文章撰写中提供的帮助和支持!

［1］ELITH J，GRAHAM C H，ANDERSON R P，et al．Novel methods improve prediction of species’distributions from occurrence data［J］．Ecography，2006，29(2):129-151．

［2］PHILLIPS S J，ANDERSON R P，SCHAPIRE R E．Maximum entropy modeling of species geographic distributions［J］．Ecol Model，2006，190(3/4):231-259．

［3］KUMAR S，STOHLGREN T J．Maxent modeling for predicting suitable habitat for threatened and endangered tree Canacomyrica monticola in New Caledonia［J］．J Ecol Nat Environ，2009，1(4):94-98．

［4］KREMEN C，CAMERON A，MOILANEN A，et al．Aligning conservation priorities across taxa in Madagascar with high-resolution planning tools［J］．Science，2008，5873(320):222-226．

［5］罗翀，徐卫华，周志翔，等．基于生态位模型的秦岭山系林麝生境预测［J］．生态学报，2011，31(5):1221-1229．

［6］张蓝月，叶向斌，刘念，等．白花油麻藤的地理分布及适生区预测［J］．广西植物，2012，32(1):27-32．

［7］吴征镒．中国种子植物属的分布区类型［J］．云南植物研究，1991，13(S4):1-139．

［8］吴征镒，周浙昆，李德铢，等．世界种子植物科的分布区类型系统［J］．云南植物研究，2003，25(3):245-257．

［9］中国科学院中国植物志编辑委员会．中国植物志:第34 卷［M］．北京:科学出版社，2013．

［10］中国科学院华南植物园．广东植物志［M］．广州:广东科技出版社，2009:1-4406．

［11］汪松，解焱．中国物种红色名录［M］．北京:高等教育出版社，2004:1-365．

［12］于永福．国家重点保护野生植物名录:第一批［J］．植物杂志，1999，5:4-11．

［13］傅立国．中国植物红皮书:稀有濒危植物:第1 册［M］．北京:科学出版社，1991:94-195．

［14］陈建民．贵州省珍惜濒危植物信息系统［D］．重庆:西南师范大学，2003．

［15］何克军，李意德．广东省国家I 级重点保护野生植物资源现状及保护策略［J］．热带亚热带植物学报，2005，13(6):519-525．

［16］乔琦．珍稀植物伯乐树的濒危机制和保育技术研究［D］．广州:中国科学院华南植物园，2010．

［17］张荣京，邢福武，萧丽萍，等．海南鹦哥岭的种子植物区系［J］．生物多样性，2007，15(4):382-392．

［18］VAN VUUREN D P，EDMONDS J A，KAINUMA M，et al．The representative concentration pathways:An overview［J］．Climatic Change，2011，109(1/2):5-31．

［19］FIELDING A H，BELL J F．A review of methods for the assessment of prediction errors in conservation presence /absence models［J］．Environ Conserv，1997，24(1):38-49．

［20］SWETS J A．Measuring the accuracy of diagnostic systems［J］．Science，1988，240(4857):1285-1293．

［21］LOZIER J D，ANIELLO P，HICKERSON M J．Predicting the distribution of Sasquatch in western north America:Anything goes with ecological niche modeling［J］．J Biogeogr，2009，36(9):1623-1627．

［22］曹向锋，钱国良，胡白石，等．采用生态位模型预测黄顶菊在中国的潜在适生区［J］．应用生态学报，2010，21(12):3063-3069．

［23］HEWITT G M．Some genetic consequences of ice ages，and their role，in divergence and speciation［J］．Biol J Linn Soc，1996，58(3):247-276．

［24］YAN Haifei，ZHANG Caiyun，WANG Fengying，et al．Population expanding with the phalanx model and lineages split by environmental heterogeneity:A case study of Primula obconica in subtropical China［J］．PLoS One，2012，7(9):e41315．

［25］LI Yong，YAN Haifei，GE Xuejun．Phylogeographic analysis and environmental niche modeling of widespread shrub Rhododendron simsii in China reveals multiple glacial refugia during the last glacial maximum［J］．J Syst Evol，2012，50(4):362-373．

［26］LUCHT W，PRENTICE I C，MYNENTR B，et al．Climatic control of the high-latitude vegetation greening trend and pinatubo effect［J］．Science，2002，296(5573):1687-1689．

［27］冯彩云，许新桥，马月萍，等．植物对气候变化生理反应研究进展［J］．世界林业研究，2012，25(4):16-21．

［28］朱耿平，刘国卿，卜文俊，等．生态位模型的基本原理及其在生物多样性保护中的应用［J］．生物多样性，2013，21(1):90-98．

［29］BOURG N A，MCSHEA W J，GILL D E．Putting a CART before the search:Successful habitat prediction for a rare forest herb［J］．Ecology，2005，86(10)，2793-2804．

［30］MENON S，CHOUDHURY B I，KHAN M L，et al．Ecological niche modeling and local knowledge predict new populations of Gymnocladus assamicus，a critically endangered tree species［J］．Endanger Species Res，2010，11(2):175-181．