王國崢,耿其芳,肖孟陽,張孟源,張云燕,王中生

南京大學(xué)生命科學(xué)學(xué)院,南京 210023

金錢松(Pseudolarixamabilis(Nelson) Rehd.)是我國特有單屬種松科植物,也是著名的第三紀(jì)孑遺植物,現(xiàn)主要分布于江蘇南部、安徽南部、福建北部、浙江、江西、湖南、湖北利川至四川萬縣交界地區(qū)[1]。最早的化石記錄出現(xiàn)于西伯利亞東北部的晚白堊紀(jì)地層中,始新世時出現(xiàn)于挪威斯匹次卑爾根島西部,美國、歐洲、亞洲中部以及日本在第三紀(jì)的不同時期都曾發(fā)現(xiàn)有該屬化石分布[2]。據(jù)化石資料推測,金錢松屬可能起源于白堊紀(jì)中世紀(jì)晚期,第三紀(jì)時曾發(fā)展為北半球中高緯度廣布屬,其中在日本和歐洲一直維持至更新世第一次冰期降臨[3]。在新生代末期冰期,它從歐洲發(fā)生地向我國南方遷移,最后定居于長江流域,形成現(xiàn)今的地理分布格局[4]。

由于金錢松特殊的進化及遷移歷史,其在松科植物系統(tǒng)發(fā)育研究中備受關(guān)注,同樣對古生態(tài)以及古氣候的研究具有重要意義[5]。從其化石記錄可發(fā)現(xiàn)金錢松的分布與地質(zhì)變化過程密切相關(guān),地質(zhì)變化導(dǎo)致氣候劇烈變化,金錢松分布范圍隨之發(fā)生改變。一些學(xué)者已從金錢松植物學(xué)特征[5]、種群生態(tài)學(xué)[6]、以及遺傳多樣性[7- 8]等方面開展了相關(guān)研究,其中高燕會[7]和劉俊[8]等利用分子標(biāo)記發(fā)現(xiàn)金錢松天然種群具有較高的遺傳多樣性,認(rèn)為遺傳因素不是金錢松的致瀕因子,其呈狹域分布的主要原因可能來自外界自然因素。

近年來生態(tài)位模型在預(yù)測、解析物種分布范圍及其變化等方面應(yīng)用廣泛,利用物種已知的分布數(shù)據(jù)和環(huán)境變量,根據(jù)一定的算法運算來構(gòu)建模型,并將運算結(jié)果投射到不同的時間和空間中以預(yù)測物種的潛在適生區(qū)[9]。目前較為常用的生態(tài)位模型較多,每種模型都可以獨立預(yù)測出潛在適生區(qū),但各模型都存在一定的偏好性[10],本文參考Miguel B. Araújo[10]的方法,采用了4種基于不同算法的模型(GARP、Bioclim、Domain和Maxent)預(yù)測金錢松的全國潛在適生區(qū),利用集合預(yù)測系統(tǒng)思想,綜合多種模型(3種以上)預(yù)測結(jié)果,盡可能降低由經(jīng)驗性選擇一種模型導(dǎo)致的假陰性或假陽性影響,同時一種模型的缺陷可能被另一種模型所彌補,從而提高預(yù)測的科學(xué)性。GARP是基于遺傳算法的規(guī)則組合進行局域環(huán)境空間建模[11]；Bioclim基于矩形框架模型,每一種環(huán)境因子被抽象成一個矩形框,絕大多數(shù)分布點(95%以上點)位于這些框內(nèi),多個矩形框構(gòu)成限制范圍,若某點位于這個范圍內(nèi),便可以認(rèn)為其適宜該物種生存[12]；Domain模型基于Gower算法,通過計算所有已知分布點之間的最大相似度,來評估預(yù)測點對于環(huán)境變量的適宜程度[13]；Maxent模型通過物種的分布數(shù)據(jù)和環(huán)境數(shù)據(jù),找出物種分布規(guī)律的最大熵,從而對物種的分布進行預(yù)測[14- 20]。

1 數(shù)據(jù)獲取和研究方法

1.1 金錢松樣本分布數(shù)據(jù)

