赤道半球屏障影响了贝洛诺斯蝴蝶（Belenois）迁徙行为的多样化

《Molecular Ecology》：An Equatorial Hemispheric Barrier Shapes the Diversification of Migratory Belenois Butterflies

　　摘要

生物地理屏障通常被认为是阻碍或严重限制物种扩散和基因流动的显著地理特征。然而，在连续的适宜栖息地中，由于生态或行为限制，也可能出现交配障碍。迁徙昆虫展现出穿越广阔地理范围的非凡能力，甚至能够跨越山脉、沙漠和海洋等显著的地形特征。然而，它们的移动并非不受限制：它们的活动受到季节性动态的影响，这些季节性动态决定了它们穿越这些地形的可行性。半球性、两个纬度半球中相反的季节性模式以及定向线索，被认为是影响迁徙昆虫多样化的潜在非生物因素。在这里，我们利用群体基因组数据来研究非洲迁徙蝴蝶（Belenois属）的多样化模式。我们发现Belenois aurota在赤道附近存在明显的系统地理分界，而在Belenois creona中，北部和南部非洲种群之间出现了种群结构的分化，这与半球性屏障所导致的迁徙分割一致。这些分化主要发生在最后一次冰盛期之前，当时发生了赤道雨林和稀树草原的收缩-扩张等重大环境变化。这加强了长期非生物因素（如半球性）在限制南北扩散中的作用这一假设。鉴于即使在肯尼亚、乌干达和坦桑尼亚的B. aurota同域种群中也未检测到基因流动，我们认为来自北半球和南半球的种群代表不同的物种，并因此将Belenois syrinx（Wallengren 1860）重新确立为南半球谱系的合法名称。我们的发现提供了昆虫迁徙分化的基因组证据，这些分化在缺乏物理屏障的情况下仍然发生，突显了半球特异性适应在驱动迁徙昆虫的生殖隔离和多样化中的作用。

1 引言

运动在物种多样化中扮演着关键但矛盾的角色。虽然扩散可以通过使物种扩展到其原生范围之外而导致隔离（Gillespie等人，2012年），但它也促进了种群间的基因流动。迁徙昆虫通常在广阔的分布范围内活动，常常跨越不同的生物群落。这些异常大的分布范围可能反映了与迁徙相关的广泛生态耐受性和复杂的性状组合，这表明与非迁徙昆虫相比，有效的基因流动障碍可能较少见，或者需要更强的选择机制（García-Berro等人，2022年）。尽管如此，迁徙昆虫物种面临的生态复杂性可能为生殖隔离的出现提供了机会，即使在没有物理屏障的情况下也是如此，这一点在鸟类和鱼类的少数案例研究中得到了证明（Turbek等人，2018年）。导致生活史差异的一个潜在驱动因素是半球间非生物条件的不对称性，例如环境线索的时间和可预测性（Chown等人，2004年）。事实上，热带地区历史上既是进化的宝库，也是物种形成的活跃场所（Bernardo-Madrid等人，2025年），以不对称的模式将谱系传播到温带地区。这一过程得益于赤道区域的独特特征：它是地球上唯一一个两个气候相似的纬度带直接接触的区域。赤道区域的广阔性不仅作为一个生态走廊发挥作用，对于分布在两个半球的亚热带物种来说，它也可以成为一个重要的生物地理屏障。在地质时间尺度上，这种半球不对称性——尤其是在主要的气候转变期间——塑造了北半球和南半球的独特进化轨迹，影响了灭绝、隔离和扩散动态（Gaston 2000年；Crame 2001年）。在迁徙物种中，这些半球差异因赤道两侧季节性和天体定向线索的逆转而进一步加剧，这可能导致不同的迁徙策略和路线。最近有基于基因组证据的研究提出，这种机制可能是迁徙昆虫多样化的潜在驱动因素，例如在非洲发现的彩斑蝶（Vanessa cardui）中，揭示了与染色体倒位相关的半球性迁徙分界（García-Berro等人，2025年）。然而，由于大多数迁徙昆虫的系统地理结构尚未得到充分描述，因此半球性（定义为两个纬度半球中相反的季节性模式和定向线索的存在）作为迁徙多样化驱动因素的更广泛相关性仍然很大程度上未经验证，需要在更多不同类群中进行研究。在这方面，非洲热带地区提供了一个探索迁徙、生态和多样化之间相互作用的理想系统。非洲在赤道上的广泛地理和栖息地连续性为测试半球性对迁徙昆虫多样化的净效应提供了理想条件。作为半球性的替代因素，赤道热带森林也可能对偏好开阔栖息地的广泛类群起到生态屏障的作用，可能导致类似的南北种群分化。然而，赤道非洲的热带森林主要分布在西部和中部地区，而在东部则较为稀少，那里有一条几乎连续的非森林栖息地走廊连接着两个半球。因此，在雨林屏障假设下预期的种群空间分布和遗传结构模式与半球性驱动的分化不同。非洲昆虫群落的生物区域化仍然是一个未解之谜，这留下了我们对历史上和当前的物理屏障及扩散走廊如何影响大陆上物种多样化和分布的基本理解上的空白。尽管如此，非洲热带地区仍可以大致划分为几个主要生态区：（i）广泛的稀树草原；（ii）低地雨林区，这是蝴蝶和其他昆虫最丰富的栖息地；（iii）山区，主要由东非的孤立山脉组成，表现出高度的局部特有性；（iv）南部非洲生态区，包括Fynbos和Succulent Karoo，虽然面积相对较小，但拥有异常的生物多样性和特有性（Allsopp等人，2016年；Talavera等人，2020年；Pringle和Woodhall，2025年）；以及（v）沙漠，只有少数蝴蝶物种全年存在，形成了对大多数蝴蝶类群的有效屏障（Larsen 1995年）。有趣的是，大多数非洲热带蝴蝶仅限于这些生态区中的一个或几个。实际上，不到1%的物种在所有非洲热带生态区中普遍存在（Larsen 1995年）。迁徙物种属于这一罕见群体，但它们倾向于避开低地雨林，而是与具有明显季节性的生态区相关联。这种模式强调了在非洲热带蝴蝶多样性中，生态特化普遍超过普遍性的趋势。尽管非洲热带蝴蝶的分类学清单已经显著增加——从第一次区域汇编的1600种（Aurivillius 1898年）增加到目前的4000多种（Williams 2020年；Sáfián和Siklósi 2025a）——但与其他生物区域相比，分子方法仍然被严重低估（French等人，2023年）。相比之下，在研究更充分的地区进行的分子调查揭示了显著的隐秘多样性，即使在传统上知名的动物群中也是如此（Dinc?等人，2011年，2015年，2021年；Vod?等人，2015年；Hinojosa等人，2019年；Li和Wiens 2023年）。随着这些努力扩展到热带地区，类似的模式可能会出现（McBride等人，2009年；Basset等人，2025年），因此预计基于系统发育物种概念的分子系统学将继续揭示非洲热带生态区中的隐藏分类多样性。在同域情况下，屏障效应表现为种群间有效基因流动的减少，可以通过多种基因组机制产生，例如染色体重排或遗传不兼容性，并且通常由进化过程（如不同的局部适应）驱动（Ravinet等人，2017年）。然而，在迁徙昆虫中，空间上的同域并不一定等同于时间上的重叠。相反，迁徙动态可能会促进时间上的异域性，即两个不同的迁徙种群在一年中的不同时间利用相同的地理区域，从而有效地减少基因流动。这种称为异时性的过程越来越被认为是生殖隔离和潜在物种形成的强大机制，特别是当物候差异受遗传控制时（Taylor和Friesen 2017年）。虽然异时性在昆虫中已被广泛记录，包括蝴蝶（例如，Eastwood等人，2023年），但其作为迁徙物种分化驱动因素的作用仍然很大程度上未被探索。Belenois属（Pieridae）的蝴蝶代表了一个研究这些问题的理想系统。该属包含28个物种，主要在非洲热带地区多样化（Williams 2020年），在东方和澳大利亚地区的代表性较小（Yata 1985年；Parsons 1998年；Braby 2000年）。其中，三种广泛分布的物种——Belenois aurota、Belenois java和Belenois creona——特别具有高度迁徙性，形成了周期性的种群爆发和迁徙群，经常影响生态系统和社会意识（例如，Taylor 1959年；Broekhuysen和Broekhuysen 1960年；Milne和Hill 1964年；Kenyi 1980年；Henning 1987年；Braby 2005年，2016年；Haldhar等人，2016年；Drame Yaye和Hamidou 2020年；Fitchett等人，2022年；Chipps和Underhill 2023年；John等人，2024年，2025年）。B. aurota的分布范围最广，覆盖非洲和南亚，而B. creona仅限于非洲，B. java主要分布在澳新地区。尽管目前还没有Belenois的完整系统发育树，但时间校准的蝴蝶系统发育估计该支系起源于约2990万年前，B. java在约560万年前从一组非洲物种（B. creona、B. subeida、B. calypso、B. raffrayi）中分化出来，但没有包括B. aurota（Kawahara等人，2023年；Zhang等人，2023年）。B. java的离散分布及其与非洲亲属的分化与最近的远距离扩散事件或范围扩张一致，而不是古老的冈瓦纳分离。在这里，我们研究了两种具有迁徙倾向的广泛分布的Belenois蝴蝶（B. aurota和B. creona）的系统地理结构，以测试半球性是否可以作为高度迁徙昆虫基因流动的屏障。首先，我们比较了从GBIF获取的年度物候模式，以评估异时性是否是这两种物种北部和南部种群之间合理的隔离机制。然后，我们检查了种群结构的基因组特征、半球谱系之间的基因渗入情况，以及距离和地理屏障的隔离效应。通过人口建模，我们测试了获得的模式是否与半球性迁徙分界和/或异时性的预期一致。最后，我们评估了检测到的谱系间遗传分化是否反映在表型变异中。

