基于MaxEnt模型的長(zhǎng)足大竹象在中國(guó)潛在分布區(qū)及其對(duì)氣候變化的響應(yīng)

楊春平，趙 霞，王嘉雯，楊 樺*，楊 偉，吳 南

（1.四川農(nóng)業(yè)大學(xué)林學(xué)院，成都 611130；2.四川省林業(yè)有害生物防治協(xié)會(huì)，成都 610081）

長(zhǎng)足大竹象（Cyrtotrachelus buqueti），屬鞘翅目Coleoptera象蟲總科Curculionoidea象甲科 Curculionidae彎頸象屬Cyrtotrachelus昆蟲。主要分布于我國(guó)四川、貴州、重慶和廣西等省市及越南、緬甸、泰國(guó)等東南亞國(guó)家分布[1-2]。主要危害慈竹（Neosinocalamus affinis）、青皮竹（Bambusa textilis）、綠竹（Bambusa oldhamii）等叢生竹幼筍[3]。長(zhǎng)足大竹象一年一代，以成蟲在地下越冬，翌年7月中旬成蟲出土，8月中旬為出土盛期，9月下旬成蟲在地上活動(dòng)終止。目前，長(zhǎng)足大竹象已成為危害我國(guó)叢生竹資源的重要害蟲，給竹產(chǎn)業(yè)帶來了嚴(yán)重影響，僅四川省每年發(fā)生面積近67 000 hm2，為害率50%～80%，嚴(yán)重可達(dá)100%，已成為制約造紙用竹林發(fā)展的重要因素[4]。

全球變化對(duì)全球各個(gè)方面有著嚴(yán)重影響，其中以全球氣候變暖最為顯著。氣候變暖對(duì)生態(tài)系統(tǒng)結(jié)構(gòu)和功能、群落組成以及物種分布等均具有影響[5-6]。氣候變化直接影響昆蟲的行為、發(fā)育、分布、生殖和生存[7-8]，對(duì)農(nóng)林業(yè)生產(chǎn)和人類健康具有潛在的重要影響[9-10]。物種分布模型（species distribution models）是將物種分布的樣本信息以及與其相對(duì)應(yīng)的環(huán)境變量進(jìn)行關(guān)聯(lián)比對(duì)，得出兩者之間的關(guān)系，以此對(duì)物種未來的潛在分布區(qū)域進(jìn)行預(yù)測(cè)[11]。近年來，國(guó)內(nèi)外常用的有最大熵模型MaxEnt、GARP模型和生物氣候模型BIOCLIM等14種物種分布模型[12-13]。由于MaxEnt模型運(yùn)行時(shí)間較短、結(jié)果穩(wěn)定、所需的樣本量小，同時(shí)其預(yù)測(cè)精確度明顯優(yōu)于其他模型，其運(yùn)用也是最為廣泛的[14]，如斑翅果蠅（Drosophila suzukii）[15]、水稻象甲（Echinocnemus squamous）[16]、華山松大小蠹（Dendroctonus armandi）[17-18]、懸鈴木方翅網(wǎng)蝽（Corythucha ciliate）[19]、和西藏飛蝗（Locusta migratoria tibetensis）[20]等昆蟲在國(guó)內(nèi)的潛在分布區(qū)。

目前，關(guān)于長(zhǎng)足大竹象的研究主要集中于生物學(xué)、行為學(xué)[21-24]，對(duì)其在大尺度范圍下的潛在分布區(qū)預(yù)測(cè)和氣候變化的響應(yīng)機(jī)制的研究較少[25-26]。由于長(zhǎng)足大竹象寄主的適生條件不同，本研究運(yùn)用MaxEnt模型，結(jié)合GIS技術(shù)，預(yù)測(cè)當(dāng)前及未來時(shí)段不同氣候情景下長(zhǎng)足大竹象在我國(guó)的潛在分布區(qū)，探究環(huán)境變量對(duì)長(zhǎng)足大竹象分布的影響，以及對(duì)長(zhǎng)足大竹象蟲害的發(fā)生預(yù)警和制定防治方案提供理論依據(jù)。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)來源

1.1.1 預(yù)測(cè)軟件來源

MaxEnt模型軟件預(yù)測(cè)長(zhǎng)足大竹象在中國(guó)的適生區(qū)，該軟件可在MaxEnt（http://www.cs.princeton.edu/-schapire/maxent/）上獲取[27]。利用 ArcGIS 軟件處理MaxEnt預(yù)測(cè)結(jié)果繪制長(zhǎng)足大竹象實(shí)際分布區(qū)及潛在分布區(qū)。