金錢松的地理分布位點主要通過中國數(shù)字植物標(biāo)本館(http://www.cvh.org.cn/)以及中國科學(xué)院北京植物所標(biāo)本館(http://pe.ibcas.ac.cn/)查閱,共獲取44個分布記錄數(shù)據(jù),準(zhǔn)確到縣。用Google Earth 轉(zhuǎn)換出分布點的經(jīng)緯度坐標(biāo),輸入Excel保存成.csv格式。

1.2 環(huán)境變量數(shù)據(jù)及預(yù)處理

環(huán)境變量共22個,包括19個氣候因子和3個地形因子,19個氣候因子來源于Worldclim (http://www.worldclim.org/),包括末次盛冰期(Last Glacial Maximum,LGM)、當(dāng)前(Current)和2070年的數(shù)據(jù),分辨率2.5 arc-minutes,將環(huán)境因子數(shù)據(jù)用ArcGIS 10.2統(tǒng)一轉(zhuǎn)化為ASCII格式以便使用。

由于環(huán)境變量之間具有多重共線性,會導(dǎo)致預(yù)測分布過度擬合[21],因此本文參考張琴[15]對環(huán)境共線性診斷的方法,對環(huán)境因子進行Spearman秩相關(guān)分析,篩選出相關(guān)性較低的環(huán)境變量(Spearman系數(shù)<0.75),并在Spearman系數(shù)大于0.75的變量中選出具有生態(tài)學(xué)重要意義的環(huán)境變量,共得到10個變量(表1)。

表1 環(huán)境變量

1.3 相關(guān)軟件

4種生態(tài)位模型所用軟件：基于GARP模型預(yù)測的Desktopgarp(Version1.1.6),基于Bioclim和Domain模型預(yù)測DIVA-GIS(Version 7.5)；基于最大熵模型預(yù)測軟件Maxent(Version 3.3.3)。

數(shù)據(jù)處理軟件有：Excel用于分布數(shù)據(jù)記錄以及數(shù)據(jù)格式轉(zhuǎn)換,DIVA-GIS(Version 7.5)用于圖層格式轉(zhuǎn)換,ArcGIS(Version 10.2)用于圖層的處理和數(shù)據(jù)的格式轉(zhuǎn)換。

1.4 模型預(yù)測與檢驗方法

1.4.1GARP模型預(yù)測

將物種已知分布點數(shù)據(jù)轉(zhuǎn)換成Desktopgarp模型軟件支持的格式,通過模型軟件中 Upload Data Points功能加載,選擇70%的分布數(shù)據(jù)作為訓(xùn)練數(shù)據(jù),其余數(shù)據(jù)作為檢驗數(shù)據(jù)[14]。環(huán)境數(shù)據(jù).asc文件通過Desktopgarp中的Dataset Manager處理轉(zhuǎn)換成Desktopgarp可識別的格式.raw,以數(shù)據(jù)集的形式加載到軟件中。

1.4.2Maxent模型預(yù)測

將已知分布點數(shù)據(jù)文件格式在Excel中轉(zhuǎn)為.csv格式文件導(dǎo)出,通過Browse加載到Maxent軟件中。再將10個環(huán)境變量ASCⅡ文件的環(huán)境數(shù)據(jù)通過Browse加載到Maxent軟件中。并進一步利用Maxent軟件對金錢松在末次盛冰期(LGM)氣候和未來氣候(2070年)條件下進行模擬分布[22],再用ArcGIS 10.2 SDMs Tool模擬金錢松從末次盛冰期到當(dāng)前的分布遷移變化[23]。

1.4.3Bioclim和Domain模型分析

在DIVA-GIS中首先添加.shp格式的訓(xùn)練數(shù)據(jù)集,再將10個環(huán)境變量文件ASCⅡ轉(zhuǎn)換成.grd格式,生成stack數(shù)據(jù)集。在Modeling-Bioclim/Domain模塊中添加stack格式的環(huán)境數(shù)據(jù)集,依次進行Bioclim和Domain模型的預(yù)測[15]。

1.5 模型預(yù)測結(jié)果檢驗

受試者工作特征曲線(Receiver Operating Characteristic,ROC)和Kappa統(tǒng)計量廣泛用于模型精度的檢驗。本文參考張琴[15]的方法,利用DIVA-GIS軟件,將10組不同的訓(xùn)練和測試數(shù)據(jù)用4種模型預(yù)測得到的40個預(yù)測分布圖生成stack數(shù)據(jù)集,再導(dǎo)入軟件創(chuàng)建的驗證點.shp文件,最后得出各模型的AUC值和Kappa值。

受試者工作特征曲線下的面積即為AUC值(Area under recriver operating characteristic curve,AUC),AUC的數(shù)值范圍為0.5—1.0,“0.5”說明預(yù)測結(jié)果為隨機分布,預(yù)測值越大表示該模型預(yù)測結(jié)果越精確。一致性檢驗統(tǒng)計量(Kappa)是對兩種方法測定結(jié)果的一致部分進行檢驗,取值范圍是[-1,1],“-1”說明兩次判斷的結(jié)果完全一致,“1”說明兩次判斷的結(jié)果完全不一致,值越大一致性越好,通用的經(jīng)驗法則是Kappa大于0.75表示一致性好,小于0.4表示一致性差[24]。

1.6 影響適生分布的環(huán)境因子分析

利用Maxent軟件中的“刀切法”(Do jackknife to measure variable importance)可以顯示出環(huán)各個境變量對預(yù)測金錢松分布的貢獻。其他3種模型(GARP、Bioclim、Domain)預(yù)測軟件不具備分析環(huán)境因子對預(yù)測分布貢獻的功能,因此將GARP、Bioclim、Domain等3種模型基于10個環(huán)境因子預(yù)測的適生區(qū)作為對照組,分別以“刀切法”分析影響較大的三個環(huán)境因子為單因子,或作為共同因子來預(yù)測金錢松的適生區(qū),進行預(yù)測結(jié)果的比較。

2 結(jié)果與分析

2.1 4種模型預(yù)測結(jié)果

從GARP模型預(yù)測圖中(圖1)可以看出金錢松在我國的適生區(qū)集中在江蘇、浙江、安徽、江西、河南、湖北、湖南、貴州、重慶、四川東部、陜西南部以及云南北部地區(qū),集中分布在24.43°—33.35°N和106.41°—123.42°E之間,而搜集的樣本數(shù)據(jù)在27.27°—34.05°N和108.41°—121.80°E之間。

Maxent預(yù)測分布圖(圖1)顯示金錢松高適區(qū)主要集中在江蘇南部、安徽南部、湖北南部、江西北部以及浙江西北部。綠色代表中度適生區(qū),依次擴散開來,最大范圍向北達江蘇省潥陽市,向南至湖南衡山磨鏡臺,最西可達重慶市萬州區(qū),最東可達浙江東天目山。

通過DIVA-GIS的Modeling功能運算得出的Bioclim模型預(yù)測圖(圖1)和Domain模型預(yù)測圖(圖1),可觀察到代表最適區(qū)的紅色區(qū)域在這兩種模型中占較小比例,Bioclim預(yù)測分布范圍較小,顯示金錢松可在浙江、湖南、安徽南部、湖北南部以及江西西北部高度適生。Domain模型預(yù)測結(jié)果分布范圍較大,總體與GARP模型預(yù)測分布范圍接近,高適與最適分布范圍與Maxent相似,但在河南新鄉(xiāng)、南陽以及云南昆明出現(xiàn)最適分布區(qū)。

2.2 不同模型預(yù)測準(zhǔn)確性評估

本文采用10組訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)對4個模型進行ROC曲線分析和Kappa一致性檢驗。GARP模型AUC平均值為0.922,Maxent模型AUC平均值為0.980,Bioclim模型AUC平均值為0.957,Domain模型的AUC平均值為0.940(見表2),均遠(yuǎn)遠(yuǎn)超過隨機模型(AUC=0.5)。并且4種模型的Kappa均值都大于0.75,表明4種模型預(yù)測的一致性均較顯著,預(yù)測精度較高。4種模型AUC值的標(biāo)準(zhǔn)偏差排序為Bioclim > Domain > Maxent >GARP。

表2 4種模型的AUC值和Kappa值

2.3 影響適生區(qū)的環(huán)境因子分析

各個境變量對預(yù)測金錢松分布的貢獻如圖2,在10個環(huán)境因子中,年均溫(Bio1)對于金錢松分布的貢獻最大,年降水量(Bio12)、最濕月降水量(Bio13)和最熱月極高溫度(Bio5)對金錢松分布的影響也比較大,溫度日較差(Bio2)和海拔對金錢松分布的影響比較小。

圖1 4種生態(tài)位模型的金錢松適生區(qū)預(yù)測圖Fig.1 Potential distribution of P.amabilis based on four ecological niche models

圖2 刀切法檢測生物氣候變量對分布的重要性Fig.2 Import of environmental variables for prediction based on jackknife test

從單因子預(yù)測結(jié)果來看(圖3),3種模型基于年均溫(Bio1)的模擬結(jié)果與各自預(yù)測適生區(qū)最為接近,但范圍較適生區(qū)均偏大。橫向?qū)Ρ?種模型分別基于Bio1、Bio12、Bio13等3種單因子的模擬結(jié)果,基于年均溫(Bio1)的預(yù)測結(jié)果相對較為接近,說明年均溫(Bio1)對這3種模型的預(yù)測結(jié)果也都有顯著的影響?；谀杲邓?Bio12)與最濕月降水量(Bio13)兩種因子的預(yù)測分布圖都與適生區(qū)差異較大,且三種模型之間對比預(yù)測分布圖差異也較大,說明這兩種因子不是最主要的預(yù)測分布影響因子。將Bio1、Bio12、Bio13等作為共同影響因子,3種模型的預(yù)測分布圖與其對應(yīng)的適生區(qū)都十分接近,說明這3種因子是預(yù)測分布結(jié)果的確存在顯著影響,尤其是Bioclim模型,二者基本完全重疊。此外,基于年均溫(Bio1)的預(yù)測分布與3因子聯(lián)合預(yù)測分布最為接近,進一步證明年均溫(Bio1)對金錢松預(yù)測分布起“框架”作用,而其余環(huán)境因子對于預(yù)測分布起到“修飾”作用。