2 材料与方法

2.1 空间时间分布模式

我们从全球生物多样性信息设施（GBIF，https://www.gbif.org/es/）获取了B. aurota和B. creona的分布数据，以探索相对丰度的纬度趋势。在过滤掉没有月份分配的观察数据后，我们获得了6607条B. aurota的观察记录（GBIF 2025a）和5707条B. creona的观察记录（GBIF 2025b），这些记录的经度范围在?23.6°至59.3°之间，纬度范围在?34.4°至33.8°之间，因此主要涵盖了非洲大陆和中东部分地区（图1）。然后使用这些观察数据计算了已建立的纬度分区的月度相对丰度（图1）。我们还使用1981-2025年的CHIRPS v3时间序列数据集（Funk等人，2015年；CHIRPS3 2025）计算了每个纬度分区的月平均降水量，作为降雨季节性的代理指标。由于我们旨在研究非洲内的迁徙动态，我们还独立获取并分析了B. aurota在亚洲的分布数据（经度范围在65.6°至93°之间，纬度范围在6.1°至40.5°之间）（2099条观察记录，图S1）（GBIF 2025c）。同样，我们使用7789条观察记录（经度范围在105.2至253.7°之间，纬度范围在?41.1°至?6.1°之间）研究了B. aurota在澳新地区的假定姐妹物种B. java的分布模式（图S1）。

2.2 采样

为了测试迁徙屏障是否由半球性引起，并估计系统地理模式，我们从非洲和亚洲的自然分布范围内收集了38个B. aurota样本和21个B. creona样本的分子数据（图1）。此外，我们还选取了七个作为外群的个体，包括B. subeida、B. zochalia和B. calypso（表S1）。这部分数据集的一部分——四个B. aurota个体和四个B. creona个体，加上一个B. gidica个体和两个B. java个体——也被用于系统基因组学研究（表S1）。这些标本保存在巴塞罗那植物研究所（IBB，CSIC，西班牙）和进化生物学研究所（IBE，CSIC-UPF，西班牙）的收藏中。

2.3 RAD测序和组装

2.3.1 DNA提取、文库制备和测序

DNA是从保存在100%乙醇中并在?20°C下储存的蝴蝶胸部中提取的。每个样本大约使用了10毫克的组织，按照Qiagen DNeasy Blood & Tissue Kit的制造商协议进行处理。所有个体的双酶消化限制性位点相关DNA（ddRAD）基因组文库都是根据Peterson等人（2012年）的方法获得的。至少300纳克的基因组材料使用EcoRI-HF和BfaI（New England Biolabs）酶在37°C下消化3小时。在30到150纳克的消化DNA中，使用T4 DNA Ligase（New England Biolabs）添加了独特的5-bp内联条形码。反应在37°C下进行30分钟，然后在65°C下进行10分钟。在整个过程中，使用Agencourt AMPure磁珠进行反应清洗，并使用Qubit dsDNA Assay Kit（Invitrogen）对DNA进行定量。从32到48个具有独特条形码的样本池中选择300 bp（264-336 bp）的片段，使用PippinPrep仪器和2%琼脂糖卡带（Sage Science，美国贝弗利，MA）进行分离。为了增加DNA的数量并添加Illumina P5条形码，使用Phusion High Fidelity Kit（New England Biolabs）进行了最后一次10轮PCR。使用Bioanalyzer仪器（Agilent Technologies）测量片段的大小分布和浓度。最终将文库按等摩尔量合并，并在瑞典斯德哥尔摩的国家基因组基础设施（NGI）上以150 bp的双端读取方式进行测序。样本池与无关文库一起分布在泳道上，并加入了10%的PhiX Control Library以增加多样性。

2.3.2 ddRAD数据处理和组装
来自66个个体的原始读取数据使用ipyrad v0.9.81（Eaton和Overcast 2020）进行解复用和从头组装。大多数参数设置保持默认值，包括序列聚类阈值85%，除了数据类型（datatype = pairddrad）、限制性突出端（restriction_overhang = AATTC, TAG）和更严格的适配器过滤（filter_adapters = 2）。然后过滤掉产生非常少位点的样本（14T865、14T355、18E246、21B363、11K929和14K409）。分别为B. aurota和B. creona生成了用于群体结构分析（无外群）和系统发育推断（有外群）的组装结果。过滤掉在少于20%的个体中出现的位点。样本统计信息见表S2。

2.4 群体结构分析
我们使用ddRAD数据集来探索B. aurota和B. creona在其整个分布范围内的群体结构。我们最初使用ipyrad中实现的ipyrad.pca工具（Eaton和Overcast 2020）进行了主成分分析（PCA），每个位点使用一个随机选择的SNP（RAD片段）来减少连锁效应。为了评估系统地理模式，我们使用了STRUCTURE v.2.3.4（Pritchard等人2000）。对于每个K值，我们在10次重复实验中进行了K值从1到5的分析，使用了一个包含每个RAD片段的一个SNP的数据集（ipyrad的.usnps文件）。总共进行了1,000,000次迭代，前200,000次被丢弃作为燃烧期。STRUCTURE的配置文件可在https://osf.io/ebjgt/获取。最佳K值是使用Evanno等人（2005）的方法估计的。我们还使用fineRADstructure v0.3.2进行了进一步的人口结构分析（Malinsky等人2018），使用在https://github.com/edgardomortiz/fineRADstructure-tools上提供的脚本生成的完整SNP矩阵。首先使用RADpainter paint命令生成共祖矩阵，然后使用finestructure命令将个体分配到群体中，进行100,000次燃烧期和100,000次进一步样本迭代进行MCMC，保留每1000个样本中的一个。根据不同个体被分配到北部和南部群体的情况，使用pixy v.1.2.6（Korunes和Samuk 2021）计算群体内的核苷酸多样性（π）和群体间的绝对核苷酸差异（dxy）。为此，通过从ipyrad .loci文件生成不变位点的VCF文件，并使用bcftools concat命令（Danecek等人2021）将其与ipyrad生成的vcf文件结合，获得了RAD位点的所有位点的VCF文件。最后，使用Bioanalyzer仪器测量了文库池的片段大小分布和浓度。文库最终以等摩尔量合并，并在瑞典斯德哥尔摩的国家基因组基础设施（NGI）上的Illumina NovaSeq6000 S4平台上进行了150 bp双端读取的测序。样本池与无关文库一起分布在泳道上，并加入了10%的PhiX Control Library以增加多样性。

