氣候變化下花櫚木在我國的潛在分布區(qū)預測

宋 穎,張港隆,賈全全,溫 強,莫曉勇,劉麗婷*

(1.江西省林業(yè)科學院，江西 南昌 330032；2.華南農(nóng)業(yè)大學 林學與風景園林學院，廣東 廣州 510642)

氣候變化是當前和未來影響生物多樣性變化的主要因素之一[1]。21世紀以來，氣候變化帶來的全球氣候變暖和極端天氣頻發(fā)，加之水文土壤及人類活動等多種因素的共同影響，物種生存環(huán)境遭受破壞，生物多樣性喪失速率加快，大量珍稀瀕危物種的生存與發(fā)展受到威脅。而在未來全球氣溫及降水的持續(xù)變化背景下，生態(tài)系統(tǒng)和生物類群在區(qū)域尺度的地理分布也必將受到影響?，F(xiàn)有研究表明，自末次盛冰期以來，劇烈的氣候變化曾導致諸多物種地理分布動態(tài)變化，生境破碎化加劇，野生種群遺傳多樣性喪失[2-3]，顯然，氣候?qū)ξ锓N的生存狀態(tài)和潛在分布具有深遠的影響，是大尺度上決定物種分布格局的關(guān)鍵因素[4]?，F(xiàn)如今人們也越來越認識到，由于社會、經(jīng)濟、環(huán)境等多種驅(qū)動因素之間的相互作用，氣候變化、生物多樣性喪失和環(huán)境污染等諸多挑戰(zhàn)并存[5]。在此背景下，一種新的社會經(jīng)濟發(fā)展情景——共享社會經(jīng)濟路徑(shared socio-economic pathways,SSPs)被逐漸用于氣候變化研究和評估。SSPs于2010年開發(fā)，以代表濃度路徑(representative concentration pathways，RCPs)為基礎(chǔ)，充分考慮相關(guān)的社會經(jīng)濟因素，可用于預測2100年不同氣候政策下溫室氣體的排放情景；其包含了4種經(jīng)濟發(fā)展路徑：SSP1——可持續(xù)發(fā)展路徑、SSP2——中間發(fā)展路徑、SSP3——區(qū)域競爭發(fā)展路徑、SSP4——不均衡發(fā)展路徑[6]、SSP5——以化石燃料為主發(fā)展路徑。SSPs在全球和區(qū)域尺度已廣泛應用于模擬氣候變化對人口與經(jīng)濟發(fā)展、氣象災害、農(nóng)業(yè)生產(chǎn)、生態(tài)系統(tǒng)等諸多領(lǐng)域的可能影響[7-8]。

花櫚木(Ormosiahenryi)為豆科紅豆屬多年生常綠喬木，在我國多分布于熱帶、亞熱帶地區(qū)的低海拔常綠闊葉林中?；澳緲渥送Π危募境Ｇ?，可作為優(yōu)良的城市綠化和庭院美化樹種。其木材致密堅硬，質(zhì)地美觀，是制作高檔家具和工藝品的原料。在藥用價值方面，其根、莖、葉和樹皮均可入藥，具有多種藥用功效[7]。有研究顯示，含黃酮C-糖苷的花櫚木葉乙醇提取物具有抗抑郁作用[8]，其精油也被證明可以減輕抑郁行為[9]，未來，花櫚木還將有望成為緩解慢性應激所致記憶障礙的潛在藥物片。作為我國珍貴樹種，花櫚木具有巨大的經(jīng)濟和生態(tài)潛力。然而，隨著社會經(jīng)濟發(fā)展的需要，由于人類對其自然資源的不合理開發(fā)利用，花櫚木的種群生存與發(fā)展已受到嚴重威脅，其適生分布區(qū)域亦呈逐漸縮小且破碎化的趨勢，另外，因其種子種皮致密堅硬，不易萌芽，自然繁殖能力較弱，植株的生長發(fā)育對生態(tài)環(huán)境要求較高，目前花櫚木的野外種群數(shù)量已經(jīng)銳減[7]，現(xiàn)已被列為國家Ⅱ級重點保護野生植物。迄今，國內(nèi)外學者已經(jīng)在種苗繁育技術(shù)、藥用價值開發(fā)、群落資源分布等領(lǐng)域?qū)澳菊归_了深入細致的研究，但還未見基于生態(tài)位模型對其進行潛在地理分布預測的相關(guān)報導。在此背景下，本研究基于生態(tài)位模型研究氣候變化對花櫚木適生范圍的影響，對有效開展該物種多樣性保護和資源培育利用具有重要的現(xiàn)實意義。