圖3 三種環(huán)境因子模擬結(jié)果Fig.3 The results of the other three models圖a1—a5分別為GARP模型基于年均溫(Bio1)、年降水量(Bio12)、最濕月降水量(Bio13)、三因子(Bio1&12&13)以及對照組結(jié)果,圖b1—b5分別為Domain模型基于年均溫(Bio1)、年降水量(Bio12)、最濕月降水量(Bio13)、三因子(Bio1&12&13)以及對照組結(jié)果,圖c1—c5分別為Bioclim模型基于年均溫(Bio1)、年降水量(Bio12)、最濕月降水量(Bio13)、三因子(Bio1&12&13)以及對照組結(jié)果

2.4 金錢松不同時期分布格局變化模擬

金錢松從末次盛冰期到當(dāng)前,分布范圍主要向南遷移,分布擴張面積達到178510.08 km2,北部分布范圍有部分收縮,達到122865.39 km2,凈分布范圍面積擴大了55644.69 km2。從當(dāng)前到2070年分布范圍呈現(xiàn)收縮趨勢,主要是南部適生區(qū)大面積收縮,達到246396.07 km2,占現(xiàn)生分布范圍的14.45%,而向北擴張68311.77 km2,僅為收縮面積的27.72%,分布范圍收縮了178084.30 km2(圖4)。

