末次間冰期以來及未來氣候情景下紅皮云杉適生分布區(qū)模擬*

王 亞 王軍輝 王福德 劉逸夫 譚燦燦 袁艷超 聶 穩(wěn)劉建鋒 常二梅 賈子瑞

（1.林木遺傳育種國家重點(diǎn)實(shí)驗(yàn)室 中國林業(yè)科學(xué)研究院林業(yè)研究所 北京 100091；2.中國林業(yè)科學(xué)研究院森林生態(tài)環(huán)境與自然保護(hù)研究所 北京 100091；3.黑龍江省林業(yè)科學(xué)研究所 哈爾濱 150081）

植物與氣候的相互關(guān)系一直是植物學(xué)和生態(tài)學(xué)研究的熱點(diǎn)（安芷生等，2005），氣候是決定物種分布的主要因素之一，物種分布格局變化能直觀反映氣候在不同時期的變化?；谖锓N分布模型（species distribution models, SDMs）和物種分布點(diǎn)的經(jīng)緯度數(shù)據(jù)與環(huán)境因子，通過數(shù)學(xué)方法可得到各環(huán)境因子的貢獻(xiàn)率，構(gòu)建物種潛在分布范圍，并以此推測未來和過去時期的適生分布區(qū)。物種分布預(yù)測在生態(tài)學(xué)、生物地理學(xué)等領(lǐng)域十分重要，已在土地利用變化、未來氣候變化情景（Faleiroet al.，2013）中得到廣泛應(yīng)用，用來研究生物多樣性損失（Bertrandet al.，2012）、瀕危物種管理和保護(hù)（Guoet al.，2016）、生物熱點(diǎn)區(qū)域確定（Walterset al., 2019）等。物種分布預(yù)測的可靠性受環(huán)境變量選取和模型參數(shù)的影響（Fourcadeet al.，2018;Sillero, 2011），環(huán)境變量之間的相關(guān)性會引入冗余信息并導(dǎo)致模擬結(jié)果偏差（Verbruggenet al.，2013）。因此，去除多余變量可消除變量間重線性關(guān)系，利于精確模擬物種的可能適生分布，并為生態(tài)保護(hù)等理論問題提供重要參考數(shù)據(jù)（Ashcroftet al.，2011）。

全球氣候變化背景下的物種適生分布區(qū)不斷變化。第四紀(jì)末次間冰期以來，全球氣候頻繁振蕩，極大影響物種種群動態(tài)變化（張愛平等，2018）。物種經(jīng)歷盛冰期和間冰期的不斷回旋，逐步形成當(dāng)前的地理分布格局（鄒旭等，2018）。物種分布模型研究表明，一些物種的分化和形成與歷史事件相關(guān)。對瀕危物種奇瓦瓦云杉（Picea chihuahuana）的最大熵模型分析發(fā)現(xiàn)，自末次盛冰期以來其適生區(qū)一直減少（Pinedo-Alvarezet al.，2019）；由于森林碎片化和氣候變化，其孤立種群在未來面臨消失的危險（Aguilar-Sotoet al.，2015）。對紫果云杉（Picea purpurea）、麗江云杉（Picea likiangensis）、青杄（Picea wilsonii）的潛在分布區(qū)預(yù)測發(fā)現(xiàn)，等溫性、最暖季降水量和年均溫是這3 個物種生態(tài)位分化的主要因子（張愛平等，2018），而且這3個物種在末次盛冰期發(fā)生擴(kuò)張。對云杉屬34 個物種的生態(tài)位分析表明，氣候?qū)υ粕夹螒B(tài)和物種分化的影響是由地理隔離、生態(tài)位保守性和生態(tài)適應(yīng)等許多生物和非生物因素介導(dǎo)的，溫度和降水分別是云杉系統(tǒng)發(fā)育和形態(tài)初級分化的主要驅(qū)動因素（Wanget al.，2017）。對東亞孑遺植物中全新世、末次盛冰期和未來的分布區(qū)預(yù)測結(jié)果發(fā)現(xiàn)，孑遺植物物種豐度和氣候避難所相對應(yīng)，中國西南部和越南北部是許多第三紀(jì)孑遺植物的長期避難所（Tanget al.，2018）。