目前，在預測、解析物種分布范圍及其變化情況等領(lǐng)域被運用較多的生態(tài)位模型有：生物氣候分析和預測系統(tǒng)模型(BIOCLIM)、生態(tài)位因子分析模型(ecological niche factor analysis,ENFA)、CLIMEX地點比較模型、基于規(guī)則集的遺傳算法模型(genetic algorithm for rule-set production,GARP)和最大熵模型(maximum entropy model,MaxEnt)等。其中，最大熵模型在已知物種分布樣本點和相關(guān)環(huán)境變量的基礎(chǔ)上，建立約束條件，通過找出物種分布規(guī)律的最大熵而對其潛在分布作出模擬，相比于其他模型，MaxEnt具有建模直觀、預測精度高、運算速度快等優(yōu)點[10-11]。目前MaxEnt模型已被廣泛用于入侵生物潛在分布模擬、物種的適生環(huán)境預測、全球氣候變化對物種分布的影響評估等方面，在生態(tài)學、保護生物學、生物地理學等多個領(lǐng)域得到開拓和應用[12-14]。

本研究嘗試在可持續(xù)發(fā)展(SSP1)、中度發(fā)展(SSP2)和常規(guī)發(fā)展(SSP5)3種氣候模式情景下，應用MaxEnt模型和ArcGIS軟件，對花櫚木在我國的潛在分布區(qū)與格局變化進行預測，以期對花櫚木的種質(zhì)資源保護與適地栽培利用提供可靠參考。

1 材料與方法

1.1 MaxEnt模型及地圖數(shù)據(jù)獲取

