濒危的玻利维亚鼠兔（Abrocoma boliviensis）的人口遗传学与谱系结构

《Conservation Genetics》：Population genetics and lineage structure of the endangered Bolivian chinchilla rat Abrocoma boliviensis

　　摘要
对濒危物种的保护遗传学研究可以帮助确定哪些种群应该成为管理计划的重点。濒危的玻利维亚毛丝鼠（Abrocoma boliviensis）目前正受到其稀有性以及人类活动导致的栖息地变化的威胁。鉴于其保护状况、有限的地理分布范围和关于其自然历史的有限信息，了解其遗传多样性的分布对于制定任何潜在的保护措施至关重要。在这项研究中，我们评估了A. boliviensis的遗传多样性和种群结构，作为对该物种进行全面评估的第一步。来自12个个体的线粒体数据显示出高水平的遗传距离、核苷酸多样性和多态性，这些都表明存在三个独立的谱系。这一发现得到了降低表示基因组数据的支持，这些数据表明这些谱系之间几乎没有基因流。这种基因流的缺乏表明这些谱系经历了不同的进化路径，至少应该被视为具有进化意义的独立单元。我们的研究强调了开展调查工作的紧迫性，这是首要行动。生成新的种群级数据对于深入理解该物种、明确其谱系的进化路径以及制定有效的保护策略至关重要。

引言
分子方法和群体遗传分析已成为保护研究中的常规做法。保护遗传学领域识别与低密度和小种群规模相关的威胁，并提供工具以更好地理解由人口统计因素驱动的过程（Willi等人，2022年）。它还可以通过确定哪些种群应该是保护工作的重点来为物种的管理和保护做出贡献。这可以通过识别进化意义上的单元（ESUs）来实现（Ryder，1986年；Casacci等人，2014年）。当种群在遗传和生态上高度分化时，就会认定这些单元，因为每个单元可能处于不同的进化路径上，并具有潜在的适应性（Funk等人，2012年；Willi等人，2022年）。例如，对濒危或其他受威胁物种的群体遗传学研究发现了哺乳动物（Degner等人，2007年；Cossíos等人，2012年；Teixeira等人，2023年）、鸟类（Quintela等人，2010年；Murphy等人，2011年）、爬行动物和两栖动物（Borges等人，2018年；Muniz等人，2018年；Castillo-Morales等人，2023年）中存在独立ESUs的证据。了解物种的遗传状况至关重要，因为当发生遗传枯竭时，小种群特别容易受到威胁，这可能危及它们在环境变化面前的生存（O’Brien，1994年）。这种现象对于栖息在自然破碎化生境中的物种尤其相关，在这些生境中，遗传漂变的影响可能会被放大（Willi等人，2022年）。

自20世纪50年代以来，玻利维亚的景观经历了多种因素（如原住民群体、殖民者和企业）和实践活动（如传统农业、机械化农业、养牛或森林利用）的影响，大部分定居点和森林砍伐发生在圣克鲁斯市和拉巴斯市附近的Yungas地区（Killeen等人，2008年）。玻利维亚东部山脉的土地覆盖已转变为农田或裸地，取代了过去的森林和牧场。这一趋势在海拔2000米以下的肥沃山谷和中等坡度地区尤为明显，这些地区已被广泛用于农业（Brandt和Townsend，2006年）。该地区的本地物种——玻利维亚毛丝鼠（Abrocoma boliviensis）仅分布于玻利维亚东部安第斯山脉的中高海拔森林中，被认为是一种稀有物种，因为自1990年首次描述以来，关于该物种的记录不到20个个体（Tarifa等人，2009年；Quinteros-Mu?oz，2015年；Hidalgo-Cossio等人，2016年；Quiroga Pacheco等人，2020年）。

