将群体基因组学与古代DNA相结合，以了解马达加斯加斑鸠的岛屿殖民历史

《Molecular Ecology》：Combining Population Genomics With Ancient DNA to Understand Island Colonisation History of the Madagascar Turtle Dove

　　摘要

马斯克林群岛（毛里求斯、留尼汪和罗德里格斯）的特点是人类首次到达的时间相对较近，这为研究近期物种殖民过程以及过去人类活动的影响提供了独特的背景。在本文中，我们结合了现代和古代DNA数据，以毛里求斯和留尼汪的马达加斯加斑鸠（Nesoenas picturata）为例，探讨了该物种的殖民历史及其对保护工作的影响。我们生成了一个参考基因组，并对现代N. picturata种群进行了重新测序，同时还获取了相关亚化石的全基因组数据。通过无模型推断、基于位点频率谱（SFS）的人口统计建模以及对种群结构（包括亚化石）的分析，我们发现N. picturata是在人类到达之前，从马达加斯加独立自然迁移到这两个岛屿的。统计结果和基于SFS的建模显示，岛屿上的种群具有较大的有效种群规模（Ne）和高遗传多样性，这与历史上关于人类导致的人口崩溃的描述相矛盾。根据拟合优度、遗传结构和多样性指标，两种最佳拟合模型无法区分：一种模型认为近期种群规模较大且迁移率较低；另一种模型则认为近期经历了严重的种群瓶颈，随后有大量个体从马达加斯加迁入。尽管如此，连锁不平衡现象更支持后一种情景，这也可能解释了高遗传多样性。现代和古代DNA数据均独立支持N. picturata原产于这两个岛屿的结论。我们的发现强调了利用多种数据来源以及多种统计方法来解析近期殖民历史的潜力。

1 引言

在过去的一个世纪里，全球化推动了人类介导的物种引入呈指数级增长（Bertelsmeier等人2025；Roy等人2024）。非本地物种被人类有意或无意地引入新区域后，常常通过捕食、竞争和栖息地破坏对本地生物多样性构成严重威胁（Freed和Cann 2009；Mooney和Cleland 2001）。因此，科学家和保护管理者经常反对这些物种的存在（Simberloff等人2013）。虽然防止非本地物种的引入是最有效的策略，但一旦它们建立起来，通常会采取根除措施来减轻其负面影响（Roy等人2024；Russell等人2017）。鉴于消除非本地物种所需的措施非常严厉，确定物种的生物地理起源对于自然环境的适当管理至关重要（Buckley和Catford 2016；Simberloff等人2012, 2013），尽管在近期时间和物种内部尺度上这可能具有挑战性。推断物种生物地理起源的关键方法是通过多种群体遗传学方法来测试特定的人口统计情景（Estoup和Guillemaud 2010）。这些方法使用不同的遗传数据分析工具，包括位点频率谱（SFS）、连锁不平衡（LD）、遗传身份和单倍型结构（Beichman等人2017；Marchi等人2021）。基于祖先重组图（ARG）的方法最近成为推断近期人口历史的强大工具，因为它们能够明确模拟单倍型谱系和重组过程。然而，这些方法通常需要分阶段的基因组和密集的采样，这些条件在非模型系统中难以满足（Brandt等人2024）。人口统计参数是通过最大似然法或贝叶斯推断从遗传数据中推断出来的，例如近似贝叶斯计算（Beaumont 2010；Beerli 2006；Estoup等人2004；Excoffier等人2013）。然而，参数估计可能不准确，特别是在关注近期世代时，而这正是研究非本地物种扩散时经常关注的对象。评估拟合优度是确保人口统计模型生物学和时间相关性的重要步骤，这需要大量的事后模拟以及整合多样的基因组信息（Johri等人2021）。通过确认过去已知时间点的谱系遗传组成，定年的古代DNA序列可能为检查近期时间尺度上的推断提供补充方法，尤其是在新环境殖民方面（Leonardi等人2017）。海洋岛屿群岛的独特特征，如地理隔离、相对较小的规模和年轻的年龄，有助于推断谱系历史（Losos和Ricklefs 2009；Martin等人2021；Warren等人2015）。这些特征结合详细的历史、生态和古生物学记录，为评估人口统计模型的生物学相关性提供了理想的背景，使岛屿群岛成为开发跨时间尺度综合遗传学方法的天然实验室。尽管全基因组数据的可用性不断增加，但在大多数情况下，解析近期的人口历史（包括人类活动的影响）仍然是群体基因组学的挑战（Santiago等人2020）。在这方面，马斯克林群岛作为马达加斯加生物多样性热点地区的一部分（Myers等人2000；Schlüter 2016；Thébaud等人2009），提供了一个特别有价值的案例研究。与其他生物多样性丰富的群岛不同，该群岛在大约400年前人类首次定居之前是无人居住的。尽管该群岛由火山形成（包括毛里求斯、留尼汪和罗德里格斯三个岛屿），但它拥有丰富的亚化石沉积物（Hume 2005），这些沉积物涵盖了人类首次定居前后的时期。这些沉积物记录了自人类到达以来本地动物群所遭受的严重压力，表明50%–60%的本地脊椎动物已经灭绝，其中包括渡渡鸟（Cheke和Hume 2008；Probst和Brial 2002）。与许多受到非本地物种严重影响岛屿系统类似（Russell等人2017），马斯克林群岛的物种灭绝部分是由引入的物种造成的（Cheke 1987；Cheke和Hume 2008），其中鸟类受到的影响尤为严重（17种特有物种受到威胁；Hume 2013）。保护剩余马斯克林鸟类的另一个重要问题是，虽然不同脊椎动物的杂交倾向各不相同，但在鸟类中杂交率很高（Grant和Grant 1997；Warren等人2012）。例如，在邻近的塞舌尔群岛，本地塞舌尔斑鸠（Nesoenas picturata rostrata）与引入的N. picturata之间的杂交导致本地物种具有与入侵者相同的表型（Rhymer和Simberloff 1996）。N. picturata原产于马达加斯加，与濒危的毛里求斯粉鸽（N. mayeri）关系密切（Cheke 2005；Hume 2011；Johnson等人2001）。来自毛里求斯和留尼汪的亚化石证据表明，这些岛屿上曾经存在与N. picturata相关的本地种群。一些具有较完整骨盆元素的亚化石被用来提出一种已灭绝的毛里求斯斑鸠——N. cicur的存在（Hume 2011）。然而，N. picturata在整个群岛的殖民历史仍不清楚。关于现存种群是否原产于毛里求斯和留尼汪，或者该物种是在人类定居后随着原有本地种群的灭绝而重新引入的，存在争议。尽管N. picturata的殖民历史存在不确定性，但在毛里求斯和留尼汪的管理措施有所不同。在毛里求斯，该物种被视为入侵物种，并在其与粉鸽重叠的区域被捕杀。在留尼汪，自1990年代经历种群瓶颈后，该物种自1983年起受到保护（Barré 1983；Probst 2002）。目前，留尼汪的保护主义者面临来自狩猎团体的压力，要求重新考虑其保护地位。我们生成了来自马达加斯加、毛里求斯和留尼汪的现代个体的全基因组，以及毛里求斯的亚化石的全基因组数据。我们研究了种群结构、遗传多样性和纯合度（ROH），并提出了使用适合观察到的SFS的共祖框架推断出的具体人口统计情景。我们通过事后建模预期的连锁不平衡模式来验证我们最可能的情景，这一模式在推断过程中并未使用。我们的主要目的是通过估计马斯克林群岛的殖民路线和时间，来确定N. picturata是原产于还是引入到毛里求斯和留尼汪的。其次，我们研究了陆地之间迁移对现代N. picturata种群的贡献，区分了人类出现前后的迁移情况。第三，我们研究了人类活动对马斯克林群岛种群遗传多样性的影响。理解过去人类活动在物种分布中的作用是全球保护生物学中反复出现的问题（Gallardo等人2015）。我们通过比较古代DNA和现代基因组DNA来源的结论，提供了一个解析近期殖民历史的案例研究。前者使用已知年龄的亚化石中的古代DNA，而后者结合了历史记录、多样的现代基因组统计和群体遗传学建模。因此，我们的结果不仅对马斯克林群岛的保护组织具有直接的管理策略意义，还展示了整合多种独立分子生态证据揭示近期物种殖民历史的潜力（García-Verdugo等人2019；Theissinger等人2023）。

