中国跳鼠总科物种系统发育与分布格局

隋璐璐 程继龙 夏霖 杨奇森*

（1 中国科学院动物研究所，北京 100101）（2 中国科学院大学，北京 100049）

随着物种灭绝速率的不断加快，生物多样性研究和保护工作正成为生物学研究的重中之重(Daviset al., 2018)，正确认识物种并明确物种的分布状况是生物多样性研究和保护的重要前提。但是，仅依赖传统形态分类学方法进行的物种分类已不利于对隐性物种多样性的发掘(Sheth and Thaker,2017)，可能导致生物多样性被低估。随着分子生物学技术的发展和系统分类学、演化生物学研究的不断深入，DNA 序列变异这一客观特征对物种鉴定和物种间系统发育关系的构建起到了极大的促进作用(Yang,2006)。同时，厘清物种间的系统发育关系，也有助于更深入认识生物多样性的内涵，为物种的保护工作提供理论基础(Pérez-Losadaaet al.,2002;Woodet al.,2019)。

物种分布模型(Species distribution model,SDM)能利用物种分布信息和当地气候变量之间的关系预测物种潜在分布区(Guisan and Zimmer‐mann, 2000; Ma and Sun, 2018; Mwakapejeet al.,2019; Zhanget al., 2019)。中国的国土面积辽阔，生物学调查难以全部覆盖，因此，SDM 可以在一定程度上指导未来的野外调查方向，为进一步探索和发现物种分布提供参考，为物种保护措施的制定提供科学基础(李娜等，2020)。在多种分布模型的模拟方法中，最大熵模型(Maximum entropy model,Maxent)可以在样本量较小的情况下，利用环境变量和物种分布位点来计算约束条件，并探索在此约束条件下最大熵的可能分布，进而预测研究区物种的生境适宜性(Merowet al.,2013)，得出优于其他模型的模拟结果(Phillipset al.,2006)。

跳鼠，因其高度特化的后肢结构和二足跳跃的运动方式而得名，广泛分布于我国中西部的诸多地区，跨越森林、草原和荒漠3种不同的典型生境，是我国生物多样性尤其是荒漠地区生物多样性的重要指示类群。同时，我国还是跳鼠物种的起源、扩散和分布中心(Zhanget al., 2013; Pisanoet al., 2015)。程继龙等(2021) 对中国跳鼠总科物种的分类系统进行了综述，认为我国共分布有跳鼠总科动物3 科6 亚科12 属22 种，是世界上跳鼠总科物种丰富度最高的国家之一。以前的研究中，无论是从分子角度(Zhanget al., 2013; Pisanoet al., 2015) 或是形态结构角度(Wuet al., 2014;Mooreet al., 2015)，其系统发育关系中都未能包括我国分布的全部跳鼠，尤其是最近几年新发现和修订的类群，如奇美跳鼠(Chimaerodipus auri‐tus) (Shenbrotet al., 2017) 和塔里木跳鼠(Dipus deasyi) (Chenget al., 2018)。因此，我国跳鼠总科新的系统发育关系亟待完善。同时，各个物种的分布状况也缺乏全面调查和整理，难以满足对该类群本底资料的掌握和认识生物多样性的需要。

本研究以中国分布的跳鼠总科动物为研究对象，通过研究物种DNA 序列上的变异，厘清物种的系统发育关系，探讨我国跳鼠总科内部在科、属和种3个分类阶元上的演化关系，为我国跳鼠总科分类提供理论基础。同时，通过物种分布模型，推测跳鼠总科物种在我国的潜在分布范围，并绘制物种丰富度格局，分析不同物种分布状况的主要影响因素，为我国跳鼠总科多样性的研究工作提供参考。

1 研究方法

1.1 系统发育关系重建

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)进行筛选，排除不确定性序列和测序质量差的序列，并尽可能降低亚种干扰，选择物种分布中心区域的序列，最终挑选出每个基因中最具代表性的基因序列。

表1 本研究应用的物种Table 1 The species used in this study

1.1.2 系统发育分析

