中國跳鼠總科物種系統(tǒng)發(fā)育與分布格局

隋璐璐 程繼龍 夏霖 楊奇森*

（1 中國科學院動物研究所，北京 100101）（2 中國科學院大學，北京 100049）

隨著物種滅絕速率的不斷加快，生物多樣性研究和保護工作正成為生物學研究的重中之重(Daviset al., 2018)，正確認識物種并明確物種的分布狀況是生物多樣性研究和保護的重要前提。但是，僅依賴傳統(tǒng)形態(tài)分類學方法進行的物種分類已不利于對隱性物種多樣性的發(fā)掘(Sheth and Thaker,2017)，可能導致生物多樣性被低估。隨著分子生物學技術的發(fā)展和系統(tǒng)分類學、演化生物學研究的不斷深入，DNA 序列變異這一客觀特征對物種鑒定和物種間系統(tǒng)發(fā)育關系的構建起到了極大的促進作用(Yang,2006)。同時，厘清物種間的系統(tǒng)發(fā)育關系，也有助于更深入認識生物多樣性的內涵，為物種的保護工作提供理論基礎(Pérez-Losadaaet al.,2002;Woodet al.,2019)。

物種分布模型(Species distribution model,SDM)能利用物種分布信息和當地氣候變量之間的關系預測物種潛在分布區(qū)(Guisan and Zimmer‐mann, 2000; Ma and Sun, 2018; Mwakapejeet al.,2019; Zhanget al., 2019)。中國的國土面積遼闊，生物學調查難以全部覆蓋，因此，SDM 可以在一定程度上指導未來的野外調查方向，為進一步探索和發(fā)現物種分布提供參考，為物種保護措施的制定提供科學基礎(李娜等，2020)。在多種分布模型的模擬方法中，最大熵模型(Maximum entropy model,Maxent)可以在樣本量較小的情況下，利用環(huán)境變量和物種分布位點來計算約束條件，并探索在此約束條件下最大熵的可能分布，進而預測研究區(qū)物種的生境適宜性(Merowet al.,2013)，得出優(yōu)于其他模型的模擬結果(Phillipset al.,2006)。

跳鼠，因其高度特化的后肢結構和二足跳躍的運動方式而得名，廣泛分布于我國中西部的諸多地區(qū)，跨越森林、草原和荒漠3種不同的典型生境，是我國生物多樣性尤其是荒漠地區(qū)生物多樣性的重要指示類群。同時，我國還是跳鼠物種的起源、擴散和分布中心(Zhanget al., 2013; Pisanoet al., 2015)。程繼龍等(2021) 對中國跳鼠總科物種的分類系統(tǒng)進行了綜述，認為我國共分布有跳鼠總科動物3 科6 亞科12 屬22 種，是世界上跳鼠總科物種豐富度最高的國家之一。以前的研究中，無論是從分子角度(Zhanget al., 2013; Pisanoet al., 2015) 或是形態(tài)結構角度(Wuet al., 2014;Mooreet al., 2015)，其系統(tǒng)發(fā)育關系中都未能包括我國分布的全部跳鼠，尤其是最近幾年新發(fā)現和修訂的類群，如奇美跳鼠(Chimaerodipus auri‐tus) (Shenbrotet al., 2017) 和塔里木跳鼠(Dipus deasyi) (Chenget al., 2018)。因此，我國跳鼠總科新的系統(tǒng)發(fā)育關系亟待完善。同時，各個物種的分布狀況也缺乏全面調查和整理，難以滿足對該類群本底資料的掌握和認識生物多樣性的需要。

本研究以中國分布的跳鼠總科動物為研究對象，通過研究物種DNA 序列上的變異，厘清物種的系統(tǒng)發(fā)育關系，探討我國跳鼠總科內部在科、屬和種3個分類階元上的演化關系，為我國跳鼠總科分類提供理論基礎。同時，通過物種分布模型，推測跳鼠總科物種在我國的潛在分布范圍，并繪制物種豐富度格局，分析不同物種分布狀況的主要影響因素，為我國跳鼠總科多樣性的研究工作提供參考。