一个物种的稀有性可以根据其数量或地理分布范围来确定，虽然并非所有稀有物种都必然面临灭绝的风险，但它们的低数量和/或受限的分布使它们更容易受到随机事件的影响（Drever等人，2012年）。在国际自然保护联盟（IUCN）濒危物种红色名录的最新评估中，A. boliviensis被列为极度濒危物种（标准B1ab(i, ii, ii)）（Bernal，2016年）。这意味着该物种的分布范围小于100平方公里，严重破碎化，或者仅在一个地点存在（A. boliviensis的分布仅限于圣克鲁斯省的模式产地），并且其栖息地的面积、范围和/或质量持续下降。此外，该物种被列入全球100种EDGE（进化独特且全球濒危）哺乳动物物种之中，这一分类适用于那些在进化上高度独特且面临灭绝威胁的物种（Gumbs等人，2023年）。

近年来，来自各种研究项目的采样工作在波托西省、科恰班巴省、塔里哈省和圣克鲁斯省新增了8个地点（Tarifa等人，2009年；Quinteros-Mu?oz，2015年；Hidalgo-Cossio等人，2016年；Quiroga Pacheco等人，2020年），显著扩展了该物种在玻利维亚东部山脉的分布范围，涵盖了两个主要的植物地理区域：Yungas和Boliviano-Tucumano森林（Josse等人，2011年）。Yungas森林从秘鲁延伸到玻利维亚中部，分布在科恰班巴省和圣克鲁斯省的部分地区，覆盖了广泛的海拔梯度（500-4000米），之后东部山脉的南部主要由Boliviano-Tucumano森林组成，其间散布着圣克鲁斯省和塔里哈省的干燥安第斯森林。这种过渡大约发生在南纬18度处。随着这些新地点的发现，A. boliviensis现在分布在由深谷和陡坡分隔的混合森林景观中，包括山地森林、干燥森林和亚热带森林（即Yungas、干燥安第斯森林和Tucumano森林）（Rex和Hanratty，1989年；Graham等人，2001年）（图1）。

材料与方法
我们使用了来自十个地点的十二个A. boliviensis标本，包括一个模式标本（图1；表S1）。所有标本均通过机构间的借阅从自然历史收藏中获得。我们分析中的所有组织样本都附有凭证标本（表S1）。

DNA提取、PCR扩增和细胞色素b基因测序
首先使用线粒体标记物（细胞色素b或Cyt-b基因）评估了A. boliviensis个体之间的关系。根据Qiagen公司的DNEasy? Blood and Tissue Kit制造说明从肌肉或肝脏组织中提取基因组DNA，并按照Salazar-Bravo等人（2023年）的协议进行PCR扩增。其中一个A. boliviensis样本是皮肤片段，首先按照de Moraes-Barros和Morgante（2007年）的协议进行清洗和补水处理。该样本被分成两部分，一部分使用Qiagen试剂盒处理，另一部分使用Zymo Research公司的Quick-DNA Fecal/Soil试剂盒处理。然后将皮肤提取物合并用于后续的PCR程序。该提取工作在一个之前未处理过哺乳动物的实验室中进行。PCR扩增结果在1%琼脂糖凝胶中观察，成功的扩增样本被送往Psomagen USA（马里兰州）进行标准Sanger测序。单个序列在DNASTAR v 5.52中组装成连续序列，并使用Mesquite v3.61软件（Maddison和Madison，2019年）的MUSCLE选项进行比对和编辑（Edgar，2004年）。最终矩阵包括十二个A. boliviensis个体的810 bp细胞色素b序列。有关所用引物和PCR循环条件的详细信息作为补充材料提供。

细胞色素b基因分析
使用IQ-TREE网络服务器（Trifinopoulos等人，2016年）通过最大似然（ML）搜索构建了细胞色素b基因树。指定了按密码子划分的方法，但允许IQ-TREE确定数据的最佳拟合替代模型。通过1000次迭代的快速自举法评估了谱系支持度。树以中间点为根。利用Abrocomidae家族其他物种的数据，通过Mega 10中实现的Kimura 2参数模型计算了种间和种内遗传变异（包括基于ML树结果分配给A. boliviensis的谱系间的距离）（见结果部分和图2）。