2 方法

2.1 参考基因组、现代基因组采样和测序

我们从一个2022年在Plaine Lièvre（毛里求斯，Black River Gorges国家公园）采集血液样本的个体中组装了一个新的N. picturata基因组。有关DNA保存和提取、长读长PacBio测序和基因组组装的详细信息，请参见S1。我们的现代重新测序基因组样本涵盖了毛里求斯（N=9）和留尼汪（N=10）的整个物种分布范围，以及马达加斯加的三个地理上不同的地区（N=3），总共22个基因组（图1）。在毛里求斯（2020年和2022年）和留尼汪（2022年）现场采集了12个个体样本。组织样本来自Field Museum of Natural History和Réunion鸟类康复中心（表S1）。有关DNA提取和基因组重新测序的详细信息，请参见S1。

2.2 亚化石采样和测序

虽然大多数亚化石可以鉴定到属级别，但从许多骨骼类型的形态学特征很难区分N. picturata和N. cicur（Hume 2011）。因此，从一开始，DNA测序就有潜力确认物种身份。为了最大化采样两种物种的机会，同时尽量减少对珍贵收藏的损害，我们获得了四个被形态学鉴定为N. picturata或N. cicur的Nesoenas亚化石，并在获得许可的情况下进行了破坏性采样以提取DNA（表S2）。其中三个来自J.P. Hume的收藏，另一个来自MNHN古生物学收藏。一个亚化石来自缺氧沼泽（Mare aux Songes遗址），另外三个来自洞穴和巨石堆（Le Pouce和Le Morne遗址；表S2）。所有古代DNA分子工作都在图卢兹的人类生物学和基因组学中心的古代DNA设施中进行，遵循严格的实验标准以避免和监测污染，包括使用一次性个人防护装备、正压环境、漂白/紫外线表面和仪器消毒以及阴性对照空白。PCR扩增、纯化和定量以及DNA捕获和测序准备都在与古代DNA设施物理隔离的建筑物内的后PCR实验室中进行。使用钻头轻轻打磨表面以去除潜在污染物后，最多50毫克的骨头被Mixel Mill 200（Retsch）微粉碎机粉碎。有关古代DNA方法的详细信息，包括提取、霰弹库制备、HyRAD探针开发、杂交捕获和测序，请参见S2。

2.3 现代样本的读取映射和变异调用

原始读取数据使用Trimmomatic v.0.39处理，去除适配子和低质量片段，使用默认设置，除了最小读取长度为100个碱基对（bp）（Bolger等人2014）。剩余的读取数据使用bwa-mem v.0.7.17映射到参考基因组（Li和Durbin 2009）。对齐结果使用SAMtools v.1.14进行排序和索引（Danecek等人2021）。使用Picard的MarkDuplicates工具（https://broadinstitute.github.io/picard/）移除了PCR重复序列。变异调用在32个超过10 Mbp的支架上进行，占基因组的83%，遵循GATK最佳实践，使用GATK v4.2.6.1（McKenna等人2010）。BAM文件使用-ERC BP_RESOLUTION标志输入到HaplotypeCaller中（DePristo等人2011；Van Der Auwera等人2013）。结果gVCF文件使用GenomicsDBImport合并，并使用GenotypeGVCF和–include-non-variant-sites标志共同调用，以便调整遗传多样性（DePristo等人2011；Van Der Auwera等人2013）。我们使用VCFtools v.0.1.16计算每个位点/个体的平均深度、质量、缺失率和近交系数，以确保数据完整性（Danecek等人2021）。使用bedtools v.2.30 intersect（Quinlan和Hall 2010）移除了重复元素（见S1进行识别）。请参阅S3部分，了解VCF文件的过滤方法。

2.4 亚化石样本的读取映射和变异调用