1 研究方法

1.1 系統(tǒng)發(fā)育關系重建

1.1.1 資料收集

從GenBank 數據庫中下載中國分布的跳鼠總科全部22 個物種(表1)12 個特定基因的全部序列，包括2 個線粒體基因：細胞色素氧化酶基因亞基I(cytochrome oxidase subunit I,COI) 和細胞色素b(cytochrome b,Cyt b)；10個核基因：α-腎上腺受體2B 基 因(alpha adrenergic receptor 2B,A2AB)、α 肽銅轉運ATP 酶(ATPase Copper Transporting Alpha,ATP7A)、乳腺癌蛋白1 基因(breast cancer protein 1,BRCA1)、乳腺癌蛋白2 基因(breast cancer pro‐tein 2,BRCA2)、大麻素受體1 基因(cannabinoid re‐ceptor 1,CNR1)、生長素受體基因外顯子10 (exon 10 of the growth hormone receptor,GHR)、光受體間視黃酸結合蛋白基因外顯子1(exon 1 of interphoto‐receptor retinoid-binding protein,IRBP)、部分重組啟動蛋白1 基因(a portion of the recombination-acti‐vating gene 1,RAG1)、重組啟動蛋白2基因(recom‐bination-activating gene2,RAG2)和血管性血友病因子基因(von Willebrand factor,vWF)。同時，選擇跳鼠總科8 個近緣物種作為外群(表1)，利用MEGA-X(Kumaret al.,2018)進行篩選，排除不確定性序列和測序質量差的序列，并盡可能降低亞種干擾，選擇物種分布中心區(qū)域的序列，最終挑選出每個基因中最具代表性的基因序列。

表1 本研究應用的物種Table 1 The species used in this study

1.1.2 系統(tǒng)發(fā)育分析

利用MEGA-X 中的Muscle (Edgar, 2004) 進行序列的比對，將序列兩端起始的部分堿基刪除后導出保存。再利用Sequence Matrix 1.8(Vaidyaet al.,2011)串聯(lián)基因，構建基因矩陣?；谒@得的基因矩陣，分別采用最大似然法(maximum likelihood,ML)和貝葉斯法(Bayesian inferences,BI)進行系統(tǒng)發(fā)育分析。運用Phylosuite(Zhanget al.,2020)軟件中的PartitionFinder2 (Lanfearet al., 2016)，基于“貪婪算法”(greedy search algorithm)和貝葉斯信息準則(Bayesian information criterion,BIC) 選擇最優(yōu)分區(qū)方式和最佳核酸替換模型。使用IQ-TREE(Guindonet al., 2010;Anisimovaet al., 2011; Minhet al.,2013;Nguyenet al.,2015)構建ML 樹，執(zhí)行5 000次超快自舉(ultrafast bootstraping)估算結點支持率。使用MrBayes 3.2.6 (Ronquistet al., 2012)構建BI樹，運行500萬代，每1 000代抽樣一次，舍棄最初的25%抽樣。最后將得到的ML樹和BI樹通過FigTree 1.4.4修飾并導出。

1.2 分布格局計算

1.2.1 數據收集及預處理

檢索和收集中國分布的跳鼠總科物種的分布數據，數據源包括全球生物多樣性信息網絡數據庫(Global Biodiversity Information Facility, GBIF)，中國科學院動物研究所國家動物博物館(National Zoological Museum of China, NZMC)，以及書籍或文獻中可以查到的分布和采集記錄(王思博和楊贛源，1983；張榮祖，1997；阿不力米提·阿布都卡迪爾，2003；Quereet al., 2009; Pisanoet al.,2015; Shenbrotet al., 2017; Chenget al., 2018,2021)。依據最新研究結果(Michaux and Shenbrot,2017;Shenbrotet al.,2017;Chenget al.,2018)，篩選出符合中國跳鼠總科22個物種劃定的分布記錄，刪除空白值、偏移值和冗余數據，最終獲得603個分布位點。其中奇美跳鼠僅包含5 個分布位點，數據較少，不能滿足物種分布模型模擬要求，暫時排除不用。