A. boliviensis的种群结构
为了进一步探索A. boliviensis的遗传多样性，使用DnaSP v6（Rozas等人，2017年）计算了核苷酸多样性（π）、单倍型数量（h）和多样性（Hd）以及多态位点。此外，在DnaSP中还计算了Fu的Fs检验和Tajima的D检验以推断物种的历史人口统计特征。

降低表示基因组水平的变异
为了确定从细胞色素b获得的遗传距离和核苷酸多样性水平是否确实代表了该物种的进化历史，而不是受到线粒体基因组独特特征（即母系遗传）的偏见，我们使用全基因组降低表示方法探讨了A. boliviensis的种群结构。本研究处理的一部分样本（n=9）按照Elshire等人（2011年）在威斯康星大学生物技术中心描述的基因分型测序（GBS）方法进行了测序。样本用nsiI和bfaI酶消化，然后使用Illumina NovaSeq 6000仪进行测序，产生成对的末端读段（2×150 bp）。GBS方法产生的非条形R2s可能会在Stacks 2中产生偏差（Rochette等人，2019年）。因此，所有下游分析都是使用德克萨斯理工大学的高性能计算中心的单端数据集进行的。

**位点筛选和SNP调用**
在解复用、大小筛选（130 bp）以及对读取数据进行质量控制检查之后，我们在STACKS 2（Rochette等人，2019年）中运行了多次denovo_map程序的迭代，修改了参数-m（决定了ustacks模块中的覆盖深度）、-M（决定了在ustacks模块中允许形成的错配数量）和-n（控制在构建cstacks目录时个体之间允许的错配数量），以找到最适合我们数据的组合用于下游分析。该过程的详细描述在补充材料中提供。大多数参数组合提供了质量相似的结果（即样本的分布和分组，以及在PCA中解释的方差量），当至少80%的个体共享读取数据时（-r = 0.8）。在VCFtools中，这些r80文件进一步使用以下选项进行了筛选：--max-meanDP 50（以移除深度值较高的位点）、--minDP 7（以确保所有基因型至少有7个读取数据）和--thin 150（每个位点只选择一个SNP）。虽然通常基于最小等位基因频率进行等位基因筛选（部分是为了去除测序错误），但我们决定不这样做，因为这个数据集的个体数量很少，在某些情况下只有一个个体代表一个遗传群，因此任何基于MAF的筛选都可能移除生物学上相关的等位基因。最后，对于下游分析，使用了Stacks设置M = 4和n = 5生成的文件，因为它们在大量变异体和避免组合同源或重复位点之间提供了更好的平衡（Paris等人，2017年）。

**群体分析**
使用上述过滤后的数据集，我们使用PLINK v1.9进行了PCA分析。然后我们使用相同的数据集再次运行STACKS的群体模块并生成了一个结构文件。群体结构使用STRUCTURE v2.3.4软件进行了评估（Pritchard等人，2000年）。根据开发者的建议，分析假设存在混合（NOADMIX = 0），提供了采样位置以帮助聚类（LOCPRIOR = 1），并指定了独立或相关的等位基因频率模型（FREQSCORR = 0和1）。提供采样位置作为先验，因为我们的数据集个体数量少且SNP数量相对较少，这限制了软件检测真实群体结构的能力。此外，Hubisz等人（2009年）开发的模型评估了采样位置是否具有信息性；只有当它们具有信息性时，软件才会使用它们。使用程序的命令行版本，分析在K值从1到9的范围内进行了运行。每个K值运行了20次，BURNIN = 200,000和NUMREPS = 500,000。为了确定最合适的K值，我们遵循了Evanno方法（Evanno等人，2005年），并通过STRUCTURE HARVESTER软件实现（Earl和vonHoldt，2012年）。我们使用Structure Plot v2.0程序生成了条形图（Ramasamy等人，2014年）。最后，使用adegenet包和Nei的距离度量标准（Jombart和Ahmed，2011年）在R v4.0.5中测试并绘制了基于距离的隔离（IBD）图。