本研究使用MaxEnt3.4.1(https://biodiversityinformatics.amnh.org/open_source/maxent/)進行模型的構(gòu)建與驗證，采用地理信息系統(tǒng)軟件ArcGIS10.2進行數(shù)據(jù)的處理與結(jié)果的可視化，地圖資源獲取自國家基礎(chǔ)地理信息系統(tǒng)(http://nfgis.nsdi.gov.cn/)。

1.2 花櫚木分布信息數(shù)據(jù)獲取及處理

通過查閱地方植物志、野外實地采集樣點，以及在線查閱全球生物多樣性信息機構(gòu)(GBIF,https://www.gbif.org/)、中國數(shù)字植物標本館(CVH，http://www.cvh.ac.cn/)、教學標本資源共享平臺(http://mnh.scu.edu.cn/main.aspx)、中國自然保護區(qū)標本資源共享平臺(http://www.papc.cn/)等信息平臺，初步得到花櫚木在國內(nèi)的599條樣點信息?；诎俣仁叭∽鴺讼到y(tǒng)(http://api.map.baidu.com/lbsapi/getpoint/)拾取補充缺少的坐標信息[15]，去除重復地域、重復坐標及不確定的樣本信息，篩選后最終得到樣點188個(圖1)。

圖1 花櫚木樣點分布

1.3 環(huán)境數(shù)據(jù)獲取及處理

本研究所采用的初始環(huán)境變量共20個，其中，當前(2000s-2020s)的19個生物氣候數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/)的中國地面氣候資料月值數(shù)據(jù)集(1951s-2020s)，由ANUSPLIN 4.36軟件中的薄板光滑樣條函數(shù)(Thin Plate Smoothing Spline Method)(樣條數(shù)為4、高程數(shù)據(jù)dem為協(xié)變量)插值而成[16]。未來(2050s，2070s)19個生物氣候數(shù)據(jù)以及1個高程數(shù)據(jù)均來源于世界氣候數(shù)據(jù)庫(/https://worldclim.org/)，坐標系為WGS84，空間分辨率為5 km；未來氣候情景采用共享社會經(jīng)濟路徑(SSPs)，根據(jù)氣候情景與社會經(jīng)濟情景的對應情況，SSP1-2.6、SSP2-4.5、SSP5-8.5等3種氣候情景模式分別對應RCP2.6(低排放)、RCP4.5(中等排放)和RCP8.5(高排放)情景。

1.4 氣候變量的選擇

由于環(huán)境變量間存在一定的相關(guān)性，為避免其多重共線性導致模型過度擬合，需排除環(huán)境變量間的空間自相關(guān)性后方可用于模型計算[17]。首先對20個環(huán)境變量使用MaxEnt進行預建模，得到所有環(huán)境變量的貢獻率；其次，使用ArcGIS對20個環(huán)境變量進行空間自相關(guān)性分析，得到相關(guān)系數(shù)矩陣。當相關(guān)系數(shù)<0.85時，對應的氣候變量得以保留；若相關(guān)系數(shù)≥0.85，則根據(jù)MaxEnt預建模得到的貢獻率，保留貢獻率較大且具有明確生物學意義的環(huán)境變量，最終篩選得到9個環(huán)境變量用于MaxEnt模型運行，包括等溫性(bio03)、溫度季節(jié)性變化系數(shù)(bio04)、最濕季度平均溫度(bio08)、最暖季度平均溫度(bio10)、最冷季度平均溫度(bio11)、降水季節(jié)變化系數(shù)(bio15)、最濕季度降水量(bio16)、最干季度降水量(bio17)和高程(elev)。

1.5 MaxEnt模型構(gòu)建與優(yōu)化

有研究表明，MaxEnt模型的預測結(jié)果與特征組合(feature combination，F(xiàn)C)、正則化乘數(shù)(regularization multiplier，RM)密切相關(guān)，使用默認參數(shù)構(gòu)建MaxEnt模型時，易出現(xiàn)模型過度擬合及轉(zhuǎn)移能力不一的情況[18]。因此本研究通過在R 4.03軟件中調(diào)用ENMeval數(shù)據(jù)包來調(diào)整RM和FC的參數(shù)值，進而優(yōu)化MaxEnt模型，得到最優(yōu)模型[19]。

在MaxEnt軟件中，共有線性(Linear，L)、二次型(Quadratic，Q)、片段化(Hinge，H)、乘積型(Product，P)和閾值型(Threshold，T)5種特征可供選擇，默認參數(shù)為RM=1，F(xiàn)C=LQHP。為優(yōu)化模型，本研究將RM設(shè)置為0.5～4，每次增加0.5來進行測試，共8個梯度[20]；采用6個特征組合：L、LQ、H、LQH、LQHP、LQHPT；運用ENMeval數(shù)據(jù)包分析上述48種參數(shù)組合對模型構(gòu)建的影響，最終根據(jù)Akaike信息量準則(Akaike information criterion correaction，AICc)、訓練AUC和測試AUC值之差(the difference between training and testing AUC，AUCDIFF)和10%訓練遺漏率(10% training omission rate，OR10)三者來評估模型性能[21]。其中，AICc是衡量統(tǒng)計模型擬合優(yōu)良性的一種標準，可以權(quán)衡模型的復雜度和擬合數(shù)據(jù)的優(yōu)良性，其優(yōu)先考慮AICc值最小(即delta AICc=0)的模型；AUCDIFF和OR10均用來衡量模型的過度擬合程度，前者與模型的過擬合程度呈正相關(guān)，后者如果大于期望值0.1，則表明模型存在過度擬合[21]。

1.6 MaxEnt模型評價

模型優(yōu)化后，將花櫚木的分布點數(shù)據(jù)以及篩選得到的9個環(huán)境變量導入MaxEnt軟件，選擇刀切法(Do jackknife to measure variable importance)功能以檢驗各環(huán)境因子對建模的貢獻率[22]，選取創(chuàng)建響應曲線(Create response curve)功能以獲得分布概率與氣候因子之間的關(guān)系。隨機選擇25%的樣本分布點作為模型測試數(shù)據(jù)，75%作為模型訓練數(shù)據(jù)，選擇bootstrap運算方法[23]重復運行模型15次，RM與FC取優(yōu)化后的參數(shù)，其他參數(shù)為模型默認設(shè)置。以受試者工作特征曲線(Receiver Operating Characteristic Curve,ROC)下的面積值(AUC 值)大小評價模型的準確性。AUC值的取值范圍為0～1，<0.6表明預測結(jié)果失敗，0.7～0.9表示模型預測精度中等，0.9以上表示模型具有較高的精確度[24]。

1.7 花櫚木生境適宜性等級劃分

將MaxEnt輸出的分布預測數(shù)據(jù)導入ArcGIS，轉(zhuǎn)換為柵格(Raster)格式后進行重分類，采用自然間斷點分級法(Natural Breaks Classification Method)將花櫚木適生區(qū)劃分為5個等級，分別為非適生區(qū)(0～0.06)、邊緣適生區(qū)(0.06～0.21)、低適生區(qū)(0.21～0.38)、中適生區(qū)(0.38～0.53)、高適生區(qū)(0.53～0.71)，并計算、評估未來生態(tài)適宜區(qū)面積變化。

2 結(jié)果與分析

2.1 模型優(yōu)化及預測準確性檢驗

基于188條花櫚木分布記錄和9個環(huán)境變量，利用MaxEnt模型模擬預測花櫚木當前潛在分布區(qū)。當優(yōu)化參數(shù)RM=2，F(xiàn)C=LQHPT時，delta AICc=0；AUCDIFF比默認參數(shù)時降低了37.5%，表明該優(yōu)化參數(shù)下模型性能較好，過擬合程度較低，可用于模型的轉(zhuǎn)移。因此本研究采用優(yōu)化參數(shù)RM=2，F(xiàn)C=LQHPT為模型構(gòu)建參數(shù)。

表1 MaxEnt模型優(yōu)化評價指標

在優(yōu)化參數(shù)下進行模型預測，采用受試者工作特征曲線分析法進行模型精度檢驗。如圖2所示，訓練集的平均AUC值和平均標準偏差分別為0.954和0.004，表示由此模型預測得出的花櫚木生態(tài)適生區(qū)具有較高的可信度。

圖2 MaxEnt模型中ROC預測結(jié)果

2.2 影響花櫚木地理分布的主導氣候因子

用模型檢測生物氣候變量對于花櫚木分布增益的貢獻率(表2)，采用刀切法評估各生態(tài)因子對花櫚木地理分布的重要性(圖3)，發(fā)現(xiàn)bio17最干季度降水量、bio16最濕季度降水量、bio11最冷季度平均溫度這3個變量對花櫚木的影響最大?？梢姡邓亢蜏囟仁怯绊懟澳镜乩矸植嫉闹鲗夂蛞蜃?。

表2 MaxEnt模型中各個環(huán)境變量的重要性

圖3 環(huán)境因子對花櫚木分布作用的刀切法評價

為進一步探究花櫚木在當前潛在適生區(qū)的氣候特征，將上述3個主導因子分別導入MaxEnt模型進行單因子建模，繪制單因子響應曲線，當花櫚木分布概率P≥0.53，適宜性等級為高適宜分布時，得到各個主導氣候因子的適宜范圍。結(jié)果顯示，bio17最干季度降水量、bio16最濕季度降水量、bio11最冷季度平均溫度的閾值分別為155～252 mm、636～1 723 mm、7.6～13.9℃。

2.3 花櫚木在當前氣候下的潛在適生分布區(qū)

將花櫚木分布點和環(huán)境變量數(shù)據(jù)導入MaxEnt模型進行模擬，得到花櫚木在當前氣候下的潛在分布區(qū)(圖4)，主要分布區(qū)范圍18°25′11″-32°26′37″N，96°36′32″-122°52′39″E，總適生面積約195.67萬km2，已覆蓋現(xiàn)存分布區(qū)域，表明模擬得到的花櫚木潛在分布區(qū)與實際分布范圍較為一致。利用 Reclassify 重分類工具進行等級劃分，得到的高適生區(qū)、中適生區(qū)、低適生區(qū)及邊緣適生區(qū)分別占總適生區(qū)總面積的 20.55%、26.46%、17.44%和35.55%。由圖4可知，花櫚木的適生范圍主要集中在我國東南地區(qū)，其中，高適生區(qū)主要位于江西、湖南、福建、廣西、浙江等地，面積約為39.38萬km2。

圖4 當前氣候條件下花櫚木在我國的潛在分布區(qū)

2.4 花櫚木在未來氣候下的潛在分布區(qū)變遷

由模型預測結(jié)果可知(表3)，與當前氣候相比，未來不同氣候情景下，花櫚木當前的適生區(qū)絕大部分仍然保留，但相較于當前分布面積總體上縮減。SSP1-2.6氣候情景下，2050s花櫚木的潛在分布區(qū)面積較當前減少4.10萬km2；在2070s，中適生區(qū)和邊緣適生區(qū)面積均縮減，此時總分布區(qū)面積為188萬km2，較當前面積變化其縮減比例為3.9%。SSP2-4.5氣候情景下，2050s的分布面積相較于當前有所縮減，其中，高、低適生區(qū)分別擴張，而中、邊緣適生區(qū)面積稍有下降。SSP5-8.5氣候情景下，2050s高、中及邊緣適生區(qū)均縮減，低適生區(qū)擴張；2070s相較于當前適生區(qū)面積變化較小，但高適生區(qū)和低適生區(qū)面積均有所增加。

表3 未來不同氣候情景下花櫚木的潛在適生分布區(qū)面積(×104 km2)

將不同氣候情景下花櫚木適生區(qū)的空間格局進行對比分析(圖5)，結(jié)果顯示，相較于當前，除SSP5-8.5(2050s)外，未來各氣候情景下花櫚木的高適生區(qū)均呈現(xiàn)擴張狀態(tài)，主要集中在安徽西北部、福建西南部、湖南中南部、湖北東南部、貴州東部等地區(qū)，臺灣北部也有零星的增加，廣東自沿海地區(qū)逐漸向內(nèi)陸縮減。在未來，花櫚木的核心適宜區(qū)于東西向穩(wěn)定，于南北向中部縮減，尤以廣東中部及江西中部收縮顯著，這些遷移的核心適生區(qū)為花櫚木應對氣候變化的敏感區(qū)，需引起重視。

圖5 不同氣候情景下花櫚木潛在適生區(qū)分布

結(jié)合適生區(qū)空間變化趨勢圖(圖6)，統(tǒng)計不同氣候變化情景下花櫚木的適生空間格局變化，得知花櫚木未來總體分布范圍相對穩(wěn)定，潛在適生區(qū)面積波動較小。新增區(qū)域主要位于湖北北部、貴州西部、云南西北部以及臺灣等地，此外，江蘇中東部及等地亦有零星拓展區(qū)域，于SSP5-8.5(2050s)氣候情景下擴張最為明顯；收縮部分主要為邊緣適生區(qū)，于云南南部、安徽西北部與江蘇中部等地，呈狹長帶狀。

圖6 不同氣候情景下花櫚木適生區(qū)的空間格局變化

3 結(jié)論與討論

本研究采用政府間氣候變化專門委員會(IPCC)提出的共享社會經(jīng)濟路徑(SSPs)，研究在未來社會經(jīng)濟發(fā)展和政策干預下，溫室氣體排放的情景變化對物種的地理分布的影響。相較于以往的模式比較計劃，引入SSPs的CMIP6能夠在人口、經(jīng)濟和技術(shù)等社會發(fā)展指標的變化的情境下，實現(xiàn)對各種減排措施的成本效益評估，同時提供了比RCPs更多元化的空氣污染物排放情景，這些改進很大程度上彌補了CMIP5中RCPs情景對社會經(jīng)濟活動不確定性[25]，從而能夠提高氣候研究對政策制定的支持力。

本研究中，MaxEnt模型檢驗AUC值達到了0.954，表明該模型預測結(jié)果具有極高的精度和可信度。花櫚木潛在適生境集中分布于我國東南部地區(qū)，其中，核心分布區(qū)主要位于江西、湖南、福建等省區(qū)，預測結(jié)果與當前實際分布相一致。花櫚木實際分布區(qū)地處大陸性中亞熱帶季風濕潤氣候，氣溫的季節(jié)變化顯著，夏季高溫多雨，持續(xù)時間長，冬季寒冷期短，最冷月平均氣溫在0℃以上。在野外調(diào)查中發(fā)現(xiàn)，花櫚木常分布于100～1 300 m的山坡、溪谷及雜木林中，多喜于在光熱充足、濕度適宜的亞熱帶常綠闊葉林生境下生存，說明模型預測結(jié)果與花櫚木的實際生長特性相吻合。

有研究表明，1960-2013年中國各地年際氣溫呈上升趨勢，年際降水除西南部趨于減少，其余大部分地區(qū)均增加，且未來還將持續(xù)這一趨勢[26]，另外，我國降水類型正不斷趨向于極端降水[68]，而氣溫、光照和水分是植物物候期變化最主要的影響因子[27]。因此，考慮到近20 a來氣候變化的顯著性，本研究采用2000-2020年中國氣象數(shù)據(jù)網(wǎng)的數(shù)據(jù)插值到研究區(qū)域并進行物種潛在分布的模擬預測。目前，MaxEnt模型在物種分布預測的研究中應用較多，有研究表明，MaxEnt模型的預測結(jié)果主要與正則化乘數(shù)(RM)、特征組合(FC)、最大背景點數(shù)(BG)密切相關(guān)；而最大背景點數(shù)，即在模型構(gòu)建區(qū)域選取的不存在分布點，MaxEnt模型會將其關(guān)聯(lián)的環(huán)境變量用于和物種存在分布點所關(guān)聯(lián)的環(huán)境變量進行對比，分析物種生態(tài)需求，進而預測物種分布[28]。本研究主要探索了RM與FC對模型的影響，雖有研究表明BG的選擇對物種的潛在分布和生態(tài)位的影響較小[29]，但下一步可考慮BG的影響，以求進一步優(yōu)化模型。另外，本研究著重探討氣候因子對物種分布的影響，未考慮到其他諸如土壤、人類干擾強度、種間相互競爭、物種擴散能力等因素，未來可充分評估各變量的綜合作用，以進一步提升模型預測效果。

氣候環(huán)境是決定物種分布的關(guān)鍵因素，物種的分布變化也是響應氣候變化最明顯的直接反應[30]。本研究中，運用MaxEnt模型與ArcGIS對花櫚木不同氣候情景下各適生等級分布區(qū)面積進行了預測，并對總體分布范圍的新增、縮減及穩(wěn)定區(qū)域?qū)崿F(xiàn)可視化表達，從而能夠更直觀了解到未來花櫚木適生區(qū)的擴張、衰退結(jié)果。值得注意的是，通過對不同氣候變化情景下適生區(qū)的空間格局進行對比分析，得知花櫚木在未來的總體空間分布與當前無太大差異，說明其對未來的氣候變化有一定的適應能力。然而，雖然目前花櫚木的適生區(qū)范圍尚屬廣泛，但通過對花櫚木野生種群資源調(diào)查發(fā)現(xiàn)，由于資源被過度開發(fā)，又因其自然更新繁衍能力較弱等內(nèi)外因素影響，花櫚木野生種群規(guī)模已經(jīng)銳減，對其制定相關(guān)保護措施仍十分重要。

在對花櫚木種質(zhì)資源保護工作中，可在中適生區(qū)適當增加花櫚木的遷入及擴散，選擇性地營造花櫚木人工群落，使之有利于花櫚木種群繁衍、擴散。對于高適生區(qū)，應就地保護與人工栽培并行，鑒于花櫚木的珍稀瀕危性，及其重要的經(jīng)濟和生態(tài)價值，需加強公眾對珍稀瀕危植物的了解，禁止濫采亂伐等行為?？梢暦N群密度大小建立自然生態(tài)保護區(qū)，加強對敏感區(qū)域的重點觀測，必要時采取補救措施，以維持其種群豐富度及遺傳多樣性。此外，在遷地保護及引種栽培過程中，除了要遵循適地適樹的原則外，還應綜合考慮病蟲害、土壤質(zhì)地、種間競爭等因素，以提高花櫚木栽培成活率。