由于古代DNA的高度片段化特性、死后积累的碱基修饰以及污染风险，我们对亚化石样本采用了不同的生物信息学流程。原始读取数据使用AdapterRemoval v.2.3.3（Schubert等人，2016年）进行了修剪和过滤。读取数据的两端包含Illumina适配器和7个碱基对的内部条形码。通过设置最小适配器重叠长度为3个碱基对（bp）和最大错配比例为1/5来修剪Illumina适配器。内部条形码通过设置–trim3p和–trim5p标志为7 bp来修剪。我们移除了长度小于25 bp的读取数据以及含有不确定碱基的低质量片段。分别处理修剪后的双端和单端读取数据，以最大化映射到的读取数量。在这两种情况下，都使用bwa-aln算法将读取数据映射到参考基因组上，最小错误率为0.1，种子长度为30（Li和Durbin，2009年）。生成的SAM文件使用bwa-sampe和bwa-samse（Li和Durbin，2009年）进行合并。移除了长度小于30 bp且映射质量得分低于30的读取数据。使用markdup工具（Samtools v.1.18；Danecek等人，2021年）移除了PCR重复序列。双端和单端BAM文件使用Samtools v.1.18 merge（Danecek等人，2021年）进行合并。使用mapDamage v.2.2.2和PMDtools v.0.60（Jónsson等人，2013年；Skoglund等人，2014年）评估DNA损伤模式。由于亚化石序列之间的重叠非常少，我们分别计算了每个亚化石样本与现代样本的基因型可能性。基因型可能性是使用ANGSD v.0.935中实现的SAMtools（?GL 1）方法计算的（Korneliussen等人，2014年）。只保留了显著的多态性位点（p < 10^-7）。测试表明，在较高深度下观察到的模式仍然保持一致。因此，为了确保我们在亚化石样本中调用到具有足够覆盖率的位点，我们将–minInd标志设置为个体总数，将–setMinDepthInd标志设置为3。

2.5 种群结构

2.5.1 现代样本

我们使用R语言中的pcadapt包v.4.3.3进行了主成分分析（PCA）来研究种群结构（Privé等人，2020年）。对于每对种群（毛里求斯-留尼汪；马达加斯加-毛里求斯；马达加斯加-留尼汪），我们使用基于Reynold估计器的自定义R脚本计算了成对FST（Reynolds等人，1983年），通过1000次排列来评估显著性，在这些排列中随机化个体采样点，并计算大于观察值的重采样值的数量。我们还使用了LEA R包中实现的sNMF聚类算法（Frichot和Fran?ois，2015年）。该程序推断祖先种群的数量K并估计个体祖先比例。我们对每个K值进行了10次重复实验，范围从1到8，并研究了交叉熵曲线以确定最可能的解决方案。

2.5.2 亚化石样本

我们使用PCAngsd v.0.99（Meisner和Albrechtsen，2018年）根据之前计算的基因型可能性对亚化石和现代样本进行了PCA。同样，我们使用ANGSD中实现的身份-状态系数（Meisner和Albrechtsen，2018年）计算了每个亚化石和现代样本的距离矩阵。

2.6 现代样本的遗传多样性和ROH

使用Python脚本（https://github.com/PierreLesturgie）从经过人口统计学过滤的VCFs中计算了位点频率谱，从中得出了核苷酸多样性（θπ）（Tajima，1983年）、Watterson的θ（θw）（Watterson，1975年）和Tajima的D（Tajima，1989年），并通过总位点数对θπ和θw进行了归一化，以获得每个位点的遗传多样性估计值。这些计算是在整个基因组上以及对于每个单独的支架上进行的。我们还从每个个体的杂合位点总数中计算了全基因组的个体观察到的杂合性（HO）。使用PLINK v.1.90中的–homozyg函数估计了ROH。该函数使用滑动窗口来识别连续的纯合片段。有关使用的参数和阈值，请参阅S4部分。

2.7 人口统计推断

我们通过比较两种已知随着进化时间推移而表现不同的方法，重建了现代样本的共祖率变化（Spence等人，2018年；Terhorst等人，2017年）。成对顺序马尔可夫共祖（PSMC v.0.6.5）使用顺序马尔可夫共祖方法，只需要一个二倍体基因组就可以推断出给定时间点的泛混合种群的有效种群大小（Ne）（Li和Durbin，2011年）。当不能假设泛混合时，PSMC仍然是一个有用的工具，可以提供关于种群结构或分化的信息（Mazet等人，2016年）。我们使用BCFtools v.1.11 mpileup及其–c标志在参考基因组上为每个个体创建了一个共识序列。有关输出转换为FASTA、相关过滤、区间大小和自助法的信息，请参阅S5部分。为了运行PSMC v.0.6.5，我们将–t标志设置为40，–r标志设置为4，–p标志设置为‘95×2+4+6’。在Mather等人（2020年）的测试提示后，我们改变了默认的时间向量（?p），因为它减少了最近时间间隔内重组事件数量较少所导致的偏差。SMC++ v.1.15.2程序也使用顺序马尔可夫共祖方法来估计共祖时间的分布，但利用SFS在最近的时间段内获得更好的准确性（Spence等人，2018年；Terhorst等人，2017年）。有关过滤和文件转换的信息，请参阅S5部分。使用estimate命令计算复合似然性和随时间变化的共祖率，将–knots标志设置为40。对于所有历史人口统计分析，我们使用了3.7年的世代时间（Seal和Bruford，1991年）和每个世代每个位点5.2×10^-9的替代率（Shapiro等人，2013年）。由于缺乏N. picturata的具体估计值，并且考虑到鸟类的基因组特性通常是保守的（Kawakami等人，2014年；Singhal等人，2015年；Zhang等人，2014年），我们使用了来自同属物种粉红鸽（Nesoenas mayeri）和家鸽（Columba livia）的研究得出的值。粉红鸽是N. picturata的姐妹物种，它们具有相似的生活史特征（Johnson等人，2001年；Shapiro等人，2013年）。

2.8 人口统计建模