紅皮云杉（Picea koraiensis）又名臭松，為松科（Pinaceae）云杉屬（Picea）植物，是中國東北重要的針葉樹種和古老的孑遺樹種，主要分布于大興安嶺北部、小興安嶺、長白山、完達(dá)山、張廣才嶺、老爺嶺等地區(qū)，地理坐標(biāo)范圍在41°40'—53°15'N、122°—131°E 之間，垂直分布范圍為250～1 100 m，屬當(dāng)?shù)刂饕炝峙c用材樹種。紅皮云杉分布于針闊混交林中，并處于針葉林到闊葉林的過渡帶，是研究植物分布響應(yīng)氣候變化的理想植物之一。紅皮云杉具有涵養(yǎng)水源、防風(fēng)固沙、水土保持等生態(tài)效益以及較高的經(jīng)濟(jì)價值，在木材生產(chǎn)和生態(tài)服務(wù)方面均有重要作用。目前，關(guān)于紅皮云杉的研究主要涉及生態(tài)學(xué)系統(tǒng)特征和生物量（張含國等，2003）、遺傳多樣性（王秋玉等， 2004）、礦質(zhì)元素含量（姜興林等，2002）等，不同時期氣候變化下的紅皮云杉地理分布格局變化尚未見報道。鑒于此，本研究基于最大熵模型（Merowet al.，2013），結(jié)合GIS 技術(shù)，利用末次間冰期（last interglacial, LIG, 120～140 ka BP）、末次盛冰期（last glacial maximum, LGM,21 ka BP）、全新世中期（mid-holocene, MH, ca.6 000 years ago）、現(xiàn)代（current, 1960—1990 年）和未來2 個時期（2050s 和2070s）的氣候數(shù)據(jù)和物種分布點(diǎn)數(shù)據(jù)，推測紅皮云杉過去、現(xiàn)代和未來的適生分布區(qū)，識別主要?dú)夂蛴绊懸蜃樱沂練夂驅(qū)t皮云杉地理分布格局形成的影響，推測種群遷移的歷史動態(tài)和生物避難所，并預(yù)測未來分布變化趨勢，以期為該物種的資源管理和多樣性保護(hù)提供科學(xué)依據(jù)。

1 材料與方法

1.1 物種分布數(shù)據(jù)收集

通過3 種方法收集紅皮云杉地理分布數(shù)據(jù)：第一，實(shí)地調(diào)查收集紅皮云杉的分布地點(diǎn)記錄，每個調(diào)查點(diǎn)之間的距離不小于20 km，記錄分布點(diǎn)的位置名稱、經(jīng)緯度和海拔；第二，檢索查詢相關(guān)文獻(xiàn)并記錄紅皮云杉分布的地理信息；第三，搜索中國虛擬植物標(biāo)本館（http://www.cvh.ac.cn/）和全球生物多樣性信息網(wǎng)絡(luò)（https://www.gbif.org/）等數(shù)據(jù)庫，收集紅皮云杉的分布點(diǎn)數(shù)據(jù)。匯總分布點(diǎn)信息，刪除地理信息不完整和重復(fù)的分布記錄，最終確定59 個紅皮云杉分布點(diǎn)。根據(jù)MaxEnt 模型的輸入文件格式，將物種分布點(diǎn)數(shù)據(jù)按照物種名、分布點(diǎn)經(jīng)度和緯度的順序進(jìn)行整理。

基于ENMTools v1.4 軟件（Warrenet al.，2010），使用Trim duplicate occurrences 函數(shù)對整理好的物種分布數(shù)據(jù)進(jìn)行冗余篩選，以避免地理聚集導(dǎo)致的空間自相關(guān)。最后共獲得55 個紅皮云杉有效分布點(diǎn)數(shù)據(jù)，用來構(gòu)建模型。從自然資源部網(wǎng)站（http://www.mnr.gov.cn/）下載中國地圖，使用ArcGIS v10.4.1 軟件將紅皮云杉分布點(diǎn)數(shù)據(jù)在地圖上展示出來（圖1）。

圖1 當(dāng)前紅皮云杉在中國的分布記錄Fig.1 Current distribution of P.koraiensis based on occurrence records

1.2 環(huán)境變量

云杉屬植物生長發(fā)育和繁殖受多種氣候因子影響，如溫度、降水等。本研究選取與溫度和降水相關(guān)的19 個環(huán)境變量進(jìn)行最大熵模型分析,其中，現(xiàn)代19個與溫度和降水有關(guān)的生物氣候變量（Bio1—Bio19）來自WoldClim version 1.4（Ficket al.，2017），精度為2.5 arc-minutes（弧分）。為避免因高度相關(guān)的氣候變量之間的相互影響而導(dǎo)致模擬結(jié)果過度擬合，利用MaxEnt v3.3.4 的刀切法（jackknife）將19 個氣候變量按貢獻(xiàn)率大小排序，篩選出貢獻(xiàn)率大于1 的變量。利用ENMTools 1.4 在19 個氣候變量之間進(jìn)行相關(guān)系數(shù)分析，一對相關(guān)系數(shù)大于|0.9|的氣候變量被視為非獨(dú)立變量。對每對顯著相關(guān)的變量，只選1 個有較大貢獻(xiàn)的變量用于分析。

1.3 過去和未來氣候數(shù)據(jù)