1.1.2 物種分布記錄

研究中2012—2016年四川省內(nèi)長(zhǎng)足大竹象分布數(shù)據(jù)通過查詢四川省林業(yè)有害生物綜合信息系統(tǒng)（sc.30120.org/fpgis/secure/Login.aspx?ReturnUrl=%2ffpgis%2f）獲得，四川省外數(shù)據(jù)通過查詢物種分布數(shù)據(jù)庫(kù)和國(guó)內(nèi)外公開發(fā)表的長(zhǎng)足大竹象相關(guān)的期刊論文獲得。本研究查詢的數(shù)據(jù)庫(kù)包括國(guó)際農(nóng)業(yè)與生物科學(xué)中心數(shù)據(jù)庫(kù)（CABI，http://www.cabi.org/）、教學(xué)標(biāo)本資源共享平臺(tái)（http://mnh.scu.edu.cn/）。除去重復(fù)、模糊的分布點(diǎn)，查詢各分布點(diǎn)的經(jīng)緯度，共獲得長(zhǎng)足大竹象分布點(diǎn)346個(gè)[20]。

1.1.3 環(huán)境變量的選擇

環(huán)境數(shù)據(jù)包括23個(gè)變量（19個(gè)生物氣候變量和月平均氣候數(shù)據(jù)），主要包括降雨和溫度，上述環(huán)境數(shù)據(jù)來自Worldclim氣候數(shù)據(jù)庫(kù)（http://www.worldclim.org/version2），空間分辨為2.5arc-minutes（約4.5 km2）（表1）。下載的數(shù)據(jù)首先利用ArcGIS軟件中空間分析工具，按中國(guó)輪廓進(jìn)行掩膜裁剪，并在Arc GIS中轉(zhuǎn)化為ASC格式。

未來時(shí)段包括2050s（2041—2060年）和2080s（2071—2090年），未來氣候情景從氣候?qū)ｉT委員會(huì)中第五次報(bào)告發(fā)布的4種濃度路徑中選擇的3種氣候情景[26]，包括 RCP2.6、RCP4.5 和 RCP8.5，在國(guó)際熱帶農(nóng)業(yè)中心（International Centre for Tropical Agriculture，CIAT）網(wǎng)站下載該數(shù)據(jù)[28]。

1.2 研究方法

1.2.1 環(huán)境變量的預(yù)處理

為確保模型的模擬精度，分析19個(gè)初始環(huán)境變量包括預(yù)測(cè)物種分布通用的19個(gè)生物氣候變量和月平均氣候數(shù)據(jù)之間的相關(guān)性，去掉高度相關(guān)但對(duì)MaxEnt模型預(yù)測(cè)結(jié)果貢獻(xiàn)不大的變量[29]。本次研究參考Worthington等[30]方法進(jìn)行潛在環(huán)境變量的篩選。

表1 生物氣候特征一覽表Table 1 List of bioclimatic characteristics

1.2.2 長(zhǎng)足大竹象的模型構(gòu)建與評(píng)價(jià)

在MaxEnt軟件中進(jìn)行模型構(gòu)建[31]，從分布點(diǎn)中隨機(jī)選取75%為訓(xùn)練數(shù)據(jù)集，25%為驗(yàn)證數(shù)據(jù)集，選擇刀切法（Jackknife），其他參數(shù)為默認(rèn)參數(shù)，運(yùn)行該模型。

本次研究模型模擬結(jié)果由受試者工作特征曲線（receiver operating characteristic，ROC曲線）下面積值（AUC值）來進(jìn)行評(píng)價(jià)[32]。AUC值的取值范圍為［0.5，1］，值越接近1表明預(yù)測(cè)結(jié)果精度越高（表 2）。

1.2.3 長(zhǎng)足大竹象適生等級(jí)劃分

長(zhǎng)足大竹象潛在分布區(qū)的等級(jí)劃分在ArcGIS中進(jìn)行，根據(jù)王茹琳等[20]與劉艷梅等[33]的劃分方法:存在概率＜0.05為不適生區(qū)；0.05≤存在概率＜0.33為低適生區(qū)；0.33≤存在概率＜0.66為中適生區(qū)；存在概率≥0.66為高適生區(qū)。

表2 AUC評(píng)價(jià)標(biāo)準(zhǔn)表Table 2 The evaluation criterion of AUC

1.2.4 幾何中心及位移計(jì)算

參考Yue T.等[34]的計(jì)算方法，計(jì)算在未來不同時(shí)期高適生區(qū)的面積變化情況和中心點(diǎn)位移，具體計(jì)算公式如下:

