基于組合模型預測蛇足石杉在中國的適生區

陳 勤，孔夢盈，孫奇芳，管畢財

(南昌大學 生命科學學院，南昌 330031)

氣候是物種分布的決定因素之一[1]。氣候變化會改變棲息環境，降低生物多樣性，增加瀕危植物滅絕的風險[2]。21世紀以來全球氣候變化復雜多樣，物種對氣候的適宜性隨之發生變化，從而影響物種分布格局，進而威脅物種的生存[3]，迫切需求將氣候變化納入物種保護的研究中[4]。

物種分布模型(species distribution models, SDMs)是基于生態位理論量化環境與物種之間關系的相關性模型[5]。它能夠揭示物種對棲息地的適宜性與環境的關系，被廣泛用于生態學研究中，特別是在預測物種的分布及其熱點[6]、管理入侵物種[7]、保護瀕危物種[8]及驗證氣候變化與物種分布的關系[9]等方面。用于構建物種分布模型的算法主要有：人工神經網絡(artificial neural networks, ANN[10])、生態位因子分析(ecological niche factor analysis, ENFA[11])、廣義線性模型(generalized linear models, GLM[12])及最大熵模型(maximum entropy, MaxEnt[13])等，但是基于單一算法構建的模型準確性有待提升[14]。基于R程序包Biomod2[15]構建組合模型，能夠提升預測結果的可靠性，日益廣泛應用于模型構建。Ren等[16]使用Biomod2調查了對中國云南亞熱帶常綠闊葉林的物種分布及生物多樣性。Ray等[17]使用Biomod2預測了橡樹林地的分布。Vacchiano等[18]使用Biomod2預測了意大利西北部入侵植物窄葉黃菀的入侵情況與潛在分布。Legal等[19]使用Biomod2研究了分布有限且不均勻的孑遺生物短角寶鳳蝶的潛在分布及其保護，這些研究均得到了比較好的預測效果。

蛇足石杉[Huperziaserrata(Thunb. ex Murray) Trevis.]屬于石松科(Lycopodiaceae)石杉屬(Huperzia)多年生蕨類植物[20]。蛇足石杉喜陰濕[21]，主要生長在海拔350～1 700 m的林下、灌叢和路旁[22]，廣泛分布于亞洲、大洋洲、美洲熱帶、亞熱帶及溫帶地區，在中國大部分地區都有分布但生物量普遍較少[23]。蛇足石杉具有重要的藥用價值，從蛇足石杉中分離出的蛇足石杉堿及石杉堿甲等生物活性物質可用于治療阿茲海默癥、重癥肌無力和記憶力下降等臨床癥狀[24-26]。但目前蛇足石杉產業化繁育尚未取得成功，材料獲取主要依賴野生資源，由于過度采集和生態環境惡化等因素影響，生物資源驟減[27]。因此，對蛇足石杉生物資源的保護刻不容緩，2021年9月，蛇足石杉被國家林草局正式列為國家二級重點保護野生植物[28]。迄今為止，對蛇足石杉的研究也主要集中在化學成分[29]、藥理作用[30]、快速繁殖[31]及分子生物學[32]等方面，為了開展對該物種的保護工作，我們需要調查收集蛇足石杉的分布情況并進行適生區的預測。因此，本研究基于蛇足石杉在中國的分布數據和氣候數據構建組合模型，預測蛇足石杉當前和未來氣候情景下在中國的適生區并計算適生區面積及質心變化，為蛇足石杉及其生物資源的保護提供理論支持。

1 材料和方法

1.1 蛇足石杉分布數據