1.2.2 環(huán)境因子選擇

依據分布數據特點和跳鼠總科動物生活習性，分別從全球氣候數據庫(WorldClim2, www. world‐clim. org) 和地理空間數據云數據庫(www. gs‐cloud.cn)下載溫度、降水、氣象和地形4 種變量，并轉換為分辨率為30 s的環(huán)境數據變量；考慮到跳鼠總科動物冬眠的習性，對于太陽輻射(sard)、水蒸氣壓(vapr)和風速(wind)3 個氣象數據變量只選取其活動季節(jié)(4—9 月)數據。最終獲得23 個環(huán)境變量，共38個柵格圖層(表2)。

表2 跳鼠總科物種分布模型使用的環(huán)境變量Table 2 Environment variables used for species distribution model of Dipodoidea

1.2.3 Maxent模型選擇

部分研究表明，Maxent 的默認參數設置是基于最大經驗調優(yōu)的結果，可能導致模型性能較差(Radosavljevicet al., 2013; Ma and Sun, 2018)。因此我們結合物種分布位點和篩選出來的環(huán)境變量，利用R 語言(R 4.0.5)的ENMeval 包(Muscarellaet al.,2014)對Maxent 進行模型優(yōu)化。設置最大背景點數為10000，要素類型(feature combination, FC)分別為L、LQ、H、LQH、LQHP 和LQHPT (L =linear, Q = quadratic, H = hinge, P= product and T =threshold)，調控頻率(regularization multiplier,RM)從0 至4，每次遞增0.5。依次匹配FC 和RM 進行模型計算，篩選出赤池信息標準(AICc)值最小的參數設置，作為最優(yōu)參數進行物種分布模型計算。

應用Maxent 3.4.1 軟件，設置隨機檢驗百分比為25%，模型重復運行30 次，重復運行類型為二次抽樣(subsample)，并使用刀切法(jackknife test)分析環(huán)境變量的重要性。以ROC 曲線(receiv‐er operating characteristic curve) 對模型的準確度進行評估。ROC 曲線與橫軸圍成的面積為AUC 值(Hanley and McNeil, 1982)，可用來衡量模型的優(yōu)劣。一般而言，AUC 值為0~0.6 表示模型預測無效；0.6 ~ 0.7 表示模型準確度較差；0.7 ~ 0.8 表示模型準確度一般；0.8 ~ 0.9 表示模型準確度良好；大于0.9表示模型準確度極好。

1.2.4 物種適宜度和豐富度

由于每個物種對環(huán)境的耐受程度不一樣，基于可查到的分布記錄的臨界值，對每個物種的適宜棲息地閾值進行劃分(Svenninget al.,2008)。然后依據閾值提取每個物種適宜棲息地圖層，疊加后得到物種適宜分布圖。對物種分布預測結果的圖層進行二值化處理，即每個柵格中該物種適宜度大于閾值則設置為1，否則為0。然后將物種二值分布圖進行空間疊加，獲得中國跳鼠總科物種豐富度分布格局。

2 結果

2.1 系統(tǒng)發(fā)育關系