過去和未來時期的氣候變量數(shù)據(jù)從Worldclim version 1.4（Ficket al.，2017）下載，其中LGM 和MH 時期的數(shù)據(jù)精度為2.5 arc-minutes（弧分），LIG 時期的數(shù)據(jù)精度為30 arc-seconds（弧秒）。大氣環(huán)流模型選擇CCSM4（community climate system model four），CCSM4 是能夠模擬過去和未來氣候變化對植物分布影響的最有效的氣候模型之一（Abdelaalet al.，2019），已成功用于類似研究（Araret al.，2020; Liet al.，2020）。利用ArcGIS v10.4.1 將所有柵格文件轉(zhuǎn)換成ASCII 格式，將LIG 數(shù)據(jù)重抽樣為2.5 arc-minutes。下載未來2050s 和2070s 兩個時間段的氣候數(shù)據(jù)，選擇RCP4.5（representative concentration pathway）和RCP8.5 兩種排放氣候場景，用于模擬紅皮云杉未來地理分布格局。

1.4 最優(yōu)模型選擇

由于要素類型（feature class，F(xiàn)C）和正則化乘數(shù)（regularization multiplier，RM）2 個變量對模型結(jié)果影響較大，為獲得更精準(zhǔn)的適生區(qū)預(yù)測結(jié)果，本研究利用MaxEnt v3.3.4 軟件、R 3.6.2 程序和kuenm 程序包（Coboset al.，2019）對這2 個變量進(jìn)行優(yōu)化。首先利用R 軟件的kuenm 程序包將每個物種csv 格式的數(shù)據(jù)點(diǎn)分布文件拆分，70%的地理分布點(diǎn)用于訓(xùn)練集構(gòu)建、30%用于測試集驗(yàn)證。基于拆分后的地理分布點(diǎn)數(shù)據(jù)和篩選后的環(huán)境變量，利用R 語言的kuenm 軟件包從MaxEnt 模型的40 個正則化常數(shù)（0.1～4.0）和10種要素類型（L、Q、P、H 的2 種組合和3 種組合）中選擇最優(yōu)模型組合，要素類型含義如下：L（linear features）為線性，Q（quadratic features）為二次型，P（product features）為乘積型，H（hinge features）為片段化。利用kuenm 程序包逐一檢驗(yàn)上述400 組參數(shù)，根據(jù)最佳模型選擇標(biāo)準(zhǔn)（Coboset al.，2019），選擇統(tǒng)計學(xué)上顯著、遺漏率（omission rate）低于閾值（0.05）且ΔAICc（delta Akaike information criterion）值不高于2 的模型，具有ΔAICc 最小值的模型被作為最優(yōu)模型。

1.5 紅皮云杉潛在分布預(yù)測

應(yīng)用MaxEnt 最大熵模型預(yù)測紅皮云杉在過去、當(dāng)前和未來氣候情景下的潛在分布。將選定的物種分布點(diǎn)數(shù)據(jù)和環(huán)境變量導(dǎo)入MaxEnt v3.3.4 軟件。采用10 倍交叉驗(yàn)證法對模型進(jìn)行訓(xùn)練和驗(yàn)證，輸出格式為Logistic，輸出文件類型為asc，最大迭代次數(shù)設(shè)置為5 000，線程設(shè)置為10，重復(fù)次數(shù)為10 次。MaxEnt建模完成后，采用受試者工作特征曲線（receiver operating characteristic，ROC）評估模型，并使用曲線下面積（area under the ROC curve，AUC）測量模型可信度（Alloucheet al.，2006），AUC＞0.8 時表明物種適生區(qū)預(yù)測精度較好，AUC＞0.9 時表明預(yù)測精度非常好。本研究中，使用ArcGIS v10.4.1 軟件創(chuàng)建當(dāng)前、過去（LIG、LGM 和MH）和未來（2050s 和2070s）的適生區(qū)預(yù)測圖。MaxEnt 生成結(jié)果導(dǎo)入ArcGIS v10.4.1 軟件，并以中國地圖為底層進(jìn)行分析。使用自然斷點(diǎn)法劃分紅皮云杉的適生區(qū)等級：非適生區(qū)（0～0.1）、低適生區(qū)（0.1～0.3）、中適生區(qū)（0.3～0.5）和高適生區(qū)（0.5～1），計算各時期的適生面積。

1.6 適生分布區(qū)質(zhì)心轉(zhuǎn)移

為明確紅皮云杉適生分布區(qū)的變化趨勢，在ArcGIS v10.4.1 軟件中加載SDMToolbox （Brown，2014）工具箱，利用Centroid changes command 命令，計算不同時期的分布質(zhì)心。將這些質(zhì)心在ArcGIS v10.4.1 軟件中利用矢量箭頭投影，表示適生區(qū)變化的遷移距離和方向。

2 結(jié)果與分析

2.1 模型選擇和精度評價