利用MEGA-X 中的Muscle (Edgar, 2004) 进行序列的比对，将序列两端起始的部分碱基删除后导出保存。再利用Sequence Matrix 1.8(Vaidyaet al.,2011)串联基因，构建基因矩阵。基于所获得的基因矩阵，分别采用最大似然法(maximum likelihood,ML)和贝叶斯法(Bayesian inferences,BI)进行系统发育分析。运用Phylosuite(Zhanget al.,2020)软件中的PartitionFinder2 (Lanfearet al., 2016)，基于“贪婪算法”(greedy search algorithm)和贝叶斯信息准则(Bayesian information criterion,BIC) 选择最优分区方式和最佳核酸替换模型。使用IQ-TREE(Guindonet al., 2010;Anisimovaet al., 2011; Minhet al.,2013;Nguyenet al.,2015)构建ML 树，执行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 环境因子选择

依据分布数据特点和跳鼠总科动物生活习性，分别从全球气候数据库(WorldClim2, www. world‐clim. org) 和地理空间数据云数据库(www. gs‐cloud.cn)下载温度、降水、气象和地形4 种变量，并转换为分辨率为30 s的环境数据变量；考虑到跳鼠总科动物冬眠的习性，对于太阳辐射(sard)、水蒸气压(vapr)和风速(wind)3 个气象数据变量只选取其活动季节(4—9 月)数据。最终获得23 个环境变量，共38个栅格图层(表2)。

表2 跳鼠总科物种分布模型使用的环境变量Table 2 Environment variables used for species distribution model of Dipodoidea

1.2.3 Maxent模型选择

部分研究表明，Maxent 的默认参数设置是基于最大经验调优的结果，可能导致模型性能较差(Radosavljevicet al., 2013; Ma and Sun, 2018)。因此我们结合物种分布位点和筛选出来的环境变量，利用R 语言(R 4.0.5)的ENMeval 包(Muscarellaet al.,2014)对Maxent 进行模型优化。设置最大背景点数为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)值最小的参数设置，作为最优参数进行物种分布模型计算。

应用Maxent 3.4.1 软件，设置随机检验百分比为25%，模型重复运行30 次，重复运行类型为二次抽样(subsample)，并使用刀切法(jackknife test)分析环境变量的重要性。以ROC 曲线(receiv‐er operating characteristic curve) 对模型的准确度进行评估。ROC 曲线与横轴围成的面积为AUC 值(Hanley and McNeil, 1982)，可用来衡量模型的优劣。一般而言，AUC 值为0~0.6 表示模型预测无效；0.6 ~ 0.7 表示模型准确度较差；0.7 ~ 0.8 表示模型准确度一般；0.8 ~ 0.9 表示模型准确度良好；大于0.9表示模型准确度极好。

1.2.4 物种适宜度和丰富度

由于每个物种对环境的耐受程度不一样，基于可查到的分布记录的临界值，对每个物种的适宜栖息地阈值进行划分(Svenninget al.,2008)。然后依据阈值提取每个物种适宜栖息地图层，叠加后得到物种适宜分布图。对物种分布预测结果的图层进行二值化处理，即每个栅格中该物种适宜度大于阈值则设置为1，否则为0。然后将物种二值分布图进行空间叠加，获得中国跳鼠总科物种丰富度分布格局。

2 结果

2.1 系统发育关系

构建系统发育关系的12 个基因串联矩阵总长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)。跳鼠总科内的各亚科之间系统发育关系明确，跳鼠亚科、五趾跳鼠亚科、长耳跳鼠亚科(Euchoreutinae)、心颅跳鼠亚科聚为跳鼠科，跳鼠科与林跳鼠科(Zapodidae)的亲缘关系更近，蹶鼠科(Sicistidae)位于跳鼠总科最基部。

图1 基于12个基因的串联矩阵构建的跳鼠总科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 物种丰富度格局

基于分布位点得到的物种丰富度结果表明，跳鼠总科主要分布在我国北方的荒漠、半荒漠地区，青藏高原东北缘的山地与高山草原地区。物种丰富度最高的分布区位于新疆北部包括昌吉州、阿勒泰地区、塔城地区、博尔塔拉州和伊犁等地，一般分布有7 种以上跳鼠总科动物。其中，塔城、额敏、裕民和托里部分地区甚至包含9~10种跳鼠总科动物(图2)。物种丰富度次高一级的分布区集中在新疆西南部、内蒙古中西部和甘肃省河西走廊地区，分布有4种以上跳鼠总科动物。其他区域的物种丰富度较低，一般只有2种以下。

图2 中国跳鼠总科物种丰富度分布格局Fig. 2 Species richness of Dipodoidea in China

2.3 物种适宜度及主要影响因素