2 結(jié)果及分析

2.1 MaxEnt模型模擬評(píng)價(jià)

圖1是當(dāng)前時(shí)段條件下的ROC曲線圖，圖中表明基于主導(dǎo)環(huán)境變量構(gòu)建的長(zhǎng)足大竹象在中國(guó)的地理分布模型的AUC值為0.985，預(yù)測(cè)準(zhǔn)確性達(dá)到極好的標(biāo)準(zhǔn)；圖2為未來不同氣候情境下模擬結(jié)果的ROC曲線圖，由圖2可知所有預(yù)測(cè)結(jié)果的AUC值都達(dá)到極好的標(biāo)準(zhǔn)。上述結(jié)果說明模型預(yù)測(cè)可信度較高，可用于分析氣候變化對(duì)長(zhǎng)足大竹象在中國(guó)分布的影響。

2.2 主要環(huán)境變量對(duì)模型預(yù)測(cè)長(zhǎng)足大竹象分布的貢獻(xiàn)率

經(jīng)過上述環(huán)境變量的篩選后，共保留了6個(gè)主要環(huán)境變量，重建長(zhǎng)足大竹象的Maxent模型。刀切法檢驗(yàn)結(jié)果顯示，對(duì)長(zhǎng)足大竹象分布貢獻(xiàn)較大的環(huán)境變量有:二月最低溫（tmin2）、最暖季降雨量（bio18）、7 月降雨量（prec7）等（圖 3、圖 4），tmin2貢獻(xiàn)率大。

2.3 長(zhǎng)足大竹象在中國(guó)的適生分布

研究預(yù)測(cè)結(jié)果表明（表3），當(dāng)前長(zhǎng)足大竹象在中國(guó)的適生區(qū)范圍為 22°～33°N，104°～122°E，位于我國(guó)亞熱帶地區(qū)，由四川東南部向我國(guó)東南部蔓延，總面積達(dá)87 983 km2，其中高適生區(qū)主要位于四川東南部，在貴州、重慶和云南等省份均有分布，面積達(dá)8 607 km2，中適生區(qū)主要分布于四川、貴州、廣西、重慶等省份，面積為15 064 km2；低適生區(qū)主要分布于廣西、廣東、湖南、湖北、浙江和云南等省份，面積為64 312 km2。

預(yù)測(cè)結(jié)果同時(shí)表明，3種氣候情景下高適生區(qū)仍主要分布于四川、貴州和重慶一帶，隨時(shí)間推移，氣候變化加劇，長(zhǎng)足大竹象高適生區(qū)向廣西一帶偏移，CO2排放越低，害蟲適生區(qū)越廣，高適生區(qū)面積也最大。

2.4 長(zhǎng)足大竹象在中國(guó)適生區(qū)面積變化及中心點(diǎn)位移

當(dāng)前及RCP2.6情景下，總適生面積和高適生區(qū)面積均先增加后減少，總適生區(qū)由87 983 km2增加至92 799 km2（2050s）后降低為89 883 km2（2080s）。高適生區(qū)面積由當(dāng)前的8 607 km2增加至的12 778 km2（2080s）；當(dāng)前至 2080 s，高適生區(qū)中心點(diǎn)由重慶榮昌（當(dāng)前）經(jīng)重慶江津（2050）位移至貴州習(xí)水（2080s）（圖5），總體上向東南方向移動(dòng)12.66 km（表4）；總適生區(qū)中心點(diǎn)由重慶綦江（當(dāng)前）經(jīng)貴州石阡（2050s）位移至貴州印江（2080s）（圖 6），總體上向東南方向移動(dòng)20.71 km（表5）。

圖1 當(dāng)前氣候情景下Maxent模型預(yù)測(cè)結(jié)果ROC曲線Figure 1 ROC curve verification of Maxent prediction under the current climate scenario

圖2 未來3種氣候情景下長(zhǎng)足大竹象Maxent模型預(yù)測(cè)的ROC曲線Figure 2 ROC curves predicted of C.buqueti by Maxent model for the future three climate scenarios

圖3 氣候因子對(duì)長(zhǎng)足大竹象的貢獻(xiàn)值Figure 3 The contribution value of climate factors to C.buqueti

圖4 單變量響應(yīng)曲線Figure 4 Single variable response curve

表3 長(zhǎng)足大竹象在當(dāng)前情景及未來氣候情景的預(yù)測(cè)面積Table 3 Predicted suitable areas for C.buqueti in current and future climatic scenarios