為精準(zhǔn)預(yù)測適生區(qū)，本研究篩選出6 個環(huán)境變量，并選擇最優(yōu)的要素類型和正則化乘數(shù)，用于最大熵模型分析。6 個環(huán)境變量為年均氣溫（Bio1）、晝夜溫差與年溫差比值（Bio3）、氣溫季節(jié)性變動系數(shù)（Bio4）、降水量季節(jié)性變化（Bio15）、最干季度降水量（Bio17）、最暖季度降水量（Bio18），最優(yōu)模型要素類型為PH（product features, hinge features），最優(yōu)正則化乘數(shù)為2.3。10 次重復(fù)訓(xùn)練數(shù)據(jù)集和測試數(shù)據(jù)集的平均AUC 分別為為0.995 8 和0.995 7，均接近于1，說明模型有較高穩(wěn)定性和精度，模擬結(jié)果可信度高（圖2）。

圖2 基于MaxEnt 模型預(yù)測紅皮云杉分布的受試者工作特征曲線Fig.2 The ROC curve predicting the distribution of P.koraiensis based on MaxEnt

2.2 環(huán)境變量篩選

根據(jù)MaxEnt 模型的輸出環(huán)境因子貢獻(xiàn)率和刀切法分析結(jié)果，比較各環(huán)境變量對紅皮云杉潛在分布區(qū)變化的貢獻(xiàn)率（表1、圖3），最暖季度降水量（Bio18）、氣溫季節(jié)性變動系數(shù)（Bio4）和降水量季節(jié)性變化（Bio15）為3 個主要變量，貢獻(xiàn)率分別達(dá)40.6%、28.9%和24.6%，年均氣溫（Bio1）、最干季度降水量（Bio17）和等溫性（晝夜溫差與年溫差比值）（Bio3）的貢獻(xiàn)率分別為4%、1.5%和0.4%。

表1 影響紅皮云杉適生分布的環(huán)境變量貢獻(xiàn)百分率Tab.1 Contribution percentage of environmental variables affecting the suitable distribution of P.koraiensis

圖3 環(huán)境變量刀切法檢驗(yàn)Fig.3 Jackknife tests for environmental variables

重要變量響應(yīng)曲線（圖4）顯示出對紅皮云杉分布影響較大的變量（貢獻(xiàn)率＞2）適生范圍：最暖季度降水量（Bio18）為332.9～527.9 mm，氣溫季節(jié)性變動系數(shù)（Bio4）為129.4～173.5，降水量季節(jié)性變化（Bio15）為87.6～111.2，年均氣溫（Bio1）為-1.4～5.6 ℃。

圖4 紅皮云杉對重要環(huán)境變量的響應(yīng)曲線Fig.4 Response curves of P.koraiensis to important environmental variables

2.3 目前、過去和未來適生分布區(qū)

紅皮云杉當(dāng)前適生分布區(qū)包括中國的大興安嶺、小興安嶺、長白山脈及其支脈等地區(qū)，涉及黑龍江、吉林、遼寧等省，其中高適生區(qū)的分布范圍與實(shí)際分布相符（圖5）。紅皮云杉當(dāng)前的總適生區(qū)面積達(dá)133.58×104km2，其中，低、中、高適生分布區(qū)面積分別為40.93×104、42.36×104和50.29×104km2（表2）。

表2 紅皮云杉在多種氣候條件下的適生面積①Tab.2 Prediction of suitable area of P.koraiensis under different climate scenarios

圖5 當(dāng)代、過去和未來氣候情景下紅皮云杉在中國適生分布區(qū)Fig.5 Suitable distribution areas of P.koraiensis in China under contemporary, past and future climate conditions

歷史氣候條件下，紅皮云杉在末次間冰期（LIG）、末次盛冰期（LGM）和中全新世（MH）的適生分布區(qū)面積分別為目前的102.59%，105.20% 和100.53%；高適生區(qū)和總適生區(qū)的變化趨勢一致，末次盛冰期種群擴(kuò)張，全新世中期及以后種群持續(xù)收縮（圖5、表2）。

未來氣候情景下，紅皮云杉的總適生區(qū)和高適生區(qū)面積均呈減少趨勢，高適生區(qū)邊緣呈片段化分布，而小興安嶺腹地、長白山腹地、老爺嶺以及完達(dá)山等地的種群較為穩(wěn)定（圖5、表2）。紅皮云杉在未來氣候情景中2 種溫室氣體排放模式下（RCP4.5-2050s、RCP4.5-2070s、RCP8.5-2050s 和RCP8.5-2070s）的總適生區(qū)面積分別為當(dāng)前的99.06%、96.88%、98.05%和95.08%。在高排放（RCP8.5）模式下，紅皮云杉適生區(qū)面積縮減更大。RCP4.5 模型預(yù)測2050s 時期高適生區(qū)縮減的區(qū)域在大興安嶺北部邊緣（塔河縣）到中段（嫩江市）、長白山西部（輝南縣和撫松縣）；2070s 時期這2 個區(qū)域持續(xù)縮減，大興安嶺北部（蒙克山地區(qū)）和中部（高峰地區(qū)）、張廣才嶺北部（勃利縣和樺南縣）和東部（海林市）的高適生區(qū)消失。RCP8.5 模型預(yù)測2050s 時期高適生區(qū)縮減的區(qū)域在張廣才嶺北部（勃利縣和依蘭縣）和東部（海林市）；2070s 時期，紅皮云杉高適生區(qū)邊緣持續(xù)碎片化，長白山西部（輝南縣和撫松縣）的高適生區(qū)消失。