跳鼠总科在我国适宜分布区域总面积约2.22×106km2，约占国土总面积的23%。其中，蹶鼠科适宜分布面积14.15×104km2(图3A)；林跳鼠科适宜分布面积11.85×104km2(图3B)；跳鼠科适宜分布面积达20.423×105km2(图3C)。21 个物种分布模型ROC 曲线的AUC 值均高于0.889，说明Maxent模拟跳鼠总科物种的潜在分布范围的模型拟合度相对较高，结果可信。环境变量对物种分布模型预测的贡献百分比结果表明，影响不同物种分布的主要环境因子不尽相同(表3)。年均温是四川林跳鼠(Eozapus setchuanus)、中国蹶鼠(Sicista con‐color)、灰蹶鼠(S. pseudonapaea) 和天山蹶鼠(S.tianshanica) 的分布预测结果的主要影响因素；年降水是塔里木跳鼠和三趾心颅跳鼠(Salpingotus ko‐zlovi) 的主要影响因素；降水季节性变化和8 月太阳辐射分别是长尾蹶鼠(S. caudata)、草原蹶鼠(S. subtilis)和大部分跳鼠科动物的主要影响因素；温度季节性变化仅为准噶尔羽尾跳鼠(Stylodipus sungorus)预测结果的主要影响因素。

图3 中国跳鼠总科各科物种适宜度分布区. A：蹶鼠科；B：林跳鼠科；C：跳鼠科Fig. 3 Suitable distribution areas of Dipodoidea in China. A:Sicisti‐dae;B:Zapodidae;C:Dipodidae

表3 环境因子对跳鼠总科物种分布模型预测的贡献百分比Table 3 Percentage contributions of environmental variables

3 讨论

3.1 系统发育关系

本文基于前人的研究基础(Lebedevet al.,2013; Pisanoet al., 2015; 程继龙等，2021)，新增在我国发现并命名的奇美跳鼠和由原三趾跳鼠努日亚种提升为种的塔里木跳鼠，同时增加新分布记录物种灰蹶鼠和准噶尔羽尾跳鼠，利用12 个基因的串联矩阵重建了我国最新的跳鼠总科物种的系统发育关系。该系统发育关系支持我国跳鼠总科动物为3科12属22种的结论(程继龙等，2021)，从分子角度为我国跳鼠总科物种分类关系提供参考。

蹶鼠科中，天山蹶鼠和中国蹶鼠位于系统发育树的最基部，与Lebedev等(2019)关于蹶鼠科研究结果一致；灰蹶鼠虽然生活在泰加林缘、灌木草原等相对湿润地带，但与其亲缘关系最近的是地理分布上最近的草原蹶鼠，且二者位于蹶鼠科系统发育的最顶端。五趾跳鼠亚科中，本研究得出的系统发育结果支持Lebedev 等(2013) 提出的传统的Allac‐taga属(包 括Allactaga、Orientallactaga、Scratu‐rus)是一个并系群的观点；同时表明东方五趾跳鼠属(由巨泡五趾跳鼠、巴里坤跳鼠和五趾跳鼠O.sibirica组成)提升为独立的属级分类阶元是合适的(Michaux and Shenbrot, 2017; Chenget al., 2021)。虽然本研究系统发育树中巨泡五趾跳鼠和巴里坤跳鼠分支结点的拓扑结构支持率不高，但巨泡五趾跳鼠和巴里坤跳鼠在头骨结构和身体形态等方面有明显的差异(Sokolovet al., 1981; Miljutin, 2008;Al‐hajeri,2021)，二者的种级阶元地位可靠。拓扑结点的支持率不高可能与东方五趾跳鼠属内物种形成时间较短，及其可能存在种间杂交有关，有待未来更深入研究(Chenget al.,2021)。

3.2 中国跳鼠总科的分布格局及影响因素

我国分布的22 种跳鼠总科动物中，生活在荒漠或半荒漠生境的跳鼠科动物有16 种，占跳鼠总科物种的72.73% (程继龙等，2021)。可能由于跳鼠总科动物的后肢结构和跳跃能力起源于森林和草原环境(Wuet al., 2014)，却在荒漠或半荒漠这种相对开阔的生境中得到进一步发展，演化出更多的物种。因此，造成我国跳鼠物种丰富度较高的地区并不是那些传统认知中生物多样性较高的区域，而是集中在我国西北部，这与我们物种丰富度的研究结果一致。物种丰富度较高的区域集中在我国新疆地区，可能与新疆独特的地理位置有关。山地与盆地交互的地理结构形式为山地森林、草原和荒漠多种生境的形成奠定了地理基础(陈曦等，2015)，也为跳鼠科和蹶鼠科动物的生存提供了生境保障；同时新疆位于欧洲—西伯利亚亚界阿尔泰—萨彦岭区、中亚亚界哈萨克斯坦区和中亚亚界蒙新区3个动物地理区系交汇处，物种丰富度也相对较高。而我国东北和四川地区森林生境较多，其他森林啮齿动物占据了较多的生态位，不利于跳鼠总科动物的分化扩张，因此物种丰富度较低。