當(dāng)前及RCP4.5情景下，總適生面積呈先增加后減少的變化趨勢(shì)，總體上由當(dāng)前的87 983 km2增加至2070s的81 137 km2；高適生區(qū)面積持續(xù)增加，由8 607 km2增加至2070s的10 663 km2。當(dāng)前至2080s，高適生區(qū)中心點(diǎn)由重慶榮昌（當(dāng)前）經(jīng)重慶大足（2050s）位移至重慶巴南（2080s）（圖 5），總體上向東南方向移動(dòng)14.44 km（表4）；總適生區(qū)中心點(diǎn)由重慶綦江（當(dāng)前）經(jīng)貴州德江（2050s）位移至重慶秀山（2080s）（圖6），總體上向東南方向移動(dòng)23.28 km（表 5）。

當(dāng)前及RCP8.5情景下，高適生區(qū)面積呈持續(xù)增加變化趨勢(shì)，總體上升由當(dāng)前的8 607 km2增加至2080s的11 593 km2；總適生面積先增加后減少，由87 983 km2增加至2070s的89 024 km2。當(dāng)前至2080s，高適生區(qū)中心點(diǎn)由重慶榮昌（當(dāng)前）經(jīng)四川合江（2050s）位移至貴州遵義（2070s）（圖 5），總體上向東南方向移動(dòng)26.03 km（表4）；總適生區(qū)中心點(diǎn)由重慶綦江（當(dāng)前）經(jīng)貴州鳳岡（2050s）位移至貴州石阡（2080s）（圖6），總體上向南方向移動(dòng)19.71 km（表 5）。

3 討論與結(jié)論

MaxEnt是現(xiàn)階段物種潛在區(qū)分布使用最為廣泛的生態(tài)位模型，現(xiàn)已廣泛使用于植物及昆蟲的預(yù)測(cè)，且預(yù)測(cè)結(jié)果較其他生態(tài)位模型更準(zhǔn)確。該軟件基于最大熵原理，以物種分布變量和環(huán)境變量為基礎(chǔ)，利用數(shù)學(xué)模型統(tǒng)計(jì)分析熵最大時(shí)物種的分布狀態(tài)[20]。因此，選用MaxEnt進(jìn)行本次研究，預(yù)測(cè)長(zhǎng)足大竹象在中國(guó)適生區(qū)分布。

圖5 未來3種氣候情景下長(zhǎng)足大竹象高適生區(qū)中心點(diǎn)遷移軌跡Figure 5 The migration trajectory of the center point of highly suitable area of C.buqueti in three future climatic scenarios

表4 高適生區(qū)中心點(diǎn)變化距離及方向Table 4 Shift distance and direction of center of highly suitable area

圖6 未來3種氣候情景下長(zhǎng)足大竹象總適生區(qū)中心點(diǎn)遷移軌跡Figure 6 The migration trajectory of the center point of totally suitable area of C.buqueti in three future climatic scenarios

表5 總適生區(qū)中心點(diǎn)變化距離及方向Table 5 Shift distance and direction of center of totally suitable area

MaxEnt預(yù)測(cè)結(jié)果表明，長(zhǎng)足大竹象當(dāng)前在我國(guó)主要適生區(qū)分布于我國(guó)西南部及南部，22°～33°N，104°～122°E，總面積達(dá) 87 983 km2，其中高適生區(qū)主要位于四川東南部、貴州以及重慶，中適生區(qū)主要分布于四川、貴州、廣西、重慶等省份，與李濤[35]等報(bào)告中我國(guó)長(zhǎng)足大竹象分布區(qū)一致，新增湖北、福建、陜西等中低適生區(qū)，長(zhǎng)足大竹象潛在分布符合我國(guó)叢生竹分布范圍[36-37]。但在湖南和云南兩省竹林廣泛分布的省區(qū)僅有少部分低適生區(qū)，原因可能與當(dāng)?shù)馗蔁釟夂蛴嘘P(guān)，湖南春季氣溫變化大，夏秋兩季常干旱；嚴(yán)寒期短，暑熱期長(zhǎng)，加上湖南亞熱帶季風(fēng)氣候?qū)俅箨懶裕瑢?dǎo)致高溫時(shí)間持續(xù)長(zhǎng)[38]。由于長(zhǎng)足大竹象在高溫條件下繁殖慢，繁殖期短，在湖南云南兩省不利于繁殖；導(dǎo)致長(zhǎng)足大竹象在兩地不適生。梁祖昂[24]在研究中表明冬暖夏濕，為大竹象提供有利生存條件，冬暖有利于成蟲越冬，夏濕有利卵孵化，這與本研究中湖南和云南不適生條件相符。刀切法檢驗(yàn)結(jié)果顯示，對(duì)長(zhǎng)足大竹象分布貢獻(xiàn)最大的環(huán)境變量是二月最低溫（tmin2）和最暖季降雨量（bio18），其次是 7月降雨量（prec7），這可能與長(zhǎng)足大竹象的越冬期和出土期有關(guān)。長(zhǎng)足大竹象以成蟲在土中越冬，5月下旬出土[39]，冬季溫度偏高，有利于長(zhǎng)足大竹象越冬，越冬死亡率降低，發(fā)生危害加重[25]。7月下旬為出土高峰，降雨的物理沖刷作用使土壤阻力降低，有利于成蟲出土，也確保叢生竹的生長(zhǎng)為長(zhǎng)足大竹象提供食物及生存空間[38-39]。