2.3.2 ddRAD数据处理和组装
来自66个个体的原始读取数据使用ipyrad v0.9.81（Eaton和Overcast 2020）进行解复用和从头组装。大多数参数设置保持默认值，包括序列聚类阈值85%，除了数据类型（datatype = pairddrad）、限制性突出端（restriction_overhang = AATTC, TAG）和更严格的适配器过滤（filter_adapters = 2）。然后过滤掉产生非常少位点的样本（14T865、14T355、18E246、21B363、11K929和14K409）。分别为B. aurota和B. creona生成了用于群体结构分析（无外群）和系统发育推断（有外群）的组装结果。过滤掉在少于20%的个体中出现的位点。样本统计信息见表S2。

2.4 群体结构分析
我们利用ddRAD数据集来探究B. aurota和B. creona在其整个分布范围内的群体结构。首先，我们使用ipyrad中实现的ipyrad.pca工具（Eaton和Overcast 2020）进行了主成分分析（PCA），每个位点使用一个随机选择的SNP（RAD片段）来减少连锁效应。为了评估系统地理模式，我们使用了STRUCTURE v.2.3.4（Pritchard等人2000）。分析针对K值从1到5进行了10次重复实验，每个K值都使用了一个包含每个RAD片段的一个SNP的数据集（ipyrad的.usnps文件）。总共进行了1,000,000次迭代，前200,000次被丢弃作为燃烧期。STRUCTURE的配置文件可在https://osf.io/ebjgt/获取。最佳K值是使用Evanno等人（2005）的方法估计的。我们还使用fineRADstructure v0.3.2进行了进一步的人口结构分析（Malinsky等人2018），使用在https://github.com/edgardomortiz/fineRADstructure-tools上提供的脚本生成的完整SNP矩阵。首先，使用RADpainter paint命令生成了共祖矩阵，然后使用finestructure命令将个体分配到群体中，进行100,000次燃烧期和100,000次进一步样本迭代进行MCMC，保留每1000个样本中的一个。根据不同个体被分配到北部和南部群体的情况，使用pixy v.1.2.6（Korunes和Samuk 2021）计算了群体内的核苷酸多样性（π）和群体间的绝对核苷酸差异（dxy）。为此，我们从ipyrad的.loci文件生成了RAD位点的所有位点的VCF文件，并使用bcftools concat命令（Danecek等人2021）将其与ipyrad生成的vcf文件结合。最后，使用pixy工具运行，设置窗口大小为1000，有效地计算了整个获得的RAD片段的π和dxy。IQ-TREE2 v.2.0（Minh等人2020）被用来推断两个物种的串联比对的最大似然（ML）树，使用了最佳拟合替换模型（-m TEST）和1000次超快自举复制来评估分支支持。

2.5 全基因组测序和系统基因组推断
为了研究所研究的Belenois物种及其地理结构化群体之间的系统发育关系，我们为11个样本生成了全基因组测序数据。这个数据集包括来自B. aurota和B. creona北部和南部群体的各两个个体，一个来自广泛分布的非洲B. gidica的代表性个体，以及两个来自澳大利亚物种B. java的个体，B. java可能是B. aurota的姐妹群或近亲（表S1）。使用Illumina TruSeq DNA PCR-Free协议准备了2×150 bp的双端文库，并在NovaSeq 6000 S4平台上进行了测序。使用OMA浏览器https://omabrowser.org/oma/export_markers（Altenhoff等人2021）为这11个测序样本检索了3892个同源标记，这些标记基于七个注释良好的鳞翅目物种的参考基因组：Bombyx mori、Helicoverpa armigera、Danaus plexippus、Heliconius melpomene、Melitaea cinxia和Papilio machaon。设置了至少0.8的最小覆盖比例和无限数量的标记。然后使用Read2Tree v0.1.5（Dylus等人2024）为这些同源标记生成核苷酸序列比对。接下来，使用Gblocks v.091b（Talavera和Castresana 2007）对每个标记应用，以去除潜在的模糊比对区域，使用默认（严格）参数。接下来，按照两种方法进行了系统发育分析。首先，我们使用IQ-TREE2 v.2.0（Minh等人2020）和最佳拟合替换模型（-m TEST）对包含2,165,769个核苷酸位置的串联比对进行了最大似然系统发育推断。使用1000次超快自举复制来评估分支支持。其次，使用ASTRAL-III（Zhang等人2018）推断了物种树，其中在IQ-TREE2 v2.0中为每个同源基因生成了基因树，并计算了局部后验概率（LPP）（Sayyari和Mirarab 2016）作为支持度量。

2.6 基于模型的群体历史情景分析（使用fastsimcoal2）
我们使用基于模型的方法，根据ddRAD数据重建了B. aurota和B. creona北部和南部群体的群体历史，特别关注估计分化时间、有效群体大小（Ne）和检测群体间的基因流。作为群体建模的输入，我们首先为每个物种生成了北部和南部群体之间的2D位点频率谱。ddRAD测序的两个特点使得重建SFS特别具有挑战性，即高比例的缺失数据和位于同一RAD位点的SNP子集的紧密连锁。因此，我们使用了easySFS（https://github.com/isaacovercast/easySFS；Gutenkunst等人2009）来最小化潜在的偏差。首先，我们预测了每个群体的样本大小，以优化个体代表性和位点恢复之间的权衡（–proj参数）。其次，选择了最大化标记数量而不引入缺失数据的预测样本大小来重建SFS。对于B. aurota，最佳预测样本大小为北部群体n=8，南部群体n=6；对于B. creona，两个群体都预测为n=6个个体。还选择了每个RAD位点最多一个SNP，以减少我们数据集中的标记连锁。根据不同个体被分配到北部和南部群体的情况，使用pixy v.1.2.6（Korunes和Samuk 2021）计算了群体内的核苷酸多样性（π）和群体间的绝对核苷酸差异（dxy）。为此，我们从ipyrad的.loci文件生成了所有位点的VCF文件，并使用bcftools concat命令（Danecek等人2021）将其与ipyrad生成的vcf文件结合。最后，使用–window_size 1000运行pixy，有效地计算了整个获得的RAD片段的π和dxy。IQ-TREE2 v.2.0（Minh等人2020）被用来推断两个物种的串联比对的最大似然（ML）树，使用了最佳拟合替换模型（-m TEST）和1000次超快自举复制来评估分支支持。

2.5 全基因组测序和系统基因组推断
为了研究所研究的Belenois物种及其地理结构化群体之间的系统发育关系，我们为11个样本生成了全基因组测序数据。这个数据集包括来自B. aurota和B. creona北部和南部群体的各两个个体，一个来自广泛分布的非洲B. gidica的代表性个体，以及两个来自澳大利亚物种B. java的个体，B. java可能是B. aurota的姐妹群或近亲（表S1）。使用Illumina TruSeq DNA PCR-Free协议准备了2×150 bp的双端文库，并在NovaSeq 6000 S4平台上进行了测序。使用OMA浏览器https://omabrowser.org/oma/exportmarkers（Altenhoff等人2021）为这11个测序样本检索了3892个同源标记，这些标记基于七个注释良好的鳞翅目物种的参考基因组：Bombyx mori、Helicoverpa armigera、Danaus plexippus、Heliconius melpomene、Melitaea cinxia和Papilio machaon。设置了至少0.8的最小覆盖比例和无限数量的标记。然后使用Read2Tree v0.1.5（Dylus等人2024）为这些同源标记生成了核苷酸序列比对。接下来，使用Gblocks v.091b（Talavera和Castresana 2007）对每个标记应用，以去除潜在的模糊比对区域，使用默认（严格）参数。接下来，按照两种方法进行了系统发育分析。首先，我们使用IQ-TREE2 v.2.0（Minh等人2020）和最佳拟合替换模型（-m TEST）对包含2,165,769个核苷酸位置的串联比对进行了最大似然系统发育推断。使用1000次超快自举复制来评估分支支持。其次，使用ASTRAL-III（Zhang等人2018）推断了物种树，其中在IQ-TREE2 v2.0中为每个同源基因生成了基因树，并计算了局部后验概率（LPP）（Sayyari和Mirarab 2016）作为支持度量。

