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

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

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

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

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

1 材料和方法

1.1 蛇足石杉分布數(shù)據(jù)

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

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

1.2 氣候數(shù)據(jù)的獲取與篩選

從世界氣候數(shù)據(jù)庫(WorldClim, https://www.worldclim.org/)下載20世紀(jì)中、后期(1960～1990)的氣候數(shù)據(jù)作為當(dāng)前氣候數(shù)據(jù)和對中國氣候模擬能力較好的CCSM4模式[33]下的21世紀(jì)下半葉(2070)的RCP2.6、RCP4.5、RCP6.0及RCP8.5這4種典型濃度路徑(representative concentration pathways, RCPs)的氣候數(shù)據(jù)作為未來氣候數(shù)據(jù)，氣候數(shù)據(jù)分辨率均為2.5 minutes。RCPs表示2100年相對于1750年的輻射強(qiáng)迫總和，向大氣排放的溫室氣體和氣溶膠會影響地球輻射平衡產(chǎn)生輻射強(qiáng)迫[34]。RCP2.6輻射強(qiáng)迫為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]。

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

表1 19個生物氣候變量

1.3 物種分布模型的構(gòu)建與評估

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

表2 10種用于模擬的算法

1.4 物種適生區(qū)及其變化預(yù)測

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

2 結(jié)果與分析

2.1 模型預(yù)測精度評估

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

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

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

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

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

2.3 蛇足石杉適生區(qū)預(yù)測

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

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

圖3 當(dāng)前蛇足石杉適生區(qū)

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

2.4 蛇足石杉適生區(qū)及質(zhì)心變化

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

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

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

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

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

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

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

A.蛇足石杉總適生區(qū)質(zhì)心變化；B.蛇足石杉中、高適生區(qū)質(zhì)心變化

3 討 論

3.1 蛇足石杉的未來分布

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

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

3.2 蛇足石杉的保護(hù)

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

4 結(jié) 論

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