通過中國數字植物標本館(http://www.cvh.ac.cn/)、中國國家標本資源平臺(http://www.nsii.org.cn/)、全球生物多樣性信息網絡(https://www.gbif.org/)和文獻查閱獲取蛇足石杉分布數據，對有準確地址和經緯度的數據進行校對，對有準確地址無經緯度的數據進行經緯度定位，剔除重復以及過于密集的分布數據，最終總共得到186條蛇足石杉分布數據，從地理國情監測云平臺(http://www.dsac.cn/)下載國家行政區矢量地圖作為底圖進行可視化表達(圖1)。

圖1 蛇足石杉在中國的分布

1.2 氣候數據的獲取與篩選

從世界氣候數據庫(WorldClim, https://www.worldclim.org/)下載20世紀中、后期(1960～1990)的氣候數據作為當前氣候數據和對中國氣候模擬能力較好的CCSM4模式[33]下的21世紀下半葉(2070)的RCP2.6、RCP4.5、RCP6.0及RCP8.5這4種典型濃度路徑(representative concentration pathways, RCPs)的氣候數據作為未來氣候數據，氣候數據分辨率均為2.5 minutes。RCPs表示2100年相對于1750年的輻射強迫總和，向大氣排放的溫室氣體和氣溶膠會影響地球輻射平衡產生輻射強迫[34]。RCP2.6輻射強迫為2.6 W·m-2，排放處于極低水平，到2100年地球表面溫度將上升0.3～1.7 ℃；RCP8.5為8.5 W·m-2，排放處于高水平，到2100年地球表面溫度將上升2.6～4.8 ℃；RCP4.5和RCP6.0分別為4.5 W·m-2和6.0 W·m-2，排放水平介于二者之間，到2100年地球表面溫度將分別上升1.1～2.6 ℃和1.4～3.1 ℃[35]。

氣候數據包括19個生物氣候變量(表1)，使用ArcGIS 10.6提取中國生物氣候變量數據，再用SDM Toolbox[36]移除相關性大于0.7的生物氣候變量，避免生物氣候變量之間相關性過高影響適生區預測準確性，將等溫性(Bio3)、溫度季節性(Bio4)、最暖月最高溫(Bio5)、最冷季均溫(Bio11)、最濕月降水量(Bio13)和降水季節性(Bio15)6個變量用于適生區預測。

表1 19個生物氣候變量

1.3 物種分布模型的構建與評估

基于Biomod2構建物種分布模型需要物種分布數據和不分布數據，不分布數據很難獲得[11]，通常生成假缺失數據作為不分布數據[5]。本研究通過Biomod2隨機生成1 000個假缺失數據作為不分布數據，使用Biomod2包含的10種算法(表2)分別構建物種分布模型。采用接收者操作特征曲線(receiver operating characteristic curve, ROC)和真實技巧統計法(true skill statistic, TSS)對模型進行評估。AUC(area under the curve)是ROC與橫坐標圍成的面積，值介于0.5～1之間[37]。TSS不依賴發生率，能很好評價預測結果[38]。使用TSS值大于0.8的算法用于構建物種分布組合模型，將AUC值作為組合模型的最終評價標準，其值越大說明預測效果越好。在模型構建過程中使用Biomod2內置評價標準對用于構建物種分布模型的6個生物氣候變量進行重要性評估。將模型預測結果輸入ERDAS IMAGINE 2015導出為TIFF格式，再使用ArcGIS 10.6對TIFF文件進行地理配準及可視化表達。

表2 10種用于模擬的算法

1.4 物種適生區及其變化預測

使用ArcGIS 10.6對組合模型預測結果進行歸一化，將預測結果轉化為0-1的形式，采用自然斷點分級法(natural breaks)進行分類，將預測結果分為非適生區和低、中、高適生區，再對分類結果進行重分類，采用兩種分類方法探究其適生區變化情況：1.將分類結果分為適生區和非適生區，適生區包括低、中、高適生區，即蛇足石杉的總適生區；2.將非適生區和低適生區合為一類，中適生區和高適生區合為一類。使用SDM Toolbox計算不同時期各適生區面積變化和總適生區質心變化情況，通過ArcGIS 10.6對結果進行可視化表達。

2 結果與分析

2.1 模型預測精度評估

基于AUC值和TSS值評估單一模型和組合模型的預測精度(表3)，除SRE算法構建的模型外，所有模型均AUC值大于0.8且TSS值大于0.7，說明單一算法構建的模型能較好預測蛇足石杉適生區，但組合模型AUC值和TSS值大于所有單一模型，預測精度更高。組合模型AUC值大于0.9，說明蛇足石杉適生區預測結果極好。

表3 組合模型及單一模型AUC值與TSS值

2.2 生物氣候變量重要性評估

由生物氣候變量重要性評估結果(圖2)可知，最濕月降水量(Bio13)是對蛇足石杉適生區影響最大的生物氣候變量，此外重要性由高到低依次是溫度季節性(Bio4)、最冷季均溫(Bio11)、最暖月最高溫(Bio5)、降水季節性(Bio15)和等溫性(Bio3)。總體上看，降水和溫度都對蛇足石杉適生范圍有重要影響。

Bio3.等溫性;Bio4.濕度季節性;Bio5.最暖月最高溫;Bio11.最冷季均溫;Bio13.最溫月降水量;Bio15.降水季節性

2.3 蛇足石杉適生區預測

蛇足石杉當前適生區主要位于中國南部(圖3)。從適宜性來看，蛇足石杉在貴州、重慶、湖南、福建及浙江等地區適宜性較高，在新疆、青海、甘肅、黑龍江及寧夏等地區適宜性很低。從適生區來看(表4)，當前總適生區面積226.99×104km2，其中：低適生區占33.69%，主要位于四川東部、廣西、湖南東部、廣東東部及江西等地區；中適生區占32.97%，主要位于云南和四川東部，在其他省區少量散布；高適生區占33.34%，主要位于貴州、重慶、湖南西部、福建及浙江等地區。4種未來氣候情景下適生區預測結果如圖4，低適生區持續擴張，在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下分別擴張31.47%、50.17%、55.01%和84.14%。中適生區和高適生區不斷收縮，中適生區在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下分別收縮6.03%、13.04%、15.55%和31.90%，高適生區在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下分別收縮17.54%、28.65%、43.23%和58.69%。總適生區變化不大，在RCP2.6及RCP4.5情景下分別擴張2.77%和3.06%，在RCP6.0及RCP8.5情景下分別收縮1.00%和1.73%。

表4 當前及未來(2070)氣候情景下的蛇足石杉適生區面積

圖3 當前蛇足石杉適生區

圖4 未來(2070)氣候情景下的蛇足石杉適生區

2.4 蛇足石杉適生區及質心變化

蛇足石杉總適生區面積凈變化不大，總變化與排放水平呈正比(圖5，表5)。在RCP2.6及RCP4.5情景下，蛇足石杉適生區面積分別擴張5.80%和8.46%，收縮3.32%和5.90%，凈擴張。在RCP6.0及RCP8.5情景下，蛇足石杉適生區面積分別擴張7.33%和10.69%，收縮8.82%和12.96%，凈收縮。在適生區擴張方面，4種情景都向云南、四川、遼寧及吉林等地區擴張；在適生區收縮方面，4種情景下江蘇和湖北等地區都存在明顯收縮情況，此外，RCP6.0及RCP8.5情景下湖南及安徽等地區適生區也明顯收縮。

圖5 未來(2070)氣候情景下蛇足石杉總適生區變化

表5 未來(2070)氣候情景下蛇足石杉總適生區面積變化

蛇足石杉中、高適生區面積隨排放水平由低到高持續收縮(圖6，表6)。在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下分別凈收縮14.34%、23.31%、31.63%和47.08%，四川東部、廣東、廣西、湖南及福建等地區的中、高適生區向低適生區和非適生區轉變，四川中部和云南北部等地區的少量低適生區和非適生區轉變為中、高適生區。

表6 未來(2070)氣候情景下蛇足石杉中、高適生區面積變化

圖6 未來(2070)氣候情景下蛇足石杉中、高適生區變化

蛇足石杉總適生區質心和中、高適生區質心均在湖南省西部，在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下質心均向西北轉移，與適生區變化情況一致(圖7)。

A.蛇足石杉總適生區質心變化；B.蛇足石杉中、高適生區質心變化

3 討 論

3.1 蛇足石杉的未來分布

蛇足石杉在中國分布廣泛，但準確的分布記錄極不均勻，在山東、江蘇及湖北等地區的精確分布數據尤其少甚至沒有，物種分布模型的準確性與物種分布數據的數量和精確程度密切相關[39]，本研究篩選精確到村級單位和有明確經緯度及地址的蛇足石杉分布數據，進行經緯度查詢、校對和空間過濾后用于構建模型，盡可能確保準確性，組合模型AUC值均大于0.9且與實際分布相符也說明預測結果可信度高, 但是分布數據有待未來通過完善標本記錄和準確目擊記錄進一步優化。此外，物種與環境之間的關系存在復雜性和不確定性，對氣候變化的響應存在滯后性[40]，增加了物種分布模型預測的難度。因此，物種分布模型預測的是適合物種生長的區域，即物種的潛在分布區域。本研究主要基于對19個生物氣候變量去相關性獲得的6個變量構建物種分布模型，而在物種實際生存中，海拔、坡度及地質等環境變量都可能會影響分布[41]，這些變量在深入研究中需要納入考量。

在二氧化碳與氣溶膠排放水平逐漸升高的RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下，蛇足石杉總適生區面積變化不大，但蛇足石杉中、高適生區劇烈收縮，收縮面積與排放水平成正比，說明高排放對蛇足石杉中、高適生區有消極影響，迫切需要合理調控排放保護蛇足石杉適生區。在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下，蛇足石杉適生區質心均向高緯度轉移，與當前大部分研究結果一致，隨著溫度升高物種適生區會向高緯度地區遷移[42-46]，與Shafer等[47]研究發現氣候變化會使得物種向各個方向遷移相符。

3.2 蛇足石杉的保護

蛇足石杉因含有藥用成分石杉堿甲，野生資源被大量采挖作為藥用材料，現在中國將蛇足石杉列為國家二級重點保護野生植物，為蛇足石杉的保護提供了政策依據。結合研究結果來看，合理控制溫室氣體排放也有利于蛇足石杉的生存及其生物資源的保護。此外，蛇足石杉高適生區主要在貴州、重慶、湖南及浙江等地區，其高適生區質心在湖南靠近貴州方向，且存在向貴州方向轉移的趨勢，建議重點加強湖南及貴州的蛇足石杉野生資源的保護工作。

4 結 論

通過構建組合模型預測珍稀植物蛇足石杉當前和未來(2070)氣候情景下在中國的適生區并計算適生區及其質心變化，結果顯示，蛇足石杉當前總適生面積為226.99×104km2，中、高適生區面積154.51×104km2，質心均位于湖南省西部。蛇足石杉未來總適生區面積變化不大，中、高適生區劇烈收縮，收縮面積與溫室氣體和氣溶膠排放水平成正比，預測未來蛇足石杉中、高適生區將向非適生區和低適生區轉變。溫度和降水都對蛇足石杉適生區有重要影響，最濕月降水量(Bio13)、溫度季節性(Bio4)、最冷季均溫(Bio11)和最暖月最高溫(Bio5)是對其適生區影響最大的氣候變量。