2.4 適生分布區(qū)質(zhì)心轉(zhuǎn)移

通過重建現(xiàn)代、過去和未來不同時期氣候條件下紅皮云杉的適生區(qū)質(zhì)心坐標(biāo)，分析質(zhì)心變化。紅皮云杉LIG 時期的分布質(zhì)心坐標(biāo)在46°14'27.38''N,127°1'46.92''E；LGM 時期分布質(zhì)心（46°10'54.01''N,127°14'49.34''E）向東南轉(zhuǎn)移21.62 km；MH 時期分布質(zhì)心（46°15'0.44''N, 127°18'54.83''E）向東北轉(zhuǎn)移8.96 km；當(dāng)前時期分布質(zhì)心（46°6'14.27''N, 126°14'28.81''E）向西南轉(zhuǎn)移104.43 km。未來氣候條件下，當(dāng)氣候情景為中排放濃度（RCP4.5）時，2050s年代分布質(zhì)心（46°25'2.47''N,127°17'54.03''E）向東北轉(zhuǎn)移105.72 km；2070s 年代分布質(zhì)心（46°15'33.81''N, 127°3'15.60''E）向西南轉(zhuǎn)移27.4 km；當(dāng)氣候情景為高排放濃度（RCP8.5）時，2050s年代分布質(zhì)心（45°49'27.43''N, 126°57'0.64''E）向東北轉(zhuǎn)移105.72 km；2070s 年代分布質(zhì)心（46°11'31.04''N,127°1'25.55''E）向西南轉(zhuǎn)移27.4 km（表3、圖6）。

表3 不同時期紅皮云杉適生區(qū)質(zhì)心坐標(biāo)和遷徙距離①Tab.3 Centroid coordinates and migration distance of suitable growth areas for P.koraiensis in different periods

圖6 不同時期紅皮云杉適生區(qū)質(zhì)心的遷移路線Fig.6 Migration route of centroid in the suitable growth area of P.koraiensis in different periods

3 討論

3.1 影響紅皮云杉分布的環(huán)境因子

根據(jù)自然選擇學(xué)說，物種形成是由于環(huán)境不同而引起的適應(yīng)性改變。環(huán)境是物種形成中最重要的因素，而全球氣候變化在環(huán)境變化中占有重要地位。本研究采用19 個環(huán)境變量，計算Pearson 相關(guān)系數(shù)，排除高相關(guān)性變量中貢獻(xiàn)率較低的變量，并通過AUC評價建模精度。本研究訓(xùn)練集和測試集的AUC 均高于0.99，說明模擬結(jié)果具有較高可靠性，當(dāng)前適生分布區(qū)的預(yù)測更接近實(shí)際情況。

氣候是決定物種分布變化的關(guān)鍵因素，其中溫度因子通過影響植物光合作用、呼吸作用來控制物種分布；而植物生長發(fā)育也受降水條件制約（Abeliet al.，2020）。本研究采用MaxEnt 自帶刀切法，不僅可提高模型預(yù)測可信度（De Cauweret al.，2014），而且能客觀、精確地評價19 個氣候變量。最暖季度降水量（Bio18）、氣溫季節(jié)性變動系數(shù)（Bio4）、降水量季節(jié)性變化（Bio15）對紅皮云杉的分布影響最大，其次為年均氣溫（Bio1）、最干季度降水量（Bio17）、晝夜溫差與年溫差比值（Bio3）等環(huán)境因子。其中，與水分相關(guān)的因子（Bio18、Bio15、Bio17）貢獻(xiàn)率為66.7%，說明與溫度因子相比，紅皮云杉適生區(qū)分布受降水影響更大，與王改妮（2018）、段國珍（2019）對紅皮云杉及其近緣種分布格局的模擬結(jié)果相似。這是因?yàn)榧t皮云杉在云杉屬物種分布區(qū)的最北端，北方河流冬季有封凍期，降水量較少，對于耐陰、喜濕潤的云杉屬植物來說，生長季降水量不足會限制生長，長期缺水則會導(dǎo)致植物的適應(yīng)性分化。與其近緣種白杄（Picea meyeri）和沙地云杉（Picea mongolica）相比，紅皮云杉適應(yīng)濕冷環(huán)境的習(xí)性十分明顯，無論是大氣濕度還是土壤水分供給狀況均高于前二者。與云杉組（Sect.Picea）內(nèi)其他類群相比，紅皮云杉是對水分需求較高的種群，其分布區(qū)內(nèi)的年均降水量僅次于云杉（Picea asperata）和青扦（Picea wilsonii）；與同域種魚鱗云杉（Picea jezoensis）（云杉屬魚鱗云杉組sect.Omorica）相比，紅皮云杉更適應(yīng)低濕的土壤環(huán)境，在排水良好的山坡，紅皮云杉的競爭力不及魚鱗云杉，分布數(shù)量較少。此外，影響擦木（Sassafras tzumu）、柳杉（Cryptomeria japonica）等植物適生區(qū)分布的主要限制因子也是降水（劉想等，2018；戴旻峻等，2022），這些植物穩(wěn)定生長在山地溝谷和水分肥沃的地區(qū)，同樣，紅皮云杉主要生長在河谷溪旁、河岸階地和山麓。綜上表明，紅皮云杉對水分的需求較高，降水是影響紅皮云杉分布和擴(kuò)散潛力的主要限制因子，因此紅皮云杉在推廣種植或異地保護(hù)時，可選擇在降水充沛地區(qū)。