**Cyt-b数据集分析**
对12个A. boliviensis个体的线粒体数据进行分析，发现存在显著的遗传距离、高核苷酸多样性和高水平的多态性，共同表明存在三个不同的支系。A. boliviensis中Cyt-b序列之间的遗传距离相当大，范围从0.3%到13.9%（表1）。但是A. boliviensis内部的最高距离（即13.9%）小于A. boliviensis与其姐妹物种A. cinerea之间的距离（平均17.4%）。这些分化水平支持将A. boliviensis的样本分类为三个独立的支系。支系A包括来自Cochabamba的所有样本、一个来自Potosi的样本以及一个来自Santa Cruz的样本（这代表了该物种在La Paz和Cochabamba中部温带山谷的分布北部），支系B包括来自Santa Cruz省最西部山区雨林边缘的样本TK161712 + TK161714 + MNK3851，支系C由Tarija省已知分布南部边缘的样本MHNC-M682代表（图2a：分别为黄色、紫色和黑色星号）。

**表1** 使用Kimura 2参数从cytochrome b序列计算出的Abrocoma boliviensis各支系之间的平均遗传距离（%）。包括与其姐妹物种Abrocoma cinerea的遗传距离值以供比较。对角线上的粗值代表支系内的距离。

**表2** 核苷酸多样性和多态性位点的值高度依赖于A. boliviensis样本是否被分组（表2）。当所有样本（不论属于哪个支系）被合并时，核苷酸多样性为0.045，存在121个多态性位点。一旦将boliviensis样本分组到支系中，就清楚地看到支系C的样本夸大了整个物种的数值。当移除这个支系（即合并支系A和B）时，核苷酸多样性为0.033，而多态性位点为68个。鉴于所有A. boliviensis样本都来自不同的地点，每个个体（Pusuq’huni的两个样本除外）在支系A和B中都有独特的单倍型，这反映在它们的高单倍型多样性上。Fs值对于A. boliviensis是正的，但只有当合并支系A和B时才显著。Tajima的D值也显示出类似的趋势（表2）。

**简化表示的基因组水平变异**
解复用后，保留的读取数据量为253,725,812个。运行denovo_map后，位点和多态性位点的数量分别为1,322,883个和1,002,230个。最后，在VCFtools中应用所有过滤器后保留的变异体数量为2,192个，每个读取数据的平均深度为17倍。

**PCA图**
基于GBS数据构建的PCA图支持了Cyt-b树中形成的组（图2b），其中PC1解释了支系C（Tarija）样本与其他样本之间51.9%的方差，PC2解释了支系B（Santa Cruz）样本中20.4%的方差。由于它们在多变量空间中的距离非常明显，除非移除其他两个支系，否则无法评估支系A内的任何结构。通过这样做，可以看出，正如Cyt-b树所描述的，Sailapata的样本与其他样本相距较远，解释了28.8%的变异（图2b）。其余样本在水平轴上几乎没有分离，但在y轴上看到了分化，解释了20.6%的变异。这些样本的分布似乎遵循从东部Sach’a Loma的样本到Pusuq’huni和Torotoro的西部样本的东-西梯度（图2a）。

**根据Evanno方法**，当K = 2且使用独立等位基因频率时，A. boliviensis的群体结构得到了更好的解释（Mean LnP = -10.202，SD = 14.1，ΔK = 342.5）。在这种情况下，一个群体由所有来自支系A的样本组成（推断祖先为1），另一个群体由支系C的样本组成（推断祖先为1）。支系B的样本显示出与支系A的混合，与支系C的混合较少（0.658 vs 0.342）（图2c）。在更高的K值下，LnP的值趋于平稳，SD增加，ΔK显著小于K = 2时的值。此外，根据我们对物种和系统的了解，将样本分配到其他群体变得更加困难。例如，在K = 3时，第三个群体被分配给CBF5455（祖先为0.382）和TK161718（祖先为0.572）的样本（图2c）。STRUCTURE HARVESTER生成的结果图作为补充材料提供。

**基于距离的隔离测试**
基于距离的隔离测试表明，A. boliviensis的结构遵循IBD模式，但其他生物学情景也可能解释其中的一部分。尽管散点图显示遗传距离和地理距离之间存在正相关，但距离矩阵之间的原始相关值（黑点）落在参考分布（直方图本身）内，表明基于距离的隔离并不显著（图3）。