構建系統(tǒng)發(fā)育關系的12 個基因串聯(lián)矩陣總長11 340 bp，包含中國分布的跳鼠總科全部22 個物種，每個物種至少4個基因片段。貝葉斯樹與最大似然樹拓撲結構一致(圖1)，從樹的根部到末端，僅巨泡五趾跳鼠(Orientallactaga bullata)和巴里坤跳鼠(O. balikunica)分支結點(BP=45,PP=0.53)的支持率較低，因此支持22 個跳鼠總科物種的分類體系。五趾跳鼠亞科(Allactaginae) 中，小五趾跳鼠(Scarturus elater)與東方五趾跳鼠屬(Oriental‐lactaga) 聚為一支，自舉支持率相對較低(BP =60)，但后驗概率較高(PP = 0.92)；大五趾跳鼠(Allactaga major) 與小地兔(Pygeretmus pumilio) 互為近緣類群(BP=80,PP=1.00)。跳鼠亞科(Dipo‐dinae)中，各結點的支持率都非常高，其中奇美跳鼠與羽尾跳鼠屬(Stylodipus) 聚為一支(BP = 100,PP=1.00)。跳鼠科(Dipodidae)內部，跳鼠亞科和五趾跳鼠亞科關系更近(BP= 100, PP= 1.00)，心顱跳鼠亞科(Cardiocraniinae) 最為原始(BP= 100,PP=1.00)。跳鼠總科內的各亞科之間系統(tǒng)發(fā)育關系明確，跳鼠亞科、五趾跳鼠亞科、長耳跳鼠亞科(Euchoreutinae)、心顱跳鼠亞科聚為跳鼠科，跳鼠科與林跳鼠科(Zapodidae)的親緣關系更近，蹶鼠科(Sicistidae)位于跳鼠總科最基部。

圖1 基于12個基因的串聯(lián)矩陣構建的跳鼠總科ML樹. 結點上的數值分別表示自舉支持率和Bayesian后驗概率Fig. 1 ML tree based on gene matrix of 12 genes. Number at the nodes correspond to bootstrap values and Bayesian posterior probabilities

2.2 物種豐富度格局

基于分布位點得到的物種豐富度結果表明，跳鼠總科主要分布在我國北方的荒漠、半荒漠地區(qū)，青藏高原東北緣的山地與高山草原地區(qū)。物種豐富度最高的分布區(qū)位于新疆北部包括昌吉州、阿勒泰地區(qū)、塔城地區(qū)、博爾塔拉州和伊犁等地，一般分布有7 種以上跳鼠總科動物。其中，塔城、額敏、裕民和托里部分地區(qū)甚至包含9~10種跳鼠總科動物(圖2)。物種豐富度次高一級的分布區(qū)集中在新疆西南部、內蒙古中西部和甘肅省河西走廊地區(qū)，分布有4種以上跳鼠總科動物。其他區(qū)域的物種豐富度較低，一般只有2種以下。

圖2 中國跳鼠總科物種豐富度分布格局Fig. 2 Species richness of Dipodoidea in China

2.3 物種適宜度及主要影響因素

跳鼠總科在我國適宜分布區(qū)域總面積約2.22×106km2，約占國土總面積的23%。其中，蹶鼠科適宜分布面積14.15×104km2(圖3A)；林跳鼠科適宜分布面積11.85×104km2(圖3B)；跳鼠科適宜分布面積達20.423×105km2(圖3C)。21 個物種分布模型ROC 曲線的AUC 值均高于0.889，說明Maxent模擬跳鼠總科物種的潛在分布范圍的模型擬合度相對較高，結果可信。環(huán)境變量對物種分布模型預測的貢獻百分比結果表明，影響不同物種分布的主要環(huán)境因子不盡相同(表3)。年均溫是四川林跳鼠(Eozapus setchuanus)、中國蹶鼠(Sicista con‐color)、灰蹶鼠(S. pseudonapaea) 和天山蹶鼠(S.tianshanica) 的分布預測結果的主要影響因素；年降水是塔里木跳鼠和三趾心顱跳鼠(Salpingotus ko‐zlovi) 的主要影響因素；降水季節(jié)性變化和8 月太陽輻射分別是長尾蹶鼠(S. caudata)、草原蹶鼠(S. subtilis)和大部分跳鼠科動物的主要影響因素；溫度季節(jié)性變化僅為準噶爾羽尾跳鼠(Stylodipus sungorus)預測結果的主要影響因素。

圖3 中國跳鼠總科各科物種適宜度分布區(qū). A：蹶鼠科；B：林跳鼠科；C：跳鼠科Fig. 3 Suitable distribution areas of Dipodoidea in China. A:Sicisti‐dae;B:Zapodidae;C:Dipodidae

表3 環(huán)境因子對跳鼠總科物種分布模型預測的貢獻百分比Table 3 Percentage contributions of environmental variables