溫度在影響紅皮云杉分布上也起到較重要作用。根據(jù)刀切法結(jié)果，年均氣溫（Bio1）對紅皮云杉分布的影響不可替代，Bio1 的適宜范圍在-1.4～5.6 ℃，最適值為4.3 ℃，在-1.4～4.3 ℃時，分布概率隨年均溫升高而增大，在4.3～5.6 ℃時，分布概率隨年均溫升高而降低，該因子可能限制紅皮云杉繼續(xù)向南擴(kuò)張。氣溫季節(jié)性變動系數(shù)（Bio4）的貢獻(xiàn)率達(dá)28.9%，Bio4 的適宜范圍在129.4～173.5，最適值為155.6。Bio4 代表溫度在一年中的變化幅度，數(shù)值越大溫度變幅越大，因此均一穩(wěn)定的生長環(huán)境可能對促進(jìn)紅皮云杉種群擴(kuò)張有積極作用。

3.2 紅皮云杉在不同氣候情景下的分布格局變化和質(zhì)心轉(zhuǎn)移

植物在與氣候的長期交互作用下形成一定的地理分布格局。通常，擁有更廣生態(tài)位的物種具有更強(qiáng)的環(huán)境適應(yīng)能力，可充分利用其周邊的自然資源（Yinget al.，2016），因此在預(yù)測紅皮云杉林分更新狀況時必須考慮其對氣候變化的適應(yīng)情況。本研究采用數(shù)據(jù)庫數(shù)據(jù)與實(shí)地調(diào)查數(shù)據(jù)相結(jié)合的方法，利用古氣候數(shù)據(jù)、現(xiàn)代氣候數(shù)據(jù)和未來氣候變化情景，構(gòu)建紅皮云杉物種分布模型，并重建LIG、LGM 和MH 3個歷史時期的地理分布格局。

研究發(fā)現(xiàn)，紅皮云杉種群適宜生境曾經(jīng)歷多次適應(yīng)性轉(zhuǎn)換，包括末次盛冰期的種群擴(kuò)張以及全新世中期及之后的種群持續(xù)收縮。末次間冰期的陸地平均氣溫比現(xiàn)代高1～2 ℃（何韜，2019），與現(xiàn)代更相似；末次盛冰期的陸地平均氣溫比現(xiàn)代低5～10 ℃，而中高緯度地區(qū)低20 ℃以上（何韜，2019）。本研究結(jié)果和Stewart 等（2010）提出的理論是吻合的，即適應(yīng)寒冷的植物在間冰期的分布范圍比盛冰期小，如白云杉（Picea glauca）（Stewartet al.，2010 ；Bannisteret al.，2001），以往對云杉屬其他物種的生態(tài)位模擬分析中也有一致的結(jié)果（Shenet al.，2019；王改妮，2018；段國珍，2019）。相對于末次間冰期，末次盛冰期氣候較冷，紅皮云杉的質(zhì)心向東南方向轉(zhuǎn)移；全新世中期氣候再次回暖，紅皮云杉的質(zhì)心向東北方向轉(zhuǎn)移。紅皮云杉的質(zhì)心轉(zhuǎn)移與氣候變暖條件下物種向溫度更低的區(qū)域遷移的規(guī)律是相符的。