在不同氣候情景下，隨時(shí)間的推移，長(zhǎng)足大竹象總適生面積均呈先上升后下降趨勢(shì)，2050s適生面積最大；高適生區(qū)面積呈上升趨勢(shì)。但不同時(shí)段及不同氣候情景下差異明顯，這可能與不同排放情景下所造成的溫度及降雨變化有關(guān)，在RCP2.6情景下形成利于長(zhǎng)足大竹象繁殖的溫度與降雨導(dǎo)致其中、高適生區(qū)面積變化最為明顯，2050s高適生區(qū)從四川重慶大部分高適生區(qū)增加到廣西貴州交界處高適生，2080s高適生區(qū)蔓延至貴州西部并像廣州中部進(jìn)行蔓延。其中，廣西地區(qū)蔓延最為明顯，由當(dāng)前大部分低適生區(qū)明顯變化為中適生區(qū)并覆蓋廣西全省。溫度是影響長(zhǎng)足大竹象生長(zhǎng)發(fā)育的重要環(huán)境因子[26]，而CO2排放的多少直接對(duì)溫度造成影響，中高情景下，CO2排放量高，溫度更高，易造成干熱天氣，不利于長(zhǎng)足大竹象繁殖。

氣候是決定物種分布的最主要因素，氣候的變化將直接影響昆蟲的生長(zhǎng)發(fā)育、繁殖、存活、寄主植物等核心生命活動(dòng)[9-10]。本研究以面積作為權(quán)重，計(jì)算不同等級(jí)適生區(qū)中心點(diǎn)位置，分析高適生區(qū)、總適生區(qū)中心點(diǎn)隨時(shí)間變化的位移趨勢(shì)，展示長(zhǎng)足大竹象對(duì)氣候變化的響應(yīng)過程[33]。對(duì)不同氣候變化情景下長(zhǎng)足大竹象適生區(qū)中心點(diǎn)的遷移變化分析表明，高適生區(qū)和總適生區(qū)在中國(guó)的中心點(diǎn)變化位于重慶、貴州和四川，總體上東南方向發(fā)生偏移。高適生區(qū)中心點(diǎn)變化一致，均向東南方向發(fā)生偏移，總適生區(qū)中心點(diǎn)位移在RCP4.5和RCP2.6下先向東北后向東南偏移，RCP8.5下向東南方向偏移。孫玉誠(chéng)等[9，40]得到的全球氣候變暖會(huì)引起昆蟲向高海拔、高緯度地區(qū)遷移的結(jié)論不完全吻合，這可能與不同物種對(duì)氣候的要求及研究選用的分布點(diǎn)不同有關(guān)。

本研究精確度以AUC值保證，由于發(fā)病率和診斷閾值對(duì)AUC值無影響，還可對(duì)不同試驗(yàn)的準(zhǔn)確度進(jìn)行綜合比較，因此AUC值是目前學(xué)者公認(rèn)的評(píng)價(jià)診斷試驗(yàn)結(jié)果的最佳指標(biāo)[41]。本研究結(jié)果表明，當(dāng)前及未來氣候模式下長(zhǎng)足大竹象的適生區(qū)預(yù)測(cè)的AUC值均大于0.95，預(yù)測(cè)準(zhǔn)確性達(dá)到極好標(biāo)準(zhǔn)，說明此次模型預(yù)測(cè)的地理分布與長(zhǎng)足大竹象實(shí)際分布擬合度較高。由于物種分布所受的限制因子多樣，不僅僅受氣候影響，預(yù)測(cè)結(jié)果有一定偏差，使用分布數(shù)據(jù)的多少也會(huì)影響模型的精度，在下一步研究中將加以改進(jìn)，提高模型預(yù)測(cè)精確度。