3 討論

3.1 系統(tǒng)發(fā)育關系

本文基于前人的研究基礎(Lebedevet al.,2013; Pisanoet al., 2015; 程繼龍等，2021)，新增在我國發(fā)現并命名的奇美跳鼠和由原三趾跳鼠努日亞種提升為種的塔里木跳鼠，同時增加新分布記錄物種灰蹶鼠和準噶爾羽尾跳鼠，利用12 個基因的串聯(lián)矩陣重建了我國最新的跳鼠總科物種的系統(tǒng)發(fā)育關系。該系統(tǒng)發(fā)育關系支持我國跳鼠總科動物為3科12屬22種的結論(程繼龍等，2021)，從分子角度為我國跳鼠總科物種分類關系提供參考。

蹶鼠科中，天山蹶鼠和中國蹶鼠位于系統(tǒng)發(fā)育樹的最基部，與Lebedev等(2019)關于蹶鼠科研究結果一致；灰蹶鼠雖然生活在泰加林緣、灌木草原等相對濕潤地帶，但與其親緣關系最近的是地理分布上最近的草原蹶鼠，且二者位于蹶鼠科系統(tǒng)發(fā)育的最頂端。五趾跳鼠亞科中，本研究得出的系統(tǒng)發(fā)育結果支持Lebedev 等(2013) 提出的傳統(tǒng)的Allac‐taga屬(包 括Allactaga、Orientallactaga、Scratu‐rus)是一個并系群的觀點；同時表明東方五趾跳鼠屬(由巨泡五趾跳鼠、巴里坤跳鼠和五趾跳鼠O.sibirica組成)提升為獨立的屬級分類階元是合適的(Michaux and Shenbrot, 2017; Chenget al., 2021)。雖然本研究系統(tǒng)發(fā)育樹中巨泡五趾跳鼠和巴里坤跳鼠分支結點的拓撲結構支持率不高，但巨泡五趾跳鼠和巴里坤跳鼠在頭骨結構和身體形態(tài)等方面有明顯的差異(Sokolovet al., 1981; Miljutin, 2008;Al‐hajeri,2021)，二者的種級階元地位可靠。拓撲結點的支持率不高可能與東方五趾跳鼠屬內物種形成時間較短，及其可能存在種間雜交有關，有待未來更深入研究(Chenget al.,2021)。

3.2 中國跳鼠總科的分布格局及影響因素

我國分布的22 種跳鼠總科動物中，生活在荒漠或半荒漠生境的跳鼠科動物有16 種，占跳鼠總科物種的72.73% (程繼龍等，2021)?？赡苡捎谔罂偪苿游锏暮笾Y構和跳躍能力起源于森林和草原環(huán)境(Wuet al., 2014)，卻在荒漠或半荒漠這種相對開闊的生境中得到進一步發(fā)展，演化出更多的物種。因此，造成我國跳鼠物種豐富度較高的地區(qū)并不是那些傳統(tǒng)認知中生物多樣性較高的區(qū)域，而是集中在我國西北部，這與我們物種豐富度的研究結果一致。物種豐富度較高的區(qū)域集中在我國新疆地區(qū)，可能與新疆獨特的地理位置有關。山地與盆地交互的地理結構形式為山地森林、草原和荒漠多種生境的形成奠定了地理基礎(陳曦等，2015)，也為跳鼠科和蹶鼠科動物的生存提供了生境保障；同時新疆位于歐洲—西伯利亞亞界阿爾泰—薩彥嶺區(qū)、中亞亞界哈薩克斯坦區(qū)和中亞亞界蒙新區(qū)3個動物地理區(qū)系交匯處，物種豐富度也相對較高。而我國東北和四川地區(qū)森林生境較多，其他森林嚙齒動物占據了較多的生態(tài)位，不利于跳鼠總科動物的分化擴張，因此物種豐富度較低。