2.6 基于模型的群体历史情景分析（使用fastsimcoal2）
我们使用基于模型的方法，根据ddRAD数据重建了B. aurota和B. creona北部和南部群体的群体历史，特别关注估计分化时间、有效群体大小（Ne）和检测群体间的基因流。作为群体建模的输入，我们首先为每个物种生成了北部和南部群体之间的2D位点频率谱。ddRAD测序的两个特点使得重建SFS特别具有挑战性，即高比例的缺失数据和位于同一RAD位点的SNP子集的紧密连锁。因此，我们使用了easySFS（https://github.com/isaacovercast/easySFS；Gutenkunst等人2009）来最小化潜在的偏差。首先，我们预测了每个群体的样本大小，以优化个体代表性和位点恢复之间的权衡（–proj参数）。其次，选择了最大化标记数量而不引入缺失数据的预测样本大小来重建SFS。对于B. aurota，最佳预测样本大小为北部群体n=8，南部群体n=6；对于B. creona，两个群体都预测为n=6个个体。还选择了每个RAD位点最多一个SNP，以减少我们数据集中的标记连锁。根据这一策略，最终的2D SFS包括B. aurota的8606个SNP和B. creona的6625个SNP。单倍性位点按比例减少，以考虑保留的分离SNP的数量并保持多态性位点的比例，这对于推断群体事件的时间特别有信息量（Kautt等人2016）。对于每个物种，我们设计了七个替代的群体历史情景，在fastsimcoal v2.8（Excoffier等人2021）中实现。前三个模型假设一个祖先群体分裂成两个群体，两者都没有经历瓶颈。模型M1假设分裂后完全的生殖隔离；M2包括群体间的持续迁移；M3模拟了每个方向的单一迁移事件，作为不频繁但反复发生的基因流的代理。其余四个模型涉及一个创始人事件，其中一小部分来自祖先群体的个体经历了瓶颈，随后有效群体大小（Ne）逐渐恢复。在模型M4和M5中，北部群体是经历瓶颈的新建立的群体，而在M6和M7中，这种情况适用于南部群体。模型M4和M5假设分裂后完全的生殖隔离，而模型M5和M7包括群体间的持续迁移。包含瓶颈的模型为物种的群体历史引入了方向性，因为瓶颈群体代表新建立的群体，而另一个群体代表祖先群体。所有模型的图形表示以及特定模型的参数在图S2中提供，实现这些模型的代码可在OSF仓库https://osf.io/ebjgt/获取。使用easySFS获得的2D-SFS作为模型拟合的输入。所有参数值的先验遵循均匀分布，除了迁移率，它遵循对数均匀先验分布，以更好地捕捉通常较低的迁移比例并允许在几个数量级上的变化。重要的是，参数估计不受先验限制的严格约束，允许模型探索超出这些范围的值，当这提高了运行的似然性时。突变率设置为每位点每代2.9e-9（Keightley等人2015），这是迄今为止唯一可用的蝴蝶的实证估计。每个模型独立运行了100次，每次运行n=100,000次共祖模拟和L=50次ECM循环。对于每个模型，保留了最大估计似然（MaxEstLhood）的重复实验进行模型比较。模型选择使用Akaike信息准则（AIC）进行。模型选择后，使用fastsimcoal v2.8（Excoffier等人2021）根据最佳运行的参数估计值模拟了50个独立的SFS数据集。然后使用每个模拟的SFS作为输入，重新运行模型，n=50,000次模拟和L=30次ECM循环。使用这50次独立运行的参数估计值计算了每个参数的均值和95%置信区间。

2.7 使用F统计量的渐渗测试
为了检测两个物种北部和南部群体之间的潜在基因流，我们使用Dtrios within Dsuite v0.5 r58（Malinsky等人2021）在我们的ddRAD数据集上估计了ABBA和BABA位点的数量（D统计量），使用默认设置。该分析通过定义三个组来计算统计量来进行，这三个组之间计算统计量。这种方法不需要预先知道测试群体之间的系统发育关系，唯一的要求是指定一个外群作为参考点。Dtrios程序测试了所有可能的群体三元组组合，并计算了包括D统计量、其p值（表示统计显著性）和f4比率统计量在内的全基因组统计量，后者提供了渐渗祖先比例的估计。在Dsuite中，D统计量始终报告为正数，因为程序自动排序P1和P2，使得D≥0。这确保了D的大小反映了P3和P2相对于P1的额外等位基因共享。我们数据集中的所有个体根据STRUCTURE分析的结果被分为三组：北部和南部，以及第三组包括来自分布两极的个体，这些个体与另一组的杂交概率最低，这是根据它们与赤道的地理距离来确定的。这样做是为了测试沿赤道接触带可能的渐渗情况，因为靠近赤道的样本最有可能发生渐渗。选择了四个来自巴基斯坦的B. aurota个体，以及一个来自南非的B. creona个体（在三个可用个体中）。我们分别使用B. calypso和B. zochalia作为B. aurota和B. creona数据集的外群。

2.8 线粒体DNA数据集
从公共仓库收集了线粒体细胞色素氧化酶I（COI）基因的标准658 bp条形码区域的序列（表S3），并使用Geneious Prime（https://www.geneious.com）进行了比对。B. aurota的公开可用数据包括来自两个半球的10个国家的145个序列。相比之下，B. creona的公共COI序列仅限于非洲南半球的国家，因此没有使用这个数据集进行后续分析。B. aurota的线粒体比对结果用于评估遗传变异的空间模式和评估线粒体-核一致性程度。使用RAxML软件（图S3），基于GTR替代模型和100次自助法复制，推断出了B. aurota的最大似然系统发育树。

2.9 遗传对比

为了测试距离隔离和主要地理障碍（山脉和赤道地区）对遗传分化的效应，我们采用了GLQ R包中实现的recluster.landscape函数进行遗传景观分析。该方法将样本之间的遗传差异矩阵与其地理坐标结合起来，构建Delaunay三角剖分。使用EEMS仓库中的bed2diffs函数（Petkova等人，2016年）计算了两个ddRAD数据集（B. aurota和B. creona）的成对遗传差异矩阵。首先使用ipyrad流程生成的VCF文件过滤掉外群样本，然后使用VCFtools v0.1.16（Danecek等人，2011年）将其转换为plink格式，接着使用PLINK v1.90（Chang等人，2015年；https://www.cog-genomics.org/plink/1.9/）进行处理。得到的文件被用作bed2diffs的输入。我们使用了版本1（bed2diff_v1），它计算每对个体在具有已调用基因型的SNP上的平均遗传差异。这种方法排除了成对比较中的缺失基因型，并且还生成了一个报告每对个体之间共享SNP数量的矩阵。对于B. aurota的mtDNA数据集，我们使用ape和seqinr包在R（版本4.3.1）中计算了遗传差异矩阵。FASTA比对结果被导入R中，转换为DNAbin格式，并使用dist.dna函数计算成对距离。遗传距离是使用‘raw’ p-distance模型估计的，该模型在计算过程中会删除缺失数据。对于三角剖分的每条边，我们假设陆地上的扩散成本低于海洋障碍，计算了最小成本路径。为此，我们从WorldClim年平均温度层（2.5弧分分辨率；https://www.worldclim.org）中获取了研究区域的陆地-海洋栅格图，将单元格分类为陆地或海洋。这个栅格图被转换为渗透性表面，通过为陆地单元格分配一系列渗透性值（0.5, 1, 2, 3, 5, 8, 12），同时将海洋单元格的值固定为1，从而表示陆地与海洋渗透性的不同比例（Luo等人，2022年）。对于每个渗透性比例，使用gdistance R包的transition和geoCorrection函数生成了一个过渡层，并使用shortpath函数计算了Delaunay三角剖分中连接的样本之间的最小成本路径。从GEBCO（https://www.gebco.net）获取了一个深度-海拔层（2.5弧分分辨率），并将其裁剪到研究区域内。对于每条最小成本路径，我们提取了：（i）路径的地理中点（纬度和经度），（ii）路径的总长度，（iii）路径上遇到的最大海拔高度，以及（iv）连接样本之间的遗传距离。最佳陆地-海洋渗透性比例被确定为使遗传距离和路径长度之间的线性相关性（Pearson's r）最大化的那个比例。使用这个表现最佳的景观模型，通过将每个路径中点与遗传距离和路径长度之间的线性回归残差关联起来，绘制了遗传对比图。这些残差在QGIS 2.18中使用逆距离加权方法插值到研究区域内。此外，我们使用mgcv R包中的gam函数拟合了一个广义加性模型（GAM），以Delaunay三角剖分段作为观测单元。遗传距离被用作响应变量，预测因子包括：（i）路径上的最大海拔高度（Max_Alt，平滑项），（ii）路径中点的纬度和经度（Midpoint_Lat和Midpoint_Long，平滑项），以及（iii）路径长度（Path_Length，平滑项）。通过包含基于路径中点之间距离的空间相关结构，考虑了空间自相关性。由于遗传距离是右偏的且非正态分布的，因此假设了Tweedie分布。模型使用限制最大似然（REML）方法进行拟合。