我们使用fastsimcoal2 v.2709（Excoffier等人，2013年；Marchi等人，2021年）进行了基于共祖的最大似然（ML）分析，以推断N. picturata的历史人口统计最可能的情景。该方法利用成对二维SFS（2D-SFS）来估计给定人口统计情景的参数。我们从经过人口统计学过滤的VCF中为每对种群生成了折叠的2D-SFS。我们测试了八个模型，逐步增加复杂性，并将每个模型与之前最可能的模型进行比较。首先，我们使用六个模型测试了几种种群树拓扑结构和殖民序列（图S1）。为了缩小模型测试的范围，我们考虑了地理因素和种群结构分析结果。具体来说，由于岛屿之间的距离和观察到的遗传结构，除了一个树拓扑模型外，所有模型都将马达加斯加设置为祖先群体。我们考虑了六个模型：（i）MOD1，马达加斯加的种群是由毛里求斯和留尼汪的种群混合而成的，以测试sNMF聚类是否可以由这一情景解释；（ii）MOD2，马达加斯加独立殖民，毛里求斯在留尼汪之前殖民；（iii）MOD3，马达加斯加独立殖民，留尼汪在毛里求斯之前殖民；（iv）MOD4，依赖性殖民，留尼汪的种群起源于毛里求斯；（v）MOD5，依赖性殖民，毛里求斯的种群起源于留尼汪；（vi）MOD6，马达加斯加独立殖民，殖民事件没有固定时间。其次，我们通过将不对称迁移纳入最可能的树拓扑结构中，创建了第七个模型MOD6_MIG。第三，我们评估了毛里求斯和留尼汪的Ne变化，以及最近从马达加斯加到两个岛屿的迁移，从而得到了最终模型MOD6_MIG_TR。引入这些额外参数是基于人口统计推断和历史记录，这些记录表明N. picturata个体从马达加斯加到毛里求斯和留尼汪的迁移发生了变化（Cheke和Hume，2008年；Probst，2002年）。每个模型独立运行了100次（最终最可能的模型运行了200次），每次都有40次期望条件最大化（ECM）循环和100,000次共祖模拟用于SFS估计。有关命令和每个参数指定范围的详细信息，请参阅表S3。我们使用赤池信息量准则（AIC）评估了模型。我们还确保了在ML参数下100次模拟的似然分布不与之前模型的分布重叠。对于最可能的模型，我们进行了参数自助法，以计算两个感兴趣的运行的ML参数的95%置信区间。对于每次运行，这个过程涉及模拟模型100次以生成2D-SFS，对每个模拟的2D-SFS进行100次独立参数估计（如同对实证数据所做的那样），并使用得到的参数分布来计算5%和95%的分位数。

2.9 最可能模型的拟合优度

2.9.1 模拟结构和多样性

在最可能的人口统计模型中，我们识别并保留了两个对比鲜明的运行结果，以评估拟合优度。我们通过使用每次运行的ML参数模拟基因组中最长的支架来独立评估每个运行的拟合优度，以估计不确定性。我们设置了每代的重组率为10^-8。使用fastsimcoal2中的以下命令选项模拟了基因型表：–n 100 –I –G –g –s0 –x –k 10,000,000。使用自定义的bash脚本，我们将基因型表转换为VCFs，并：（i）计算了折叠的SFS，从中获得了θπ、θw和Tajima的D；（ii）过滤掉了MAF低于0.05的SNPs，并通过PCA、sNMF聚类和成对FST探索了种群结构。

2.9.2 逆瞬时共祖率（IICR）和LD统计

我们使用fastsimcoal2和recordMRCA标志，在给定的人口统计参数下，模拟了一对单倍体染色体（T2）的107个独立的最近共同祖先（TMRCA）时间值。我们使用最可能模型的每个参数自助法的ML参数重复了这个过程。我们建立了一个包含600个对数间隔时间点的向量，时间范围从20到4.5×10^5代之前，使用自定义的python脚本获得了T2的累积密度函数F(T2)。有关如何获得IICR和LD衰减的信息，请参阅S6部分。我们遵循Ragsdale和Gravel（2019）的建议，将理论LD衰减与观察到的LD衰减进行比较（https://momentsld.github.io/moments/）。有关LD统计的计算、观察数据的置信区间和LD衰减的期望值的信息，请参阅S6部分。

3 结果

3.1 数据质量和过滤

3.1.1 现代样本

整个数据集的平均测序深度为23×，对于长度超过10 Mbp的支架，这一深度增加到36×，每个样本的原始读取数量在253,784,758到639,821,009之间（表S1）。过滤后，保留了16,466,416个SNPs用于种群结构分析。人口统计学过滤后，每个种群保留了10,897,013到12,739,591个SNPs。

3.1.2 亚化石样本

样本的内源DNA含量在捕获前介于0.03%到1.5%之间，捕获后介于0.11%到5.8%之间（表S2）。每个亚化石的映射读取数量在31,024到2,719,895之间（表S1）。MiniSeq鸟枪法测序的DNA链长度分布在35到55之间，这与亚化石的情况相符。使用mapDamage2和PMD-tools对未经过USER处理的文库进行测序后得到的真实性谱没有显示出显著的DNA损伤，这对于亚化石来说是典型的。亚化石序列之间的重叠部分，在3×的深度下仅为1030 bp。每个亚化石与现代数据集共享的SNPs数量在2811到680,209之间（表1）。只有在确认使用了足够的亚化石骨骼来获取所需的DNA后，才能对其他亚化石（50-100 mg）进行放射性碳测年（由牛津放射性碳加速器单元进行）。经过校准后，Sub1的年代为903-1024 calAD，Sub4为1691-1954 calAD，Sub13为1689-1954 calAD（表S2）。由于从现场返回时使用了硬化剂保护亚化石，Sub3未能获得年代数据。表1显示了在≥3×深度下亚化石和现代样本之间的SNPs数量。

3.2 种群结构

PCA显示了三个不同的簇，分别对应于来自三个采样岛屿的个体：马达加斯加、毛里求斯和留尼汪（图2A和S2）。PC1解释了数据集中16.7%的总变异，明确地将三个岛屿区分开来，表明马斯克林群岛内部（毛里求斯岛和留尼汪岛之间）的差异比任何一个马斯克林岛屿与马达加斯加之间的差异都要大。毛里求斯岛内的变异在PC1和PC2上都有体现。基于p距离的层次聚类得出了类似的结果（图S3）。与PCA一致，最高的成对FST值出现在留尼汪岛和毛里求斯岛之间（FST=0.15，p=0.001；图2B）。马达加斯加和毛里求斯岛之间的FST值略高（FST=0.11，p=0.006），而马达加斯加和留尼汪岛之间的FST值较低（FST=0.09，p=0.004；图2B）。

使用N. picturata基因组进行的人口结构分析：(A) PCA。轴代表前两个主成分（PC）。标签标识了样本205E和225E。(B) 人群间Reynold's FST的热图。星号表示p<0.05。使用sNMF进行聚类显示马达加斯加斑鸠种群之间存在显著的结构，K=2时能最好地解释数据中的结构（图S4）。当每个种群的样本量相等时（图S5），或者使用其他聚类方法时（图S6和S7；详见S7节），这一结果仍然成立。根据K=2时的sNMF聚类，来自毛里求斯的个体和留尼汪岛的个体具有最高的混合血统比例。我们包括K=3，因为PCA显示有三个簇。在sNMF下，虽然留尼汪岛的个体具有相同的血统，但毛里求斯的个体表现出混合血统，这与PCA分析中观察到的较大变异性一致（图3A）。