圖4 基于Maxent和ArcGIS金錢松不同時期的分布格局變化Fig.4 Comparison of the distribution changes of P.amabilis based on ArcGIS and Maxent

3 討論

3.1 金錢松適生范圍預(yù)測

近年來物種分布模型發(fā)展迅速,多個可用的方法模型相繼被開發(fā)[21]。本文通過模型檢驗得出四個模型都可以較為準(zhǔn)確地預(yù)測出金錢松的潛在適生區(qū)(AUC>0.9),但GARP預(yù)測結(jié)果的AUC均值最小,精度較差,這可能是金錢松樣本數(shù)據(jù)較少導(dǎo)致。楊會楓[19]用不同樣本容量檢驗GARP預(yù)測結(jié)果精度,發(fā)現(xiàn)受到樣本容量干擾,只有當(dāng)樣本容量超過一定閾值,AUC均值才會趨于穩(wěn)定(大于0.9)。Domain和Bioclim預(yù)測結(jié)果相比其他模型適生區(qū)范圍較小,這可能與模型受樣點信息有關(guān),金錢松分布數(shù)據(jù)主要來自于標(biāo)本館,樣本數(shù)據(jù)存在一定偏好性,科研工作者一般根據(jù)自己研究所需進行采樣或者單純?yōu)闃?biāo)本館收集標(biāo)本,樣本信息較為零散,缺乏系統(tǒng)性和代表性；其次金錢松樣本數(shù)據(jù)較少可能導(dǎo)致生態(tài)位空間縮小的現(xiàn)象[26],并且由于某些地區(qū)金錢松種群遺傳多樣性較為豐富[5- 6],其物種在不同區(qū)域的生態(tài)位也可能發(fā)生漂移[26]。