2.10 翅膀的几何形态测量分析

从两个种群中选取了33只B. aurota和20只B. creona个体的前翅（表S1）进行比较形态测量分析。前翅在飞行过程中主要产生推力，功能上受到更多限制，而后翅主要负责转向，表现出更大的随机变异。因为我们的目标是研究扩散的形态相关性，所以我们关注前翅这一更稳定且功能相关的结构。选择样本时考虑了前翅的完整性，以确保可以可靠地评估细胞周围和翅缘的固定标志点；有4只B. aurota个体因缺少标志点而被排除在外。所有翅膀都在标准条件下使用Canon EOS 250D相机拍摄，然后使用R包StereoMorph中的digitzeImages函数在所有样本的腹面数字化标志点，遵循之前在其他蝴蝶物种中进行的 vein 交点方法（例如，Platania等人，2020年）。使用geomorph R包中的gpagen函数对标志点数据进行了广义Procrustes分析（GPA），以消除位置、尺度和方向的非形状变异，并将样本叠加在一个共同的坐标系统中（Rohlf和Slice，1990年）。然后将对齐的Procrustes坐标投影到一个线性切线空间中，以获得用于后续分析的标志点坐标（Adams等人，2023年）。使用mixOmics R包中的plsda函数进行了部分最小二乘判别分析（PLSDA），以根据样本的起源半球（北部或南部）进行区分。选择PLSDA而不是其他方法（例如，线性判别分析或典型变量分析），是因为它对异方差性和多重共线性不那么敏感，这些在几何形态测量数据集中很常见，并且明确寻求预测因子之间的相关性与事先定义的样本组之间的变异之间的最佳平衡（Mitteroecker和Bookstein，2011年）。由于翅膀形态在不同性别之间往往存在差异，而且在迁徙物种中（例如，携带发育中卵的雌性）性别特定的飞行限制可能特别重要，因此根据性别和半球的组合将样本分为四组。首先使用这种分组进行了初步的PLSDA，以可视化整体分化模式。随后，进行了一系列留一法交叉验证的PLSDA，以测试是否可以根据剩余个体制建的模型将每个样本盲目分配到正确的组中，使用mixOmics包中的perf函数。判别准确性量化为正确分配到其起源半球的雄性和雌性的百分比。

3 结果

3.1 季节性出现模式

在定义的纬度区间内，每月的相对丰度显示了与两种物种的迁徙活动一致的动态年度物候学特征（图1）。在赤道纬度（5°N–5°S），B. aurota和B. creona都显示出与长雨季和短雨季时间相对应的双峰丰度峰值。相比之下，在只有一个雨季的纬度区间，出现了单峰模式，反映了与当地气候制度相关的季节性繁殖活动。在B. aurota中，在南非最南端的纬度区间（25°–35°S），3月（澳大利亚秋季）之后观察到丰度显著下降（图1，左侧面板），种群数量一直保持较低直到10月（春季）。在撒哈拉以北的北部带（5°–20°N），种群数量在旱季下降，这与迁徙-移民模式一致。相比之下，B. creona在整个年度内的相对丰度更加稳定，大多数纬度区间的月度波动较小（图1，右侧面板）。这表明其迁徙行为更受地域限制或不太严格。

3.2 ddRAD组装和变异调用

去多重复制后的总读取数为126,103,056，每个个体的平均深度覆盖约为2,137,339个读取（最小296,768，最大12,176,130个读取；表S2）。B. aurota数据集包含229,411个SNP，其中75.33%的位点缺失，总共3,749,805个组装的核苷酸（包含外群）；或者188,579个SNP，其中74.70%的位点缺失，总共3,602,309个组装的核苷酸（不包含外群）。B. creona的相应值为242,645个SNP，其中69.93%的位点缺失，总共5,004,207个nt（包含外群）；或者180,860个SNP，68.23%的位点缺失，总共4,726,229个nt（不包含外群）。

3.3 种群结构

B. creona的种群内核苷酸多样性（π）为0.0103，B. aurota为0.0120。分别考虑B. aurota的北部和南部种群，这些值分别为0.0115和0.0129。北部和南部种群之间的平均核苷酸差异（dxy）分别为0.0147和0.0179。主成分分析（PCA）、STRUCTURE分析和fineRADstructure共祖分析一致地揭示了B. aurota和B. creona的南北遗传结构模式（图2）。在B. aurota中，个体聚类为两个不同的组，主要沿着PC2分离，解释了总方差的3.9%，而PC1解释了4.9%（图2）。STRUCTURE分析在K=2时显示了北部和南部种群之间的明显遗传分离，这种分离进一步得到了成对遗传相关性热图的支持，该热图显示了种群内的强相似性和区域间的明显遗传不连续性（图3）。在B. creona中，PCA也显示北部和南部种群沿着PC2分离（PC1和PC2分别解释了总方差的7.3%和6.2%），这一模式进一步得到了STRUCTURE和共祖分析的证实，两者都揭示了南部和北部谱系的聚类。然而，在这个物种中，来自埃塞俄比亚的一个个体在STRUCTURE和共祖图中显示出潜在的混合，表明东非存在持续的基因流（图3）。

3.4 人口统计历史

fastsimcoal2中模拟的人口统计情景揭示了B. aurota和B. creona建立半球分界线背后的大致相似的历史。对于这两种物种，最简单的模型（M1）提供了最佳的整体拟合（表S4）。该模型描述了一个祖先种群分裂为两个种群（每个半球一个），没有明显的人口收缩，并且在分裂后有完全的生殖隔离（图4）。在B. creona中，一个包括分化后双向迁移的替代模型（M3）显示出略高的似然性，但总体支持度有限（AIC权重=0.154）（表S4）。尽管这些人口统计历史大致一致，但两个物种在估计的有效种群大小和分化时间上有所不同。在两个半球中，B. aurota推断出的Ne值都高于B. creona（图4b）。在B. aurota中，北部种群推断出的Ne值更大（Ne=1,438,732；95% CI：1,414,346–1,463,118），而南部种群为Ne=1,301,798；95% CI：1,276,336–1,327,259）。相比之下，在B. creona中，北部种群的Ne值略低（Ne=980,396；95% CI：956,661–1,004,131），而南部种群为Ne=1,030,443；95% CI：1,009,049–1,051,837）。尽管每个物种内部的半球差异微妙，但北部和南部种群之间的置信区间不重叠，表明在有效种群大小上存在半球差异。此外，在分化之前，B. aurota的祖先种群规模（Ne）比B. creona的更大。图4（在图形查看器中打开，PowerPoint格式）：

(a) Belenois的系统基因组树。该拓扑结构与基于WGS数据集的ASTRAL-III物种树推断结果一致，并且与使用IQ-TREE2进行的串联分析结果相符。节点标签显示了ASTRAL局部后验概率（LPP；0–1），后跟IQTREE2超快自助法支持率（%）。

(b) 为B. aurota和B. creona推断出的最佳人口模型。一个祖先种群在北半球和南半球之间发生了纬度分裂，分化过程中没有出现种群瓶颈现象。显示了50次独立模拟的SFS的平均参数估计值，95%置信区间用括号表示。有效种群规模（Ne）针对单倍体基因组，分化时间（t）以距今的代数表示。分化时间的估计进一步表明了这两个物种之间的时间差异。B. aurota的半球分裂发生在距今570,701代（95%置信区间：561,698–579,704代），而B. creona的分化时间较晚，为439,514代（95%置信区间：429,161–449,866代；图4b）。假设热带非洲的世代时间为6–8周（例如，Kenyi 1980年的数据），这些估计值表明B. aurota的分化时间大约在95,000–71,000年前，B. creona的分化时间大约在75,000–55,000年前。

3.5 系统发育和系统基因组学推断