**讨论**
对线粒体和基因组核区的分析表明，本研究中评估的Abrocoma boliviensis样本代表了三个不同的支系。这通过高遗传距离值（表1）、几乎没有混合（图2）以及谱系之间的长分化时间得到证实（即支系C大约在5.4 Ma分化，支系A和B大约在2.8 Ma分化）（Arenas-Viveros，2024年）。A. boliviensis内部的遗传结构，即支系A、B和C，在Cyt-b基因树以及使用GBS数据进行的PCA和STRUCTURE分析中都很明显（图2）。此外，我们发现这些遗传群落占据了具有不同生态和植物地理特征的地区：支系A和B出现在Yungas或干燥安第斯森林主导的地区，而支系C出现在Boliviano-Tucumano和干燥安第斯森林的镶嵌地带。总的来说，这表明这些谱系已经处于不同的进化路径上，因此至少应该被识别为不同的ESUs（进化单元）。也可以将每个支系提升到亚种的层次，或者将支系A和B提升到亚种水平，并将来自玻利维亚南部的个体（即支系C）归类为新物种。然而，我们采取了一个更为保守的立场，认识到仅凭遗传分化本身并不足以作为分类特征。此外，关于每个支系成员的形态学和生态学的比较研究尚缺乏，在这些信息被包括之前，很难对该物种的分类得出有充分支持的结论。例如，另一种安第斯高地特化物种vicuna（Vicugna vicugna）根据线粒体研究被划分为来自湿润和干燥Puna的亚种，但这种分类也得到了形态学和遗传特征的支持（Marín等人，2007年）。此外，还需要纳入来自支系B和C的额外数据，并确定Chuquisaca省（即支系B和C之间的地区）样本的缺失是否具有生物学意义，或者是采样偏差的结果。这一考虑很重要，因为该地区其他需要保护的物种，如安第斯猫（Leopardus jacobita）的分布自然上是碎片化的，建议在未来的保护行动中考虑到这一点（Cossíos等人，2012年）。虽然我们的结果得到了线粒体和核数据的支持，但在解释结果时必须考虑我们研究中样本数量较少的问题，特别是从人口统计分析来看。玻利维亚毛丝鼠的核苷酸多样性（π）和多态性位点的数量无疑受到是否将ESUs分组或单独评估的影响。仅包括来自Tarija的样本（即支系C）就使多态性位点的数量翻倍（121 vs 68），并且核苷酸多样性增加了36%（0.045 vs 0.033）（表2）。与两种高地特化物种的山地 Degu（Octodontomys gliroides，π = 0.0059和24个多态性位点）和安第斯猫（π = 0.0011–0.0027和15个多态性位点）相比，无论支系如何分组，A. boliviensis都具有更高的核苷酸多样性和多态性位点数量。这项研究很有趣，因为这两种物种的分布范围都非常广泛，涵盖了阿尔蒂普拉诺（Altiplano）的广大地理区域——安第斯猫（Andean Cat）的分布范围从秘鲁中部延伸到玻利维亚北部，而山地鼠兔（mountain degu）的分布范围则涵盖了智利和阿根廷（Cossíos等人，2012年；Rivera等人，2016年）。相比之下，A. boliviensis的已知分布范围仅限于玻利维亚东部。这表明，一个仅限于安第斯山脉某一个山坡的物种可能比那些分布在大陆更广阔区域的物种拥有更高的遗传多样性和多态性。然而，当分别计算A和B两个分支的数值时，核苷酸多样性（分别为π = 0.012和0.003）和多态性（分别为28和4）都显著下降。这表明每个进化支系（ESU）都携带着其他种群所没有的独特多样性，进一步支持了这样的假设：对其中任何一个分支的干扰都可能导致该物种发生严重的遗传损失；因此，需要针对它们分别采取监测、研究和保护措施（Borges等人，2018年）。