使用N. picturata基因组的sNMF聚类：(A) K=2（对应最低的交叉熵）和(B) K=3。标签标识了样本205E和225E。

3.2.2 次化石样本
使用次化石和现代样本的数据进行了PCA分析（图4；见表1了解变异数量）。Sub1和Sub13分别被描述为N. cicur或N. picturata（表S2），它们与现代毛里求斯谱系聚类在一起。这些次化石与现代个体之间的成对IBS距离也表明它们与毛里求斯种群具有高度遗传相似性（图S8和S9）。Sub3和Sub4在形态上被鉴定为N. cicur或N. picturata（表S2），它们更接近现代马达加斯加谱系。成对IBS距离显示Sub3和Sub4在遗传上与所有现代谱系都不同，并不显示出与马达加斯加个体的较短距离（图S10和S11）。

3.3 现代样本的遗传多样性和ROH
全基因组的估计表明，马达加斯加的遗传多样性高于毛里求斯和留尼汪（表2），且它们的多样性分布不重叠（图S12）。火山岛种群的Tajima's D值为正，表明高频变异在中性常数人口模型下有轻微的过度代表（表2）。相反，马达加斯加的Tajima's D值为负，表明低频变异的过度代表程度也略高于中性模型预期。

通过计算平均个体全基因组的HO值，我们发现小岛屿种群的遗传变异较低，毛里求斯和留尼汪的平均值为每kb 5.74和5.68个杂合子（het/kb）（图5）。相比之下，马达加斯加种群的平均值为6.82个杂合子/kb。火山岛种群的全基因组HO方差较大，尤其是在毛里求斯，其范围在5.42到5.88个杂合子/kb之间（图5）。

3.4 人口统计推断
PSMC和SMC++得出的种群动态曲线不同（图7），我们通过突变率对其进行缩放以提供大致的时间框架。火山岛种群中的IICR动态在时间和幅度上相似。PSMC表明三个岛屿的个体具有相似但不重叠的IICR轨迹，而SMC++明确表明马达加斯加至少在10万年前就从毛里求斯和留尼汪分离。如果用Ne变异来解释，SMC++的推断显示毛里求斯和留尼汪在1万到10万年前经历了急剧下降，随后是快速扩张期，进入稳定阶段。基于单基因组的PSMC推断突出了毛里求斯的两个个体（205E和225E）的不同趋势。这两个个体在sNMF和ROH分析中也可以识别出来，但在种群内部结构分析中它们并未从自己的样本组中分离出来（图S15）。

3.5 人口建模
我们使用fastsimcoal2测试了八种人口统计情景。首先测试了MOD1到MOD6的树拓扑和殖民序列（图S1）。MOD2被认为是最可能的模型（表S3）。该模型涉及两次独立的从马达加斯加到马斯克林群岛的殖民，其中序列是固定的，第一次是到毛里求斯，第二次是到留尼汪。两个岛屿的估计殖民时间相似，分别为119,036年和110,612年前。当比较MOD2和MOD6时，尽管MOD2的序列是固定的，但它们的似然分布重叠（图S16）。鉴于殖民时间的相似性以及MOD2和MOD6之间的似然差异可以忽略不计，我们保留了MOD6作为最可能的树拓扑，因为它允许独立估计殖民事件。接下来，我们在MOD6_MIG模型中加入了不对称基因流。正如预期的那样，允许基因流显著增加了模型的似然性，并降低了AIC（图S16）。最后，我们在MOD6_MIG_TR中增加了模型的复杂性，以考虑：(i) 使用SMC++推断的毛里求斯和留尼汪的合并率变化，这个时间足够久远，可以排除近期时间偏差；(ii) 历史证据表明，最近从马达加斯加到毛里求斯和留尼汪的个体迁移很常见（Cheke和Hume 2008；Probst 2002）。MOD6_MIG_TR的AIC最低，且似然性显著高于其他任何模型（表S3）。我们运行了200次MOD6_MIG_TR模型并分析了结果参数分布。这项分析揭示了两种不同的情景：一些运行估计毛里求斯和留尼汪的当前种群规模较大，迁移率较低（NMAU和NREU），而其他运行则表明种群规模较小，迁移率较高（图S17和S18）。因此，虽然一些运行表明当前种群规模较大，但其他运行则表明发生了显著的近期瓶颈效应。估计较大种群规模的运行倾向于具有较高的似然性，这与历史记录表明的从马达加斯加到毛里求斯和留尼汪的频繁迁移以及近期种群下降相矛盾（Cheke和Hume 2008；Probst 2002）。我们展示了两个示例运行，它们体现了当前种群规模和迁移率之间的对比关系，分别命名为HNLT和LNHT（图8）。

最可能的人口统计情景（MOD6_MIG_TR）的表示，按参数估计进行了缩放：(A) HNLT；(B) LNHT。在HNLT中，NMAU和NREU的估计值分别为102,403和43,602；而在LNHT中，这些估计值降至540和629（表3）。从殖民到现在，HNLT和LNHT中的Ne变化模式明显不同。在HNLT中，Ne在T3和T4时间点后减少，而在LNHT中则相反，在较近的过去增加。这些近期变化发生在T5和T6时间点，HNLT估计为207和366年前，LNHT估计为307和370年前，这些时期与马斯克林群岛的人类定居时间相符。