所有系统发育分析都一致地显示了B. aurota和B. creona之间存在明显的南北分裂。从串联的ddRAD比对中推断出的ML树完全支持这两个物种之间的南北分裂（图2a,b）。使用ASTRAL III对WGS数据集进行的基因树推断得出了与IQ-TREE2串联分析相同的拓扑结构，所有节点都获得了最大的LPP和自助法支持率（图4）。基于COI数据为B. aurota构建的ML系统发育树也支持南北分裂（图S3）。

3.6 渗透测试

为了测试B. aurota和B. creona地理结构化种群之间的潜在基因流，我们对这两个物种进行了D统计分析（ABBA–BABA测试）（表S5）。在B. aurota中，比较结果显示南北个体之间没有显著的基因流证据（D统计量=0.0164，Z分数=0.461，p=0.645）（表S5）。ABBA和BABA模式的计数几乎相等，f4比率=0.0145也表明没有显著的等位基因共享增加信号。在B. creona的情况下，尽管在埃塞俄比亚发现了一个具有混合血统的个体（图2a），但由于该个体无法被归类为属于任一种群，因此渗透测试没有显示出南北个体之间存在渗透现象（D统计量=0.0427，Z分数=1.374，p=0.169）（表S5）。尽管ABBA模式的计数略高于BABA模式（分别为386.8和355.1），但这一信号在统计上并不显著。f4比率（0.0324）也较低，表明混合程度有限或没有混合。

3.7 遗传对比

初步分析将遗传距离与在不同陆地-海洋渗透率（k）下计算的路径长度相关联，结果在不同标记和物种间略有不同。对于B. creona，核标记的最高相关性出现在k=3时；而对于B. aurota，核标记的最高相关性出现在k=5时，线粒体DNA（mtDNA）的最高相关性出现在k=3时（表S6）。广义加性模型显示遗传距离与路径长度之间存在正相关关系，这与隔离距离效应一致。此外，对于B. creona的RAD数据和B. aurota的COI数据，遗传距离也与沿最低成本路径遇到的最大海拔高度显著相关。在所有物种中，纬度都表现出高度显著的非线性效应（有效自由度[edf] > 1，表S7）。当可视化时，这种效应呈现出钟形模式，最大值集中在赤道纬度附近（图3）。这种空间结构还通过将遗传距离的回归残差与空间距离进行插值得到了进一步的支持，这揭示了三个物种之间的一致的地理模式（图3）。

3.8 缺乏翅膀形状的表型分化

GPA后的地标坐标图显示样本之间的变异较小。在两个物种中，这种小变异的大部分可以通过性别来解释，如第一个PLSDA轴所示（图2c）。变量相关性图的检查显示，大多数地标参与了区分雄性和雌性的第一个组分，而且B. creona的雌性个体也倾向于更大（与质心大小相关）。PLSDA的散点图显示来自两个半球的样本之间几乎没有分化，这一点通过第二个PLSDA组分中涉及的形状变量的稀缺性得到了证实。因此，采用留一法（leave-one-out procedure）来识别属于正确半球的雄性和雌性样本，结果显示B. aurota的正确分配比例为52%，B. creona为60%，非常接近随机预期的50%。

4 讨论

4.1 Belenois aurota和B. creona的差异性迁徙动态

B. creona和B. aurota都被广泛认为是强迁徙性蝴蝶，有大量的观察证据支持这一地位，包括关于大规模群集、一致方向飞行以及典型的迁徙行为季节性人口波动的报告（Williams 1930, 1958; Williams 2020; Chowdhury et al. 2021）。然而，对于许多通常被认为是迁徙性的蝴蝶，尤其是那些栖息在非洲热带地区的蝴蝶，关于它们季节性范围变化的详细实证数据仍然很少。需要强调的是，在昆虫个体层面进行迁徙行为表型分析仍然具有挑战性，需要精细的实验方法。同样，量化种群层面的迁徙结果需要长期的大规模监测计划，而这在南方半球尤其罕见（Talavera et al. 2023）。我们在非洲进行的基于纬度的时空筛查揭示了这两个物种在种群动态上的基线差异。在乌干达，观察到B. creona在六月向南迁徙，这与盛行风向相反，并且与其他迁徙性蝴蝶物种同步（Sempala 1970）。六月也是这两个物种在赤道纬度处数量最多的时期（图1）。然而，B. creona在这个月内的密度分布范围更广（大约-15°到20°），而B. aurota的密度峰值仅限于赤道带（图1）。这种对比表明，B. creona可能遵循不太具体的纬度方向，并进行广泛的移动，可能反映了比B. aurota更受区域限制的、可选的迁徙模式。相比之下，B. aurota的迁徙范围从非洲南部跨越半个地球到亚洲，但明显缺席于西地中海地区——这一模式与其他几种泛热带蝴蝶相同。在赤道纬度（约0°），密度显示出双峰分布（图1）。在乌干达，繁殖活动记录在七月和八月，发生在漫长的雨季之后（Kenyi 1980），十一月到二月可能有第二个高峰，跟随短暂的雨季（Safari Ecology Blog 2012）。这种模式与该地区的双峰降雨模式一致，可能反映了季节性的异时性，即南部谱系在五月到八月占主导地位（此时最南端的密度下降），而北部谱系在十一月到二月变得更加普遍（图1，图S4）。在V. cardui中也观察到了类似的模式，尽管在那种情况下，北部遗传谱系进一步将其分布扩展到欧洲，实际上在二月到八月期间从非洲热带地区消失（Talavera et al. 2023; García-Berro et al. 2025）。尽管北部个体之间缺乏精细的遗传结构，但目前尚不清楚大多数B. aurota个体是否参与每年在非洲赤道和中东之间的周期性迁徙，这将意味着穿越撒哈拉沙漠的迁徙严格限制在东部通道。观察数据显示，在北部热带地区（5°N–20°N），B. aurota的密度较低，主要分布在四月到十月之间，这些地区以草原景观为主。这种季节性缺失可以通过以下几种情况来解释：（i）迁徙连续性，即大多数非洲热带个体在春季迁移到中东，然后在秋季返回；（ii）向赤道地区的南部或高海拔地区撤退；（iii）这两种类型的移动的结合。我们假设二月至三月之间的种群减少促使了向北和向南的移动。在这种情况下，一部分北部种群向南移动到赤道地带，因此可能比纯粹的异时性模型更频繁地遇到南部谱系。这一观点得到了肯尼亚（Laikipia）在五月至六月短暂雨季期间发现的六个北部谱系个体的支持（图S4）。两个B. aurota谱系在赤道地区的某种程度的时间重叠可能发生在两个雨季之后。例如，我们的数据集包括一个来自乌干达的样本（17G874；0.85°S，30.27°E），在十二月采集，它在遗传上与南部谱系聚类，而在肯尼亚同时采集的样本（十一月和十二月，>0.74°S）则与北部谱系聚类（图2）。解决这一接触区域的季节性分离程度需要大量的采样工作，这对于理解两个B. aurota谱系之间显著分化背后的进化机制至关重要。虽然季节性异时性可能有助于这种分化，但我们的当前结果表明，导致基因组不兼容性的选择力量可能比在严格异时性下的漂变作用更强。在赤道地带之外，相对丰度通常呈现单峰模式（图1，图S1）。B. aurota在十二月到二月期间定期出现在南半球，通过南非向北或东北方向迁徙，与南半球雨季的高峰期相吻合（例如，Taylor 1959; Broekhuysen and Broekhuysen 1960; Milne and Hill 1964; Terblanche 2020）。在北半球，特别是在亚洲，B. aurota的移动似乎不那么严格地受纬度限制（图S1），可能受到季风驱动的季节性的影响，这在南亚的其他迁徙性蝴蝶中也有很好的记录（Bhaumik and Kunte 2018）。例如，在印度拉贾斯坦邦（约28°N），对Capparis寄主植物上的B. aurota幼虫进行了为期两年的监测研究，发现幼虫数量在十二月达到高峰，而在五月到八月完全消失（Safari Ecology Blog 2012）。B. aurota的一个显著生态特征是它依赖于两个半球内的干旱和半干旱落叶灌木丛——特别是在赤道地带之外——在那里它利用木本常绿Capparaceae植物作为幼虫寄主植物（Haldhar et al. 2016; Terblanche 2020; Dramé Yaye and Hamidou 2020）。这种依赖性表明，除了幼虫寄主植物的可用性之外，温度或季节性寄主使用竞争等因素可能是触发迁徙的关键因素。例如，在卡拉哈里地区——被认为是一个关键的繁殖区域——B. aurota种群在四月到十月期间繁殖（Terblanche 2020），但在冬季几乎消失。这种季节性缺失可以通过以下几种情况来解释：（i）迁徙连续性，即大多数非洲热带个体在春季迁移到中东，然后在秋季返回；（ii）向赤道地区的南部或高海拔地区撤退；（iii）这两种移动类型的结合。我们假设二月至三月之间的种群减少促使了向北和向南的移动。在这种情况下，一部分北部种群向南移动到赤道地带，因此可能比纯粹的异时性模型更频繁地遇到南部谱系。这一观点得到了在肯尼亚（Laikipia）的mtDNA数据集中发现的六个北部谱系个体的支持，这些个体出现在五月至六月短暂雨季期间（图S4）。两个B. aurota谱系在赤道地区的某种程度的时间重叠可能发生在两个雨季之后。例如，我们的数据集包括一个来自乌干达的样本（17G874；0.85°S，30.27°E），在十二月采集，它在遗传上与南部谱系聚类，而在肯尼亚同时采集的样本（十一月和十二月，>0.74°S）则与北部谱系聚类（图2）。解决这一接触区域的季节性分离程度需要大量的采样工作，这对于理解两个B. aurota谱系之间显著分化背后的进化机制至关重要。虽然季节性异时性可能有助于这种分化，但我们的当前结果表明，导致基因组不兼容性的选择力量可能比在严格异时性下的漂变作用更为重要。在赤道地带之外，相对丰度通常呈现单峰模式（图1，图S1）。B. aurota在十二月到二月期间定期出现在南半球，通过南非向北或东北方向迁徙，与南半球雨季的高峰期相吻合（例如，Taylor 1959; Broekhuysen and Broekhuysen 1960; Milne and Hill 1964; Terblanche 2020）。在北半球，特别是在亚洲，B. aurota的移动似乎不那么严格地受纬度限制（图S1），可能受到季风驱动的季节性的影响，这在南亚的其他迁徙性蝴蝶中也有很好的记录（Bhaumik and Kunte 2018）。例如，在印度拉贾斯坦邦（约28°N），对Capparis寄主植物上的B. aurota幼虫进行了为期两年的监测研究，发现幼虫数量在十二月达到高峰，而在五月到八月完全消失（Haldhar et al. 2016）。B. aurota的一个显著生态特征是它依赖于两个半球内的干旱和半干旱落叶灌木丛——特别是在赤道地带之外——在那里它利用木质常绿Capparaceae植物作为幼虫寄主植物（Haldhar et al. 2016; Terblanche 2020; Dramé Yaye and Hamidou 2020）。这种依赖性表明，除了幼虫寄主植物的可用性之外，温度或季节性寄主使用竞争等因素可能是触发迁徙的关键因素。例如，在卡拉哈里地区——被认为是一个关键的繁殖区域——B. aurota种群从春季到初秋繁殖（大约十月到三月）（Terblanche 2020），但在冬季几乎消失，这可能是由于温度较低。