无论分组方式如何，Fu的Fs值和Tajima的D值均为正值，且只有当将A和B两个分支合并时，Fs值才具有统计学意义，这暗示了可能存在种群瓶颈效应或平衡选择的作用。从结果来看，A. boliviensis的种群可能经历了遗传瓶颈事件，结合其有限的分布范围，这凸显了进一步开展保护工作的必要性。遗憾的是，目前没有关于这些鼠兔类物种的长期种群趋势数据可供参考；尽管我们的研究受到样本规模的限制，但这些结果仍能大致反映它们的现状。

现存鼠兔类物种的进化历史可以追溯到大约1000万年前，当时安第斯山脉及其周边低地仍在经历造山运动，这些变化改变了地貌并促进了物种的分化（Gregory-Wodzicki，2000年；Garzione等人，2008年；Graham，2009年）。A. boliviensis的祖先在大约600万年前开始走上了自己的进化路径（Arenas-Viveros，2024年），当时东部山脉发生了第二次抬升（Garzione等人，2008年）。自大约1000万年前首次抬升以来，安第斯山脉的东坡就形成了其特有的起伏地形。随着山脉继续抬升，这种崎岖的地貌在中高海拔地区形成了物种迁移的障碍（Graham，2009年）。此外，更新世晚期从安第斯山麓延伸到巴西卡廷加斯（Caatingas）的连续分布的热带雨林因其他生境类型的出现而变得支离破碎（例如，A. boliviensis在玻利维亚南部的分布区域中的图库马诺-玻利维亚森林（Tucumano-Boliviano forests）（Mogni等人，2015年）。这些历史上的迁移障碍，加上对不同环境的适应（即北部和南部样本所处的不同生态区），可以解释A. boliviensis各分支之间的高分化程度以及几乎不存在基因交流的现象。除了历史过程外，玻利维亚鼠兔栖息地的广泛变化和开发（Brandt和Townsend，2006年；Killeen等人，2008年）也可能成为导致物种不连续性的原因；即使目前某些种群状况良好，但随着未来环境的变化，这种情况也可能发生变化。

本研究中描述的进化支系（即A、B和C分支）可能由于过去的过程而彼此隔离，并且可能因现代地貌的变迁而继续处于隔离状态。识别和保护这些进化支系的主要目的之一是为了维持遗传多样性，从而在面对环境变化时最大化物种的进化和适应潜力（Funk等人，2012年）。这在安第斯山脉东部尤为重要，因为除了人为因素导致的地貌变化外，过去植被的变化（主要是安第斯草原中的Polylepis林地）还受到温度、降水量和燃烧模式变化的影响（Williams等人，2011年），而这些因素都因气候变化而加剧。与其他生活在阿尔蒂普拉诺开阔环境中的鼠兔类物种不同，A. boliviensis主要与森林区域相关（Patton和Emmons，2015年），因此地貌的变化对其种群的未来构成了威胁，除非采取措施管理和保护其原生栖息地。

除了濒危和特有性之外，玻利维亚鼠兔还非常稀有，但这在其保护工作中可能具有优势。在保护多样性方面，稀有物种是选择保护和管理地点的最佳指标之一（Lawler等人，2003年）。事实上，专门为保护稀有物种而设立的保护区已被证明能更有效地保护生物多样性（Drever等人，2012年）。此外，稀有物种的存在和维护有助于维持群落和生态系统的功能结构（Mouillot等人，2013年；Leit?o等人，2016年）。所有这些因素，加上为A. boliviensis的保护工作赋予的经济价值（例如通过生态旅游和保护捐款），都有助于说服相关部门和机构开展系统的监测和研究计划（Drever等人，2012年）。

总之，基于本研究的结果以及对A. boliviensis了解的有限程度，我们强调了迫切需要加强调查工作。一旦获得更详细的种群数据，人口统计学、生态学和遗传学研究将有助于更好地了解该物种及其分支的进化轨迹，包括是否需要对其进行分类学上的重新调整。这些信息将有助于制定有效的管理和保护措施，包括确定监测区域、建立保护区及生态走廊，以及制定减少人类活动对地貌影响的缓解计划。