3.5 人口建模
我们使用fastsimcoal2测试了八种人口统计情景。首先测试了MOD1到MOD6的树拓扑和殖民序列（图S1）。MOD2被认为是最可能的模型（表S3）。该模型涉及两次独立的从马达加斯加到马斯克林群岛的殖民，其中序列是固定的，第一次是到毛里求斯，第二次是到留尼汪。两个岛屿的估计殖民时间分别为119,036年和110,612年前。当比较MOD2和MOD6时，尽管MOD2的序列是固定的，但它们的似然分布重叠（图S16）。鉴于殖民时间的相似性以及MOD2和MOD6之间的似然差异可以忽略不计，我们保留了MOD6作为最可能的树拓扑，因为它允许独立估计殖民事件。接下来，我们在MOD6_MIG模型中加入了不对称基因流（MOD6_MIG_MIG）。正如预期的那样，允许基因流显著增加了模型的似然性，并降低了AIC（图S16）。最后，我们在MOD6_MIG_TR中增加了模型的复杂性，以考虑：(i) 使用SMC++推断的毛里求斯和留尼汪的合并率变化，这个时间足够久远，可以排除近期时间偏差；(ii) 历史证据表明，最近从马达加斯加到毛里求斯和留尼汪的个体迁移很常见（Cheke和Hume 2008；Probst 2002）。MOD6_MIG_TR的AIC最低，且似然性显著高于其他任何模型（表S3）。我们运行了200次MOD6_MIG_TR模型并分析了结果参数分布。这项分析揭示了两种不同的情景：一些运行估计毛里求斯和留尼汪的当前种群规模较大，迁移率较低（NMAU和NREU），而其他运行则表明种群规模较小，迁移率较高（图S17和S18）。因此，虽然一些运行表明当前种群规模较大，但其他运行则表明发生了显著的近期瓶颈效应。估计较大种群规模和较低迁移率的运行倾向于具有较高的似然性，这与历史记录表明的从马达加斯加到毛里求斯和留尼汪的频繁迁移以及近期种群下降相矛盾（Cheke和Hume 2008；Probst 2002）。我们展示了两个示例运行，它们体现了当前种群规模和迁移率之间的对比关系，分别命名为HNLT和LNHT（图8）。两次运行都表明存在基因流动，迁移率在HNLT中范围为1.88×10^-6（mREU?MAU）到4.17×10^-5（mMAU?REU），在LNHT中范围为1.56×10^-5（mMAD?MAU）到8.28×10^-5（mREU?MAD）。在这两次运行中，迁移时间的估计都对应于人类在马斯克林群岛定居的时期。从马达加斯加到毛里求斯的迁移（T7）在HNLT中被估计发生在17万年前，在LNHT中被估计发生在27万年前。同样，从马达加斯加到留尼汪的迁移（T8）在HNLT中被估计发生在29.6万年前，在LNHT中被估计发生在33.3万年前（表3）。

3.6 最可能模型的拟合优度

3.6.1 模拟结构和多样性
为了确定两个MOD6_MIG_TR情景中哪一个最合理，我们通过使用每次运行中100次参数自助法获得的ML参数估计值，模拟了一个10,934,710个碱基对的支架（支架1的长度），来检验HNLT和LNHT的拟合优度。两次运行都准确地重建了观察到的种群结构。来自HNLT和LNHT的模拟的PCA一致地显示了三个种群簇（图S19）。与观察到的成对Fst值一致（图2），最高的估计分化发生在毛里求斯和留尼汪之间（两次运行中均为0.15；图S20）。sNMF簇在HNLT和LNHT之间几乎相同（图S21）。与观察数据一样，当K=3时，种群并没有分成不同的组。此外，马达加斯加的种群始终显示出比火山岛更高的遗传多样性（图S22）。观察到的和模拟的基因组平均统计量之间的微小差异，特别是在Tajima's D中，可以归因于支架之间和支架内部的固有变异（图S12；特别是对于支架1，图S23），这很可能反映了基因组的复杂性（这很难通过模拟来再现）。总的来说，两次运行都成功地再现了观察到的遗传结构和多样性模式，使得这些分析无法在合理性方面区分两次运行。

3.6.2 IICR和LD统计量
我们进一步使用HNLT和LNHT的参数自助法模拟了IICR（图9）。马达加斯加的IICR在两次运行中都相对稳定且较高。相比之下，毛里求斯和留尼汪的IICR轨迹在两次运行中显示出显著差异。在HNLT中，IICR表明在整个岛屿历史中种群规模保持稳定，除了大约1万年前的一次祖先迁移。在LNHT中，毛里求斯和留尼汪的IICR曲线显示出大约50万年前的明显瓶颈，导致当前种群规模减少。LNHT的IICR显示，与HNLT相比，毛里求斯和留尼汪在最近时期的方差更大，这与瓶颈强度的估计一致（图S24）。当我们使用PSMC和SMC++重新分析HNLT和LNHT的模拟数据时，得到了类似的轨迹，但LNHT的最近瓶颈没有被这两种算法中的任何一种检测到（图S25）。

4 讨论
非本地物种的扩散是全球保护关注的问题（Simberloff等人，2013年）。设计有效的保护策略来管理非本地物种通常受到对物种是本地物种还是非本地物种的不确定性阻碍。我们的研究利用了已知的人类首次到达马斯克林群岛的时间（400年前，Cheke和Hume，2008年），以及现代和古代DNA的测序，来研究马斯克林群岛上一种受保护计划关注的鸟类的状况，包括人类活动在最近过去的贡献。我们使用该数据集来分析遗传结构和多样性，估计ROH，并推断人口历史，发现了支持马达加斯加龟鸽殖民马斯克林群岛的具体情景的多个证据（表4）。我们进一步使用全面的人口遗传学方法来解决人口建模中的歧义。我们的关键发现，在地方层面上具有重要的保护意义（见S9），即N. picturata是毛里求斯和留尼汪的本地物种。两条独立的证据表明，N. picturata在人类到达之前很久就已经在马斯克林群岛建立了种群，表明其殖民是自然的，而不是由人类介导的。表4。支持N. picturata种群历史的证据一致性，包括具有多重支持的结论（+），可能解释数据的推断，以及对（+）和反对（-）的矛盾信号的解释。

结论可能解释数据
解释具有矛盾信号的支持和反对

支持的种群历史
情景
情景

1
2
3
4
5

马斯克林群岛的本地物种
马斯克林群岛的两次独立殖民
最近从马达加斯加的迁移
人类增加了马斯克林群岛种群的遗传多样性
高迁移率到马斯克林群岛，以及人类到达后先增加后减少的近期瓶颈（LNHT）

分析

PCA和成对IBS距离的定年化石
+
fastsimcoal2、模型的似然性和参数的ML
+
+
+
+
+
SMC++
+
+
+
使用sNMF进行聚类
+
+
+
使用ADMIXTURE进行聚类
+
+
+
PCA和热图
+
+
观察到的与理论LD衰减的比较
+
+
使用PLINKa计算的ROH