溫度適當(dāng)升高可以加速植物生長發(fā)育，并影響種群數(shù)量變化。隨著全球溫度持續(xù)上升，在未來氣候變化2 種溫室氣體排放情景下，紅皮云杉的適生區(qū)面積均呈現(xiàn)持續(xù)減少趨勢，高度適生區(qū)呈片段化分布，全球氣候變暖不利于紅皮云杉適生區(qū)擴(kuò)張。對附子（Aconitum carmichaeli）（韋潔敏，2020） 、 金錢松（Pseudolarix amabilis）（張文秀，2020）、桫欏（Alsophila spinulosa）（楊啟杰等，2021）等植物的潛在分布預(yù)測分析也發(fā)現(xiàn)，隨著未來全球氣候變暖，這些物種適生區(qū)面積會持續(xù)減少，與本研究結(jié)論一致。未來氣候變化情景中，與RCP4.5 排放路徑相比，RCP8.5 排放路徑下的適生區(qū)面積縮減速度較大，紅皮云杉在高排放情景下適生區(qū)面積變化更顯著，對氣候變化響應(yīng)明顯。隨著未來氣候變暖，紅皮云杉總適生區(qū)質(zhì)心先向高緯度遷移，后逐漸向低緯度回落，出現(xiàn)這種情況的主要原因是紅皮云杉潛在適生區(qū)在向西南方向海拔更高的長白山地區(qū)擴(kuò)張。全球變暖影響植物的多樣性、分布和遷移模式，大多數(shù)物種遷移到高海拔和高緯度地區(qū)（Parmesan，2006）。木梨（Pyrus xerophila）、四合木（Tetraena mongolica）未來的潛在分布范圍有向西北方向的高緯度地區(qū)遷移的趨勢（劉超等，2020；段義忠等，2019）。未來氣候情景下，紅皮云杉的高適生區(qū)縮減，邊緣呈片段化分布，縮減區(qū)域包括大興安嶺的塔河、嫩江、蒙克山和高峰，長白山西部的輝南和撫松，張廣才嶺的勃利、樺南、依蘭和海林等地，建議對這些區(qū)域的種群加強(qiáng)就地保護(hù)、建立自然保護(hù)區(qū)、減少人為破壞與干擾。小興安嶺腹地和長白山腹地的種群生境較為穩(wěn)定，是建立紅皮云杉異地保護(hù)的最佳候選區(qū)域。

3.3 紅皮云杉避難所推測

當(dāng)全球經(jīng)歷第四紀(jì)冰期的漫長考驗(yàn)時，眾多生物依靠生物避難所的保護(hù)而得以幸存，并在冰期后重新擴(kuò)散。識別生物避難所有3 個主要途徑：一是基于生物地理模式，通過收集目標(biāo)生物體的古生物學(xué)、生態(tài)學(xué)和遺傳學(xué)數(shù)據(jù)，量化生物避難所的時空緯度；二是調(diào)查物種避難所棲息地的產(chǎn)生過程；三是通過模擬物種分布格局對環(huán)境變化的響應(yīng)來識別避難所（胡忠俊等，2013）。冰期避難所理論表明，避難所在環(huán)境波動期間的生態(tài)條件相對穩(wěn)定，可促進(jìn)物種多樣性的積累。因此物種具有較高水平遺傳多樣性的地理分布區(qū)域，可以作為尋找避難所的依據(jù)之一（石甜，2021）。

紅皮云杉形成于第三紀(jì)中晚期（4.12～7.64 Ma, 百萬年）(Shaoet al.，2019)，說明最先引起紅皮云杉群體遺傳分化的因素是第三紀(jì)的地質(zhì)和氣候事件。在第三紀(jì)中晚期，青藏高原出現(xiàn)快速抬升（8～10 Ma），加速了亞洲大陸內(nèi)部的干燥度和亞洲季風(fēng)的產(chǎn)生（Anet al.，2001）。第三紀(jì)末發(fā)生地殼運(yùn)動，使得大興安嶺沿東側(cè)斷層翹起；第四紀(jì)初期的構(gòu)造運(yùn)動又引起火山噴發(fā)，在小興安嶺地區(qū)形成大片玄武巖熔巖臺地，同時有輕微翹起，沿北西方向斷裂隆起（樊祺誠等，2015；劉建強(qiáng)等，2017）；第四紀(jì)中晚期至全新世，長白山發(fā)生多次火山噴發(fā)，長白山主體于全新世形成后，地殼相對穩(wěn)定（王團(tuán)華，2006）。紅皮云杉的分化時間段發(fā)生于第三紀(jì)末至第四紀(jì)期間，東亞季風(fēng)氣候增強(qiáng)和地殼變動的時間基本一致。與其分化時間接近的還有第三紀(jì)植物楓香（Liquidambar formosana）（10.30 Ma）（孫榮喜，2017）、水青樹（Tetracentron sinense）（9.6 Ma）（Sunet al.，2014）、青岡棟屬（Cyclobalanopsis）（9.07 Ma）（Denget al.，2018）等。