4.2 B. aurota的半球性迁徙分裂

B. aurota的种群结构在东非赤道附近显示出明显的遗传分裂（图2），目前没有明显的地理或生态障碍。两个遗传谱系在肯尼亚、坦桑尼亚和乌干达之间的相对狭窄的纬度带内接触。根据mtDNA数据集（表S3，图3），南部谱系的最北代表个体在乌干达被发现（17G874；0.85°S，30.27°E），而北部谱系的最南代表个体在坦桑尼亚被发现（CBG-A00313-B06；2.9°S，35.35°E），表明在东非至少有2°的纬度范围（约220公里）的重叠区域。在两个季节中，除了南部海岸外，没有在肯尼亚发现南部谱系的样本。坦桑尼亚的大多数样本以及一个来自乌干达的样本与南部谱系聚类，而肯尼亚的大多数样本——除了南部海岸的样本——属于北部谱系（图2和图3）。坦桑尼亚北部的三个地点——Ngorongoro、Serengeti和Busega——在mtDNA数据集中同时包含了两个谱系的个体，确定了这一广阔的草原区域作为潜在的接触带（图3）。在坦桑尼亚总共28个样本中，有3个（约11%）是北部起源的样本，在十月到三月期间与南部谱系的优势个体同时被采集（表S3；图S4）。这表明两个谱系在赤道地带确实是共存的，尽管两个谱系的相对丰度模式可能随季节变化，因此这些地点可能存在部分时间分离。系统基因组学分析进一步证实了B. aurota的北部和南部种群之间的深刻分化。人口建模表明，这种分裂大约发生在距今570,000代之前（BP）。假设每年有6-8代（大约9.5-7.1万年前），这种分化主要发生在最后一次冰盛期（LGM）之前，当时在上一个间冰期结束后，非洲的干旱程度逐渐增加（大约11.5万年前；Stokes等人1997年；Owen等人2018年）。那时，平均气温比现在低约2°C（Chevalier等人2021年）。有趣的是，估计的分化时间与这一更广泛的干旱趋势中的几个短暂湿润期重叠（Lézine等人2017年；Owen等人2018年）。这种巧合挑战了历史气候波动可能导致B. aurota的半球种群之间出现生殖隔离的可能性，因为在它们分化的大部分时间里，半干旱的落叶灌木丛在赤道地区广泛分布。此外，在这段时间内没有明显的物理障碍出现，除了今天仍然存在的中非雨林。由于推断4万年前古环境条件的固有难度（Van Andel和Tzedakis 1996年），以及没有明显的地理或气候障碍阻碍扩散和基因流动，B. aurota的生殖隔离可能是由于生物和/或遗传因素造成的。需要进一步的全基因组分析来解开这种分化的机制。我们观察到B. aurota的Ne估计值与直接测量值之间存在一些差异，因为推断的Ne在北半球较高（图4b），而π在南半球较高。这与中性数据的经典预期不符（Kimura 1971年；Watterson 1975年），尽管这可能由几个技术因素解释。首先，由于ddRAD组装中缺失数据的比例很高，SFS的尺寸必须大幅缩小。这样的样本量较小，不如Excoffier等人（2013年）通常用于基准测试的样本量，可能会增加人口模型拟合的不确定性，并导致运行之间的变异性增加。此外，这里应用的人口模型相当简单，主要关注估计分化时间、种群规模和潜在的基因流动，而不是复杂的人口情景。因此，它们可能无法完全捕捉到两个种群自分化以来所经历的历史事件的复杂性，这可能会影响参数估计，从而影响测试情景的结果。因此，应谨慎解释B. aurota不同种群之间的Ne差异。尽管不完全可比，但对Belenois的直接遗传变异估计似乎落在蝴蝶报告的平均范围内（Mackintosh等人2019年）。Belenois java是一种形态相似且高度迁徙的物种，仅限于澳新大陆和东南亚（印度尼西亚），被认为是B. aurota的近亲（图4）。尽管缺乏关于B. java季节性迁徙的直接追踪数据，但出现记录表明其迁徙模式具有纬度结构（图S1；Braby 2000年）：温带和南亚热带地区的种群通常在澳大利亚春季向南迁移，而北亚热带和热带地区的种群在同一时间向北迁移（Braby 2016年）。然而，虽然方向性似乎因地区而异，但尚未记录到返回运动（Dingle等人1999年）。B. java完全局限于南半球的事实加强了这一观点，即该组的迁徙行为可能是半球特异性的，并且在进化上是保守的。实际上，B. aurota（在亚洲）和B. java之间的分化可能是由半球限制塑造的。尽管南亚和东南亚与澳新大陆之间的地理环境复杂——以马来西亚和印度尼西亚的群岛为主——但这些破碎的景观不太可能对强迁徙物种构成显著的扩散障碍。B. aurota和B. java的祖先可能曾经具有跨半球的分布。实际上，对B. java的出现数据进行纬度筛选显示，1月份在印度尼西亚有高峰记录，随后几个月急剧减少（图S1）。这种模式为考虑印度尼西亚和澳大利亚之间的一定程度的连通性提供了可能性。

4.3 B. creona的遗传结构证据