注：下面呈现的替代种群历史包括：1. 首次引入马斯克林群岛；2. 马斯克林群岛的单一殖民；3. 群体分化后没有岛屿谱系之间的混合；4. 人类减少了马斯克林群岛种群的遗传多样性；5. 低迁移率到马斯克林群岛，以及人类到达后先减少后增加的近期种群增长（HNLT）。

结论和推断——3和4——通过PLINK估计的ROH得到了支持。使用bcftools和不同的参数设置来估计ROH，既不支持也不反对3和4，尽管推断出的ROH的总体模式是相同的（S8，图S13）。首先，来自毛里求斯的两个Nesoenas亚化石与现代N. picturata种群聚类在一起，其中一个（Sub1）来自Mare aux Songes，经放射性碳测定为903-1024 calAD（表S2），表明现存的毛里求斯种群至少有1000年的历史，早于人类到达。鉴于Sub1和Sub13的PCA分别基于3911和21,794个变异（表1），与现代种群的聚类不能通过样本数量不足或缺乏变异来解释。虽然Sub13的放射性碳日期不能证实人类到达之前的推断，但它也不与其矛盾。我们注意到，如果我们从人口推断中得到相反的结果，即现存的N. picturata种群是完全新引入的，那么Sub13将提供关于引入事件时间的重要信息。我们还识别出在PCA中显示出介于毛里求斯和马达加斯加当前种群之间的中间分化的亚化石（Sub3和Sub4）。某些亚化石与现代种群的分化可能代表了一个已灭绝的、不同的谱系，如N. cicur，或者可能是由于基因组覆盖度低或年龄相关偏差导致的基因型置信度不足。正如预期的那样，从仅30-50毫克来自保存质量低的热带环境（沼泽、洞穴和巨石堆）的亚化石骨骼中获取的全基因组数据是具有挑战性的。遇到的最大挑战是，尽管优化了用于杂交的温度，我们获得的捕获特异性仍然很低，因此样本之间的重叠也很低。尽管如此，捕获仍然成功地将内源性DNA的产量提高了五倍，并在现代和古代样本之间获得了多达680,000个共享变异。因此，我们的研究表明，古代DNA可以确认非模型殖民物种在100-1000年的时间尺度上的身份，这段时间跨越并超过了人类干预发生的时期，而在这期间人口推断通常具有挑战性。其次，从现代基因组中进行的人口推断对于确定殖民和迁移的时间至关重要。多次运行最可能的基于共祖的模型一致估计出毛里求斯和留尼汪的殖民时间早于人类到达，即使考虑了置信区间（表3），为N. picturata是马斯克林群岛的本地物种提供了第二个独立的证据（表4）。在HNLT中，殖民事件估计发生在数十万年前，而LNHT的估计接近一百万年前。此外，推断的共祖率轨迹（特别是来自SMC++的）表明，火山岛种群在人类首次到达马斯克林群岛之前很久就已经从马达加斯加分化（图7）。推断的IICR取决于潜在的人口历史，并可以反映种群分化（Lesturgie等人，2022年）。我们的结果还展示了殖民事件的性质，三项分析（聚类、PCA和fastsimcoal2；表4）支持毛里求斯和留尼汪分别从马达加斯加独立殖民。PCA和FST值显示，两个火山岛之间的遗传分化更大（FST=0.15），相距164公里，而与马达加斯加之间的距离≥665公里（Thébaud等人，2009年；图2）。虽然遗传距离和地理距离之间的缺乏相关性与其他一些研究鸟类殖民多个岛屿的研究（例如，Spurgin等人，2014年）相反，但在印度洋岛屿鸟类中并不罕见（Warren等人，2013年）。不太可能是一个混合的马达加斯加种群影响了我们对独立殖民的推断，因为一个模型认为马达加斯加种群是由毛里求斯和留尼汪种群混合而成的（MOD1）在生物地理学上不太可能（涉及从小型遥远岛屿殖民大型陆地），并且得到的支持较少。相比之下，毛里求斯和留尼汪独立殖民的模型显示出显著更好的AIC分数。在测试N. picturata的殖民历史情景时，这项研究强调了使用综合方法验证人口模型生物学相关性的必要性。尽管共祖框架有助于比较复杂的情景，但它不能保证拟合优度，因此后验验证是必不可少的。我们通过最大化观察到的和预测的SFS之间的相似性来确定了一个可能的模型。包括最近从马达加斯加的迁移（MOD6_MIG_TR）的模型获得了最低的AIC和最高的似然性，为这一情景提供了第一线证据（表4）。有趣的是，这个模型（MOD6_MIG_TR）为两个火山岛的当前种群规模和迁移率提供了两组不同的参数估计。HNLT表明当前的种群规模很大，从马达加斯加的迁移率可以忽略不计。相反，LNHT表明存在显著的瓶颈，强度约为950倍，并与高迁移率相关（表3）。鉴于推断的瓶颈时间（307-370万年后）和迁移事件（270-370万年后），以及历史上关于减少和故意重新引入的报告（Cheke和Hume，2008年；Probst，2002年），瓶颈和迁移很可能是由人类介导的。尽管HNLT获得了更高的似然性（一条支持线；表4），但它与关于引入和种群减少的历史报告不一致。为了确定最生物学相关的情景，我们使用额外的总结统计量（如IICR和LD指标）来评估拟合优度。来自HNLT和LNHT的ML参数值的模拟成功再现了毛里求斯和留尼汪观察到的遗传结构和多样性（图S16-S19）。相比之下，随时间变化的IICR模拟显示了与瓶颈一致的最新人口轨迹。从观察数据中推断IICR可能无法揭示这样的变化，因为顺序马尔可夫方法可能无法检测到共祖率的近期变化，特别是当遗传多样性很大时（Palsb?ll等人，2013年）。确实，PSMC和SMC++对LNHT模拟数据进行分析时未能捕捉到最近的瓶颈现象（见图S25）。此外，基于SFS的统计指标（如Tajima的D值和遗传多样性）以及IICR方法也无法有效区分HNLT和LNHT，因为SFS统计指标的表现相似，且IICR的推断对最近的人口变化（过去约20,000年内的变化；Li和Durbin 2011；Mather等人2020）缺乏敏感性。因此，我们还研究了LD衰减的模式，这种方法有可能揭示更近期的事件，因为重组发生的频率高于突变（Boitard等人2016）。通过分析理论和观察到的LD衰减模式，我们发现这些模式可以区分HNLT和LNHT，尽管两种方法都无法完全解释所有观察到的数据，但在较高的重组距离下，LNHT能更好地预测观察到的LD衰减情况，这反映了近期的人口变化（见表4）。我们注意到，这一结果一方面支持LNHT，另一方面也支持HNLT，即使LNHT与历史记录更为吻合（Cheke和Hume 2008；Probst 2002）。我们的研究表明，仅依赖遗传多样性指标进行物种保护可能会产生误导。事实上，毛里求斯和留尼汪岛上的遗传多样性对于殖民种群来说是非常异常的（Wittwer等人2023）。这种多样性体现在全基因组的HO值上，其范围为5到7 het/kb（见图5），远高于许多经历了近期种群缩减的物种（例如，长须鲸的HO值为0.5–2 het/kb；Nigenda-Morales等人2023）。尽管如此高的遗传多样性可能表明种群规模较大且稳定，但我们的综合分析表明，在这种情况下，观察到的多样性更可能是由LNHT中识别的近期迁移事件所塑造的，这与人类活动增加了马斯克林群岛种群遗传多样性的观点一致。尽管使用基因型似然值对亚化石的群体杂合度进行解释时应谨慎（见S10），但估计结果与这种人类介导的迁移情景相符（见图S26），其中人类出现之前的N. picturata亚化石（Sub1）的杂合度低于人类出现之后的亚化石（Sub13）。对近期经历瓶颈的种群的研究表明，少数个体的迁入可以显著增加遗传多样性（?kesson等人2016；Consuegra等人2005；Mills和Allendorf 1996）。除了LD衰减模式外，还有其他证据（见表4）支持从马达加斯加进行的近期迁移，这些证据来自基于重组的人口统计分析（PSMC）和基于PLINK的ROH估计，这些分析显示毛里求斯的两个个体（205E和225E）表现出与其种群中其他个体不同的遗传特征（见图6和7），这些特征不能归因于序列深度不足（见表S1）。由于使用不同方法估计ROH得到的结果不同，因此在解释时需要谨慎（见S8和图S13）。尽管如此，使用PLINK分析时，这些个体的SROH和NROH值较低（接近马达加斯加的个体；见图6），这与它们是不同谱系混合后产生的后代一致（Ceballos等人2018；Foote等人2021；Pemberton等人2012）。SNMF聚类分析（K=2，具有最低的交叉熵）显示，这些异常个体在毛里求斯种群中的混合血统比例最高，这为马达加斯加和毛里求斯谱系之间的近期混合提供了第五条证据（见表4）。总体而言，这些结果反映了人口过程（如人类介导的迁移）对遗传多样性形成的重要影响。总之，我们的发现支持将N. picturata重新定义为同时原生于毛里求斯和留尼汪岛的物种。了解过去人类活动对物种分布的影响是全球保护工作中的一个重要问题（Gallardo等人2015）。直到上个十年，确定人类活动影响的一个主要挑战是难以利用遗传数据解析近期的人口历史。我们的研究表明，亚化石数据可以补充和验证全基因组模型在近期时间尺度上的人口推断。我们发现这两种独立的推断方法得出了一致的结果，进一步增强了关于人类殖民前的证据。我们的研究还强调了使用独立于推断过程的总结统计量进行严格的事后验证的必要性，以减少接受生物学上不相关模型的可能性。未来，保护遗传学研究将从将多种总结统计量整合到统一分析中受益。