Maxent模型預(yù)測結(jié)果比其他模型更為精確(AUC平均值為0.980),相關(guān)文獻也證明Maxent在大、小樣本數(shù)據(jù)條件下均能很好地預(yù)測物種適生區(qū)[14- 20]。但不同模型預(yù)測結(jié)果可以相互補充[21,27- 29],如在Maxent預(yù)測分布圖上,河南鄭州、信陽大別山以及云南昆明均屬于金錢松的低適生區(qū),而Domain模型中這一區(qū)域則表現(xiàn)為高適區(qū),經(jīng)查閱文獻發(fā)現(xiàn)河南鄭州、信陽大別山及云南昆明等地區(qū)均有人工成功引種金錢松的記錄[28- 29],可確認(rèn)是金錢松適生區(qū)。

3.2 金錢松分布范圍變化

金錢松為高大落葉喬木,為深根系樹種[5],其一年生幼苗生長旺盛期為4—5月和8—9月,這兩個階段的月均溫和月降水量接近,為金錢松幼苗生長的最適條件[29]。本文發(fā)現(xiàn)年均溫(Bio1)與年降水量(Bio12)是金錢松潛在適生區(qū)預(yù)測的主要影響因子,也決定了金錢松分布格局的形成以及遷移方向。末次冰期氣候條件惡劣[23,30],而在金錢松分布區(qū)的浙江天目山、大明山、安徽黃山以及湖南、江西、湖北三省交界處的羅霄山脈是其良好的避難所[31],這些地區(qū)山峰地形多變,小生境水熱條件優(yōu)越,有效的降低了冰期氣候?qū)χ参锓植嫉挠绊?。從末次間冰期(LGM)到當(dāng)前氣候條件下金錢松預(yù)測分布范圍變化可知,因冰期之后氣溫回暖,金錢松的分布格局逐漸向溫?zé)岬哪戏綌U張。從當(dāng)前到2070年金錢松的預(yù)測分布范圍變化可知,未來其分布會一定程度向北擴張,但在南方卻大面積收縮,總縮小面積達到178084.30km2。隨著全球氣候變化加劇[32],溫室效應(yīng)將大大影響植物分布格局以及其物候期[33-35],而植物分布變化是個較為緩慢的過程,全球氣溫上升速度遠(yuǎn)遠(yuǎn)超過金錢松自身遷移的能力,因此必須適時考慮人為輔助金錢松的北向遷移[33]。

3.3 金錢松潛在適生區(qū)與金錢松保護

金錢松作為瀕危物種,自然分布范圍有限,且生境受到各種威脅,建立金錢松自然保護區(qū)與種子園可以有效保護金錢松[5],自然保護區(qū)可以最大程度防止人為破壞,種子園可以有效保障種質(zhì)資源和遺傳水平,而二者的建立與物種適宜生境密切相關(guān)[36- 38]。本文結(jié)合多種生態(tài)位模型對金錢松的潛在適生區(qū)進行預(yù)測及其分布格局變化分析,對金錢松保護提出以下建議：除現(xiàn)有浙江西天目山金錢松自然保護區(qū)外,還可在江西銅鼓縣、湖南張家界和衡陽這幾處當(dāng)前金錢松高適分布區(qū)域內(nèi)建立自然保護區(qū)與種子園,保護與育苗、培種相結(jié)合,最大程度保護現(xiàn)有種群；在目前金錢松中等適宜區(qū)以及未來氣候條件下高適分布區(qū),如安徽北部、河南南部、湖北東南部等地區(qū)可以考慮人工引種,進行原土栽植協(xié)助金錢松北移實驗,未來可將其作為園林樹種,發(fā)展為園林綠植,而4種模型都預(yù)測為低適區(qū)或不適區(qū)的地方則不適于選擇種植金錢松。