植物在冰期會退縮到避難所，耐寒植物在冰期往往會擴(kuò)張且在間冰期發(fā)生種群收縮，故生物避難所還應(yīng)包括間冰期避難所（Stewartet al.，2010）。紅皮云杉屬于耐寒植物，通過分子證據(jù)、適生分布區(qū)模擬結(jié)果以及孢粉化石證據(jù)可推測紅皮云杉的間冰期避難所。分子證據(jù)顯示，在紅皮云杉不同地理種群中，長白山腹地的臨江種群（π= 0.002 35）和天橋嶺種群（π =0.002 32）均具有較高的核苷酸多樣性（nucleotide diversity, π）(Wanget al.，2023)；此外，由于多個歷史時期的物種分布重疊區(qū)往往是植物潛在避難所（Carnavalet al.，2008），結(jié)合對末次間冰期和末次盛冰期的適生分布區(qū)模擬結(jié)果推測，長白山腹地可能是紅皮云杉第四紀(jì)間冰期避難所，且在該種群進(jìn)化過程中發(fā)揮了重要作用。從生境特征上看，長白山地區(qū)屬于紅皮云杉分布區(qū)內(nèi)海拔最高的區(qū)域（海拔100 0～274 9 m），較大的海拔落差為紅皮云杉等山地植物適應(yīng)氣候變化提供了充足的緩沖環(huán)境，且該區(qū)域末次間冰期時的氣候濕冷（孫廣友等，2012），為紅皮云杉提供了生長所需的充沛水分。此外，小興安嶺腹地也可能是紅皮云杉的避難所之一。分子證據(jù)顯示小興安嶺腹地的種群（沾河）核苷酸多樣性較高（π= 0.002 31）（Wanget al.，2023）；此地又是多個時期適生分布區(qū)的重疊區(qū)域；從孢粉數(shù)據(jù)來看，松遼盆地發(fā)現(xiàn)大量中更新世晚期的云杉屬孢粉化石（石寧，1996），與末次間冰期的時間吻合；小興安嶺處于華力西褶皺帶，第三紀(jì)時山體北部與東北方的結(jié)雅河低地、西南方的松嫩平原聯(lián)成一片，在上新世末更新世初沿北東向和北西向斷裂抬升為山嶺，與結(jié)雅河低地及松嫩平原分離（王歡等，2018），其復(fù)雜多變的山地峽谷地貌為紅皮云杉提供了豐富的避難場所。因此，結(jié)合分子證據(jù)、適生區(qū)模擬、孢粉化石和地形地貌，推測紅皮云杉在間冰期的生物避難所在小興安嶺腹地和長白山腹地。

由于氣候數(shù)據(jù)是根據(jù)數(shù)學(xué)方法估計所得，與古代實(shí)際氣候存在差異，且未考慮非生物因素的影響，因而MaxEnt 模型在適生分布區(qū)預(yù)測方面存在一定局限性。然而，通過與其他物種分布模式的比較，MaxEnt具有較高的可靠性和穩(wěn)定性（Phillipset al.，2006；Elithet al.，2010；朱耿平等，2016）。在模擬物種當(dāng)代適生分布區(qū)時，其結(jié)果和物種實(shí)際分布范圍具有高度一致性（高文強(qiáng)等，2016；趙光華等，2021；王東升等，2022）。該方法可以呈現(xiàn)物種在不同時期的地理分布格局，判斷環(huán)境因素對物種分化和形成的影響，推測物種的生物避難所，預(yù)測物種在未來氣候下的分布范圍變化趨勢并為物種保護(hù)提供依據(jù)（李國慶等，2013；朱耿平等，2013）。

4 結(jié)論

MaxEnt 模型預(yù)測紅皮云杉適生分布區(qū)結(jié)果可靠，訓(xùn)練集和測試集的AUC 均大于0.99。與溫度因子相比，降水對紅皮云杉的地理分布影響更重要，在其推廣種植或者異地保護(hù)時，可選擇降水充沛地區(qū)。溫度因子可能限制紅皮云杉向南擴(kuò)張，均一穩(wěn)定的生長環(huán)境對促進(jìn)其種群擴(kuò)張有積極作用。由于第三紀(jì)末至第四紀(jì)期間的地質(zhì)和氣候事件，紅皮云杉的適生區(qū)經(jīng)歷多次轉(zhuǎn)換，包括末次盛冰期種群擴(kuò)張、全新世中期及之后的種群持續(xù)收縮。紅皮云杉的適生分布區(qū)質(zhì)心轉(zhuǎn)移與氣候變暖條件下物種向溫度更低的區(qū)域遷移的規(guī)律相符。全球氣候變暖不利于紅皮云杉的種群擴(kuò)張，未來情景下紅皮云杉的適生區(qū)呈縮減趨勢。應(yīng)考慮優(yōu)先保護(hù)大興安嶺北部到中部、長白山西部和張廣才嶺北部等高適生區(qū)的邊緣種群。根據(jù)分子證據(jù)、適生區(qū)模擬，結(jié)合孢粉化石和地形地貌，推測小興安嶺腹地和長白山腹地是紅皮云杉的間冰期避難所。