种群结构和系统基因组分析也证实了B. creona的北部和南部种群之间的分化。B. creona的北部和南部支系之间的平均核苷酸差异低于B. aurota（0.0147对比0.0179），尽管核苷酸多样性有中等差异（0.0103对比0.0120）。这些对B. creona的遗传多样性的直接估计大致反映了两个种群推断的Ne值（图4b）。然而，B. creona的种群分化似乎比B. aurota更近期，发生在大约45万代之前（假设每年6-8代的情况下相当于7.5-5.6万年前；图4b）。与B. aurota不同，B. creona在东非的赤道地区显示出部分混合的迹象。这由来自埃塞俄比亚的一个标本表明，其基因组祖先中南部谱系的比例较大（图2b）。尽管在该地区的B. creona采样仍然有限，但在肯尼亚的赤道线以北发现了南部谱系的个体，这在B. aurota中没有观察到。这表明B. creona的接触区可能更靠北，需要更全面的采样来解决。这一历史基因流动的信号也与人口模型M3一致——根据模型选择标准，M3是第二好的模型（表S4）——它与最佳拟合模型M1具有相同的基本结构，但包括了一个双向迁移事件作为分化后偶发接触的代理（图4b）。对M3的支持表明，M1中推断的严格生殖隔离可能过于简化，B. creona可能经历了一段较长的接触期。虽然需要更广泛的采样来评估混合的程度和结构，但目前的证据指向B. creona的纬度梯度，而不是像B. aurota那样明显的赤道分界。这种检测到的B. creona的种群结构模式通常与两个广泛分布的亚种相符：北部的B. creona creona和南部的B. creona severina（Sáfián和Siklósi 2025b, 2025c, 2025d）。与B. aurota一样，分化时间与撒哈拉以南非洲间冰期后干旱趋势中的间歇性湿润期相吻合（Lézine等人2017年；Owen等人2018年；Nazari等人2025年），尽管具体的环境驱动因素尚不清楚。

4.4 形态和遗传分化是脱节的

B. aurota和B. creona在遗传谱系之间的翅膀形态上都没有明显的分化（图2c），表明尽管存在遗传分化，形态特征仍然保持一致。然而，在B. creona中，遗传分裂——尽管比B. aurota更近期（图4）——大致与传统亚种分类一致，后者将B. creona creona识别为南部物种，B. creona severina识别为北部物种（Williams 2020）。这种分类区分与翅膀颜色的微妙差异有关，B. creona creona被描述为颜色较浅。然而，这一特征的可靠性值得怀疑，因为变化是渐进的且高度主观的，通常无法提供明确的诊断价值（Van Son 1949）。定向选择可以影响迁徙蝴蝶的飞行形态（Dockx 2007；Chan等人2023），但如果在两个半球中高效的飞行性能都很重要，翅膀形态可能受到稳定选择的影响，这与迁徙和定居帝王蝶的选择压力不同（Dockx 2007）。

4.5 迁徙物种中的隐秘物种形成

在迁徙类群的多样化中，一个核心问题是姐妹物种倾向于隐秘地分化还是明显地分化。鉴于它们在长距离扩散和克服主要生物地理障碍方面的显著能力，可以认为异域物种形成可能比同域物种形成更为常见。迁徙可能导致独特的扩散事件，如果扩散的种群在新范围内成功建立，就可以在异域条件下实现隔离（Palahí等人2025）。然而，在连续的景观中，周期性迁徙被认为会抑制稳定空间梯度的形成，除非有强大的基因组或非生物障碍减少基因流动（Sanchez-Donoso等人2022；Lundberg等人2023；García-Berro等人2025）。半球障碍是否能够一致地驱动迁徙蝴蝶的多样化，或者在整个迁徙动物中也是如此，仍然是一个开放且引人注目的问题，这需要跨类群的广泛比较研究。B. aurota的案例为之前在V. cardui中报告的案例增添了另一个例子，其中遗传分化与半球分布和反向迁徙模式一致（García-Berro等人2025）。在V. cardui中，这种分化与染色体倒位有关，这与不同的迁徙策略相关，但其他迁徙物种对也可能反映类似的过程。V. cardui的分布进一步支持了这一点：除了南部非洲外，该物种在其他南半球大陆上几乎不存在。尽管其祖先曾殖民澳大利亚，但物种分化导致了其姐妹物种V. kershawi的出现。在南美洲，尽管多次尝试殖民（Suchan等人2024），V. cardui尚未建立定居种群——这表明半球限制可能影响了这些结果。实际上，最近在澳大利亚建立的帝王蝶观察到了季节性线索的异常解释（Tenger-Trolander等人2023）。尽管归因于迁徙行为的丧失，但半球性可能发挥了作用。这些过程可能不仅限于昆虫。每个半球的独立迁徙进化也可能影响脊椎动物的分化，包括迁徙鸟类（Dingle 2008；Areta等人2021；Helm和Muheim 2021）和海洋哺乳动物如鲸鱼（Jackson等人2014；Pérez-Alvarez等人2021）。B. aurota、B. creona和B. java的案例为探索基因组分化如何映射到昆虫的迁徙行为提供了有希望的机会。了解特定基因组区域（例如，倒位、调控元件）是否与季节性时间、方向或范围使用相关，可以为研究迁徙的遗传结构提供窗口（Liedvogel等人2011；Liedvogel和Lundberg 2014；Merlin等人2019）。然而，半球分化并不是普遍的结果：例如，基于线粒体DNA的对Lampides boeticus的研究（Lohman等人2008）在其广泛的跨半球范围内没有发现种群结构，而在高度分散的蛾类如Utetheisia pulchella中也观察到了类似的混合现象（Zhang等人2024）。

4.6 非热带地区的隐秘物种形成

似乎有理由断言，非洲热带地区的昆虫遗传多样性在很大程度上尚未得到充分研究。只有少数蝴蝶物种具有明确的系统地理框架，如Danaus chrysippus（Smith和Owen 1997；De-Kayne等人2025），或V. cardui（Suchan等人2024；Reich等人2025；García-Berro等人2025），而其他物种如Hypolimnas misippus（Orteu等人2024）或Byciclus anynana（Brakefield等人2009；Nowell等人2017）已成为生态学、遗传学和发育研究的模型物种。Pseudopontia spp.的案例代表了非洲隐秘蝴蝶物种形成的一个例子，其中可能在热带地区出现了五个谱系，有些处于同域状态（Mitter等人2011）。同样，Cymothoe egesta和C. confusa这对隐秘物种在喀麦隆的一个狭窄区域内被发现处于二次接触状态（McBride等人2009）。然而，这些例子是大多数非洲热带昆虫缺乏大规模种群遗传数据的例外。在调查良好的温带地区（如欧洲和北美）进行的比较评估通过系统地理筛选揭示了大量的隐藏多样性，隐秘谱系表明真正的物种丰富度可能被低估了多达10%（Dinc?等人2021；D'Ercole等人2021）。最近的全球评估进一步表明，热带地区的昆虫隐秘多样性水平可能更高（Li和Wiens 2023）。鉴于非洲热带地区拥有异常的昆虫多样性和栖息地完整性的加速威胁，推进我们对分类边界和系统地理结构的理解对于昆虫分类学和保护生物学来说是一个紧迫的优先事项。

5 分类学

Belenois syrinx（Wallengren 1860），重新确立为法定名称。Pinacopteryx syrinx Wallengren 1860（34）。Pieris syrinx（Wallengren）。—Kirby（1871, 457）。Belenois syrinx（Wallengren）。—Butler（1872, 56）[属于‘Cronea组’]。Pinacopteryx gidicae var. syrinx（Wallengren）。—Wallengren（1872, 44）。Pinacopteryx severina Stoll, 1781.—Wallengren（1875, 90）。[权威名称错误地给出为‘Cram’，实际上是指Cramer, 1781]。Pieris mesentina var. syrinx（Wallengren）。—Aurivillius（1879, 44）。Pieris mesentina Cramer，[1780]。—Trimen和Bowker（1889, 60），Aurivillius（1898, 408），Verity（1905–1911, 128）。Belenois aurota（Fabricius, 1793）。—Bridges（1988, 292），Ackery等人（1995, 200）。模式产地：‘Ad Swakop Africae’[纳米比亚]。模式标本：正模♀，hel