不同生态型的跳鼠动物分布的主要影响因素不同。蹶鼠科和林跳鼠科的分布主要受年均温、降水季节性变化和8 月太阳辐射的影响，其中中国蹶鼠、灰蹶鼠和天山蹶鼠都分布于中高海拔的森林或林缘草地生境，降水和日照都相对充足，推测其受年均温的影响较大；而草原蹶鼠与其他4种蹶鼠不同，分布于生境相对干旱的草原生境(Lebedevet al.,2020)，推测其受降水的影响较大；长尾蹶鼠分布预测结果主要受8月太阳辐射的影响，可能与其分布于我国东北高纬度地区有关。四川林跳鼠是我国特有种，其分布区域与中国蹶鼠有部分重叠，也分布于降水相对充沛的森林生境，受年均温的影响较大。跳鼠科动物分布于荒漠或半荒漠生境(Lebe‐devet al., 2013; Zhanget al., 2013; Pisanoet al.,2015; Michaux and Shenbrot, 2017)，8 月 太 阳 辐 射和降水，包括年降水和降水季节性变化，成为影响其分布的主要因素。8 月太阳辐射的影响可能由于跳鼠总科动物多数在9、10月进入冬眠，第二年的4、5 月醒来。太阳辐射对植物营养物质的累积(Monteith, 1977) 和植物分布(Waltonet al., 2005)产生影响，8 月太阳辐射可能通过影响植被生长，间接影响跳鼠总科动物越冬所需营养物质的累积，成为影响分布的重要环境因素。而准噶尔羽尾跳鼠与同属的羽尾跳鼠相比，分布范围更窄，活动时间和繁殖时间也相对缩短(Michaux and Shenbrot,2017)，更集中在温度相对较高的时间段内，可能温度的季节性变化对其分布的影响更大。

3.3 生物多样性的研究与应用

物种是生物多样性研究的基础单位，明确物种的系统发育关系和分布格局是认识物种的基础和关键。分子系统发育研究的不断发展，极大地丰富了跳鼠总科的多样性(Quereet al.,2009;Pisa‐noet al., 2015; Shenbrotet al., 2017; Chenget al.,2018; Bannikovaet al., 2019; Lebedevet al., 2019;Chenget al.,2021)。我国跳鼠总科的系统发育和分布格局的研究也表明，荒漠地区并不如想象中的那么“贫瘠”，其物种丰富度相对较高(李娜等，2020)，同样应该成为生物多样性研究的关注地。同时研究结果也揭示，我们对跳鼠总科多样性的认识仍不够完善，稀有物种的本底调查工作不够深入，如蹶鼠科中灰蹶鼠、中国蹶鼠的样本量严重缺失；心颅跳鼠亚科中五趾心颅跳鼠(Cardio‐cranius paradoxus) 的样本量和针对其所进行的研究工作不足；奇美跳鼠仅在宁夏西吉县南部有狭窄分布，种群数量较低(Shenbrotet al.,2017)，缺乏对其进一步的研究工作。应用物种分布模型对物种丰富度和适宜度进行研究，可以为本底调查提供参考，有效促进生物多样性的研究工作。

对于不同生境类型的物种，我们关注的环境因素侧重点应有所不同。对于森林生境，我们更应该关注温度变化对物种的影响；而对于生活在草原或荒漠生境的物种，则更应该关注降水和太阳辐射的影响。啮齿动物的种群波动是环境变化的重要指标(Utreraet al.,2000)，密切关注跳鼠总科不同谱系、生境类型和分布区域的种群变化，有效控制种群波动，对保护生物多样性和维持生态系统的平衡与稳定具有重要意义。基于得到的物种丰富度和适宜性分布格局，可以对野外工作进行指导，填补并扩充当下的物种数据，为物种本底调查提供参考。

?