作者贡献：
B.H.W.提出了最初的项目构想。B.H.W.、C.T.、S.M.、G.A.、L.O.、J.P.H.和L.C.获得了资金支持。N.D.和B.H.W.撰写了手稿，S.P.H.、C.T.、R.L.-J.和G.A.也参与了手稿的撰写和编辑工作。N.D.在S.M.和R.L.-J.的帮助下完成了所有统计和群体基因组分析工作。N.D.和A.J.F.V.提取了大部分现代样本。S.P.H.、C.T.和A.J.F.V.负责亚化石的实验室工作。S.M.、R.L.-J.和S.P.H.监督了生物信息学分析工作。M.B.负责PacBio测序。D.E.K.和T.B.S.完成了基因组组装。J.P.H.提供了亚化石样本，并就其他收藏中的亚化石提供了建议。所有作者都批准了最终的手稿提交。

致谢：
我们感谢编辑和三位匿名审稿人的建设性评论和建议。我们感谢SEOR提供留尼汪岛的样本并参与该项目。我们感谢毛里求斯国家公园和保护服务局以及毛里求斯野生动物基金会对我们的工作提供的宝贵帮助和便利。我们还要感谢Céline Bon和Angelika Borowiecka在生物信息学流程方面提供的建议。我们感谢巴黎MNHN古生物学收藏馆允许我们对Sub13进行破坏性采样，以及菲尔德自然历史博物馆提供的组织样本。我们感谢Genotoul生物信息学平台（Bioinfo Genotoul；http://bioinfo.genotoul.fr/）和MNHN PCIA集群提供的计算资源。这项研究得到了法国国家研究署（ANR-20-CE02-0009（Suscept-Ext）和DEAL Réunion的支持。高通量PacBio测序工作由居里研究所的ICGex NGS平台完成，该平台得到了法国国家研究署（ANR-10-EQPX-03（Equipex）和ANR-10-INBS-09-08（France Génomique Consortium）项目的资助，以及ITMO-Cancer Aviesan（Plan Cancer III）和SiRIC-Curie项目（SiRIC Grant INCa-DGOS-465和INCa-DGOS-Inserm_12554）的支持。PacBio数据的管理、质量控制和初步分析由居里研究所的生物信息学平台完成。开放获取出版物的资金由COUPERIN CY26提供。

资金支持：
这项工作得到了法国国家研究署（ANR-20-CE02-0009、ANR-10-EQPX-03、ANR-10-INBS-09-08）、留尼汪岛环境、规划和住房部门以及居里研究所（INCa-DGOS-465、INCa-DGOS-Inserm_12554）的支持。

利益冲突：
作者声明没有利益冲突。

数据可用性声明：
全基因组鸟枪法项目（参考基因组读段和组装）已存放在GenBank中，编号为JBJUWA000000000。其他样本的原始序列读段已存放在PRJNA1400786名下。本研究的好处在于我们能够通过上述方式在公共数据库上共享数据。我们与提供遗传样本的岛屿的保护组织建立了合作关系，我们的研究也解决了他们在马达加斯加龟 Dove 保护管理方面的一些关切。