不同生態(tài)型的跳鼠動物分布的主要影響因素不同。蹶鼠科和林跳鼠科的分布主要受年均溫、降水季節(jié)性變化和8 月太陽輻射的影響，其中中國蹶鼠、灰蹶鼠和天山蹶鼠都分布于中高海拔的森林或林緣草地生境，降水和日照都相對充足，推測其受年均溫的影響較大；而草原蹶鼠與其他4種蹶鼠不同，分布于生境相對干旱的草原生境(Lebedevet al.,2020)，推測其受降水的影響較大；長尾蹶鼠分布預測結果主要受8月太陽輻射的影響，可能與其分布于我國東北高緯度地區(qū)有關。四川林跳鼠是我國特有種，其分布區(qū)域與中國蹶鼠有部分重疊，也分布于降水相對充沛的森林生境，受年均溫的影響較大。跳鼠科動物分布于荒漠或半荒漠生境(Lebe‐devet al., 2013; Zhanget al., 2013; Pisanoet al.,2015; Michaux and Shenbrot, 2017)，8 月 太 陽 輻 射和降水，包括年降水和降水季節(jié)性變化，成為影響其分布的主要因素。8 月太陽輻射的影響可能由于跳鼠總科動物多數在9、10月進入冬眠，第二年的4、5 月醒來。太陽輻射對植物營養(yǎng)物質的累積(Monteith, 1977) 和植物分布(Waltonet al., 2005)產生影響，8 月太陽輻射可能通過影響植被生長，間接影響跳鼠總科動物越冬所需營養(yǎng)物質的累積，成為影響分布的重要環(huán)境因素。而準噶爾羽尾跳鼠與同屬的羽尾跳鼠相比，分布范圍更窄，活動時間和繁殖時間也相對縮短(Michaux and Shenbrot,2017)，更集中在溫度相對較高的時間段內，可能溫度的季節(jié)性變化對其分布的影響更大。

3.3 生物多樣性的研究與應用

物種是生物多樣性研究的基礎單位，明確物種的系統(tǒng)發(fā)育關系和分布格局是認識物種的基礎和關鍵。分子系統(tǒng)發(fā)育研究的不斷發(fā)展，極大地豐富了跳鼠總科的多樣性(Quereet al.,2009;Pisa‐noet al., 2015; Shenbrotet al., 2017; Chenget al.,2018; Bannikovaet al., 2019; Lebedevet al., 2019;Chenget al.,2021)。我國跳鼠總科的系統(tǒng)發(fā)育和分布格局的研究也表明，荒漠地區(qū)并不如想象中的那么“貧瘠”，其物種豐富度相對較高(李娜等，2020)，同樣應該成為生物多樣性研究的關注地。同時研究結果也揭示，我們對跳鼠總科多樣性的認識仍不夠完善，稀有物種的本底調查工作不夠深入，如蹶鼠科中灰蹶鼠、中國蹶鼠的樣本量嚴重缺失；心顱跳鼠亞科中五趾心顱跳鼠(Cardio‐cranius paradoxus) 的樣本量和針對其所進行的研究工作不足；奇美跳鼠僅在寧夏西吉縣南部有狹窄分布，種群數量較低(Shenbrotet al.,2017)，缺乏對其進一步的研究工作。應用物種分布模型對物種豐富度和適宜度進行研究，可以為本底調查提供參考，有效促進生物多樣性的研究工作。

對于不同生境類型的物種，我們關注的環(huán)境因素側重點應有所不同。對于森林生境，我們更應該關注溫度變化對物種的影響；而對于生活在草原或荒漠生境的物種，則更應該關注降水和太陽輻射的影響。嚙齒動物的種群波動是環(huán)境變化的重要指標(Utreraet al.,2000)，密切關注跳鼠總科不同譜系、生境類型和分布區(qū)域的種群變化，有效控制種群波動，對保護生物多樣性和維持生態(tài)系統(tǒng)的平衡與穩(wěn)定具有重要意義。基于得到的物種豐富度和適宜性分布格局，可以對野外工作進行指導，填補并擴充當下的物種數據，為物種本底調查提供參考。

?