Run and output files from: Holocene population expansion of a tropical bee coincides with early human colonisation of Fiji rather than climate change

There is substantial debate about the relative roles of climate change and human activities on biodiversity and species demographies over the Holocene. In some cases, these two factors can be resolved using fossil data, but for many taxa such data are not available. Inferring historical demographies of taxa has become common, but the methodologies are mostly recent and their shortcomings often unexplored. The bee genus Homalictus is developing into a tractable model system for understanding how native bee populations in tropical islands have responded to past climate change. We greatly expand on previous studies using sequences of the mitochondrial gene COI from 474 specimens and between 171 and 3,928 autosomal (DArTSeq) SNP loci from 19 specimens of the native Fijian bee, Homalictus fijiensis (Perkins & Cheesman, 1928), to explore its historical demography using coalescent and mismatch analyses. We ask whether past changes in demography were human- or climate-driven, while considering analytical assumptions. We show that inferred changes in population sizes are too recent to be explained by past climate change. Instead we find that a dramatic increase in population size for the main island of Viti Levu coincides with increasing occupation by humans and their modification of the environment. We found no corresponding change in bee population size for another major island, Kadavu, where human populations and agricultural activities have been historically very low. Our analyses indicate that molecular approaches can be used to disentangle the impacts of humans and climate change on a major tropical pollinator and that stringent analytical approaches are required for reliable interpretation of results. 

Methods
Sampling sites and methods

Collections were made throughout Fiji between 2010 and 2017 from multiple localities but with the greatest number of samples from the largest island of Viti Levu (n = 309) and then the island of Kadavu (n = 109) (Table S1). Samples were collected from 3 to 1,328 m asl by sweep netting both native and introduced plants, and from nesting aggregations along roadsides. For each collection site, latitude, longitude and elevation were recorded using GPS devices (primarily using a Garmin 550). Once collected, bees were immediately transferred into individual vials of ³98% ethanol. Vials were kept at ~5˚C and ethanol was replaced within a week of collection to lessen DNA degradation.

Geographic information systems

In order to explore whether patterns in historical demography were related to subaerial land mass over time, we used bathymetric maps to examine how subaerial landmass and connectivity have changed since the last glacial maximum. Bathymetric data were downloaded in R version 3.6.2 using the package marmap version 1.0.4 (Pante & Simon-Bouhet, 2013). The marmap package was also used to produce maps and calculate subaerial landmasses presently and at the last glacial maximum. 

COI data generation

We subjected a subset of COI data from a previous study (Dorey et al., 2020b) to different analyses to answer novel hypotheses about past population demography of a single species, rather than relationships between many species. Tissue samples for DNA extraction were obtained by removing a single hind leg from each specimen. Samples prior to 2015 were sequenced at the Canadian Centre for DNA Barcoding (CCDB) at the Biodiversity Institute of Ontario (Groom et al., 2013). For these samples DNA amplification used the universal primer pair LepF1 and LepR2 (Groom et al., 2013; Hebert, Penton, Burns, Janzen, & Hallwachs, 2004). For all other samples DNA extractions and PCR amplifications were completed at the South Australian Regional Facility for Molecular Ecology and Evolution (SARFMEE) and DNA sequencing and purification carried out at Macrogen Inc. (Korea). DNA extractions at SARFMEE were performed using a Gentra Puregene® DNA Purification kit (Gentra Systems Inc.) according to the manufacturer’s protocol. PCRs amplified a 710 bp fragment of the mtDNA (COI) gene using the primers LCO1490 (forward) and HCO2198 (reverse). The 25 µL PCR reactions comprised the following reagents: Sterile H2O (15.9 µL), MRT buffer (5 µL), 1 µL (5 µM) of LCO1490, 1 µL (5 µM) of HCO2198, Immolase Taq (0.1 µL) and mtDNA from specimen (2 µL). The thermocycling regime comprised of one cycle at 94˚C for 10 minutes, then five cycles at 94˚C for 60 seconds, 46˚C for 90 seconds, 72˚C for 75 seconds, followed by 35 cycles at 94˚C for 60 seconds, 51˚C for 90 seconds, 72˚C for 75 seconds, followed by 72˚C for 10 minutes and then 25˚C for 2 minutes. 

Sequences were checked against the NCBI database using BLAST (blastn and blastx) to screen for non-HomalictusDNA. Forward and reverse sequences of each H. fijiensis specimen were aligned and checked for stop codons and/or nucleotide assignment errors using chromatograms examined with Geneious version 10.2.2 (Kearse et al., 2012). Any sequences with one or more base pairs that could not be reliably determined were excluded from the dataset. The H. fijiensis alignment was trimmed to 630 bp to remove primers and avoid spurious results that could arise from missing data in mismatch and Bayesian skyline coalescence analyses (Grant, 2015; Joly, Stevens, & Jansen van Vuuren, 2007).A total of 474 H. fijiensis sequences were analysed from across the entire Fijian archipelago including 309 sequences from the largest island of Viti Levu.

SNP data generation

We subjected the raw SNP data from a previous study (Dorey et al., 2020b) to more rigorous filtering and analyses that resulted in substantially changed subsets of the initial dataset that are more relevant to the present questions. The thorax and front legs were taken from 19 Viti Levu females from each of five species: H. fijiensis, H. tuiwawae, H. ostridorsum, H. groomi and H. sp. S (Dorey et al., 2020a; Dorey et al., 2019). We used the solid-state method Diversity Arrays Technology in Canberra, Australia (DArTseq™) (Jaccoud, Peng, Feinstein, & Kilian, 2001) to perform restriction site-associated DNA sequencing. DArTseq combines complexity reduction with a next generation sequencing platform in a conceptually-similar method to double digest RADseq (Shams et al., 2019). The restriction enzymes PstI and Hpall were used.

A total of 62,426 SNP loci were called across all five sequenced species. We used the R package DArTR version 1.3.4 (Gruber, Unmack, Berry, & Georges, 2018) to filter our data. The original SNP dataset was filtered to only include H. fijiensis. Monomorphic (non-variable) sites were then removed, leaving 7,719 loci. We then filtered these data to remove all missing data (4,046 loci remaining), for repeatability (percentage of scores that are repeated in the technical replicate dataset; 3,928 loci remaining) and to remove secondaries (multiple linked SNPs per fragment; 3,768 loci remaining). Genome-wide SNPs can suffer from large numbers of linked loci and this linkage can break assumptions of independence for many analyses, and bias results (O'Leary, Puritz, Willis, Hollenbeck, & Portnoy, 2018). Hence, we analysed a wide variety of linkage disequilibrium (LD) filtering criteria. We filtered for linkage, removing loci with r2values below 0.9, 0.7 and 0.2, retaining 1,811, 1,646 and 171 loci, respectively. Linked loci were removed sequentially in order of most- to least-linked connections to retain loci that might otherwise be removed (Script S1). We used the latter five filtering levels in analyses.

Haplotype analyses of COI data

A minimum-spanning network (Bandelt, Forster, & Röhl, 1999) of our complete COI dataset was created using PopART version 1.7 (Leigh & Bryant, 2015). Geneious version 10.2.2 (Kearse et al., 2012) was used to examine unique haplotypes and amino acid sequences.

Mismatch analyses of COI data

Mismatch analyses, and extended Bayesian skyline plots (EBSPs), make several assumptions about the data provided, including that: (i) a random sample is drawn from the population, (ii) the population is panmictic, and (iii) largely neutral markers were used (Grant, 2015). To examine panmixia, we used Arlequin version 3.11 (Excoffier, Laval, & Schneider, 2005) to examine pairwise FST values of COI haplotypes between all islands and combined island datasets that were not significantly different (p > 0.05). We then carried out mismatch analyses to explore whether past demographic changes could be explained by population expansion towards the present, graphing observed pairwise nucleotide differences with those expected under a recent population expansion, with 2,000 simulations used to generate an expected distribution of nucleotide differences. A unimodal distribution in a mismatch graph can be consistent with a sudden population expansion, whilst multimodal distributions can suggest past population bottlenecks or demographic structure (Rogers & Harpending, 1992).

Extended Bayesian skyline plots (EBSP) of COI data

We employed PartitionFinder version 2.1 (Lanfear, Frandsen, Wright, Senfeld, & Calcott, 2016) to determine the most appropriate model of molecular evolution for all COI datasets. Because the existence of population structure violates assumptions of panmixia (Heller, Chikhi, & Siegismund, 2013) our all-islands dataset was not used in our final demographic Bayesian analyses. For the island groups with sample sizes >50 (see Appendix) — Viti Levu and Kadavu — we used extended Bayesian skyline plots of COI data sets in BEAST version 2.6.3 (Bouckaert et al., 2019a; Drummond, Suchard, Xie, & Rambaut, 2012) to infer changes in historical demography. We restricted demographic EBSP analyses to the third codon position, where most synonymous mutations occur (Ermolaeva, 2001). We applied a strict molecular clock and the best-fit PartitionFinder 2 model, an HKY and an HKY+Γ substitution model for the Viti Levu and Kadavu populations, respectively. Four independent runs for each analysis were performed, to confirm stationarity. For the Viti Levu population, each run consisted of 4 chains with heating, carried out for 300 million iterations, resampling every 30,000th iteration using the BEAST package CoupledMCMC version 1.0.2 (Altekar, Dwarkadas, Huelsenbeck, & Ronquist, 2004). Multiple chains were required to properly sample across multiple possible optima in phylospace. For the Kadavu population, we used single chains with 500 million iterations sampled every 100,000th iteration. The log files from each run were examined in Tracer version 1.7.1 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018) and a burnin of 10% was used, which was always after stationarity had been achieved (effective sample sizes all exceeding 200). Log and tree files were combined using LogCombiner version 2.5.0 (Bouckaert et al., 2019a). The EBSP log files were analysed with the plotEBSP script in R version 3.5 (Bouckaert et al., 2019a).

The estimated mutation rate of 1.09x10-7 per site per generation, based on only the 3rd codon position from the whole mitogenome of Caenorhabditis elegans (Denver, Morris, Lynch, Vassilieva, & Thomas, 2000) was applied to all of our EBSP plots to infer an approximate time scale in generation units. This directly-estimated mutation rate is appropriate for inferring recent demographic changes and is broadly consistent with other empirical values (Gratton, Konopiński, & Sbordoni, 2008; Papadopoulou, Anastasiou, & Vogler, 2010). The AT bias of C. elegans (70.3% from 21 C. elegansCOI sequences; BOLD (2020)) is similar to that of our H. fijiensis COI fragment (74% from 474 H. fijiensis COI sequences). We converted the EBSP time scale from generations to chronological time by assuming four Homalictus generations per year, following Groom et al. (2014). However, we also explored the effects of assuming three or five generations per year, which is analytically equivalent to assuming a faster or slower per-generation mutation rate. 

For our reconstruction of island ancestral states, we analysed all three codon positions using within-island unique COI sequences and the BEAST 2 package CoupledMCMC with eight chains and an EBSP tree prior. Four independent runs were performed, to confirm stationarity. We used two outgroup species — H. groomi and H. sp. O. The molecular data were allocated to a single partition to which we applied a single HKY+I substitution model and an uncorrelated relaxed clock model. "Island" was included in the analysis as a discrete trait given a symmetric change model and a strict clock (more complex models prevented this partition from converging). Each run was 100 million iterations, resampling every 20,000th iteration. Log files were examined in Tracer and a burnin of 10% was used. The maximum clade credibility tree as well as posterior support values were produced in TreeAnnotator version 2.6.3 (Bouckaert et al., 2019b) using median node heights. The tree was visualised in FigTree version 1.4.4 (Drummond, 2016).

Extended Bayesian skyline plots and Ne using SNP data

We employed PartitionFinder 2 to determine the most appropriate model of molecular evolution for each SNP dataset corresponding to the different levels of linkage filtering. For our SNP datasets, we ran EBSP analyses in BEAST 2 using the following models: for our LDR2=0.2 we used a K80 model, for LDR2=0.7 and LDR2=0.9 we used HKY+Γ and for the remaining datasets (without and with secondaries — multiple SNPs on a single fragment — included) we used a GTR+Γ model. We used a relaxed log normal clock for all SNP EBSP analyses (Lemey, Rambaut, Welch, & Suchard, 2010). All runs except for that with secondaries were executed for 100 million iterations, sampling every 50,000thiteration, and were repeated four times to confirm convergence in EBSP results. For the run that kept secondaries, four heated-chain runs were carried out for 100 million iterations, resampling every 10,000th iteration using the BEASTpackage CoupledMCMC. The log files from each run were examined in Tracer and then combined using LogCombinerwith a burnin of 20%, which was always after stationarity had been achieved (effective sample sizes all exceeding 200). The EBSP log files were analysed using plotEBSP (Bouckaert et al., 2019a). 

Supplementary methods summary

We examined haplotype sample sizes required for robust demographic inference, using rarefaction analysis in EstimateS version 9.1.0 (Colwell, 2013) (see Appendix). We also undertook nested sampling (NS package in BEAST 2) and used DIYABC-RF version 1.0.12 in R (Collin et al., Submitted; Durif & Collin, 2021) to explicitly compare alternative possible demographic patterns (see Appendix).

学界围绕全新世以来气候变化与人类活动对生物多样性及物种种群动态的相对作用尚存广泛争议。部分类群可借助化石数据厘清二者的影响，但多数类群并无相关化石数据。通过分子手段推断类群的历史种群动态已成为常用方法，但相关分析方法多为近年提出，其局限性往往未得到充分探讨。熊蜂属（Homalictus）正逐步成为探究热带岛屿本土蜂类种群对历史气候变化响应的易操作模式系统。本研究大幅拓展了前期研究的范围，针对斐济本土熊蜂——斐济熊蜂（Homalictus fijiensis，Perkins & Cheesman, 1928），利用474份标本的线粒体基因细胞色素c氧化酶亚基I（COI）序列，以及19份标本的171至3928个常染色体多样性芯片技术（DArTseq）单核苷酸多态性（SNP）位点数据，通过溯祖分析与错配分布分析探究其历史种群动态。研究在考量分析假设的前提下，探讨种群动态历史变化是由人类活动还是气候变化驱动。结果显示，推断的种群规模变化时间过晚，无法用历史气候变化解释。相反，维提岛（斐济最大岛屿）的种群规模显著扩张，与人类在此区域的定居及环境改造进程高度吻合；而在人类活动与农业开发历史上均极低的坎达武岛，并未观测到蜂类种群规模的相应变化。本研究表明，分子手段可用于厘清人类活动与气候变化对热带主要传粉昆虫的影响，且需采用严谨的分析方法才能得到可靠的结果。

## 研究方法

### 采样地点与方法
2010至2017年间，研究团队在斐济全境多个点位开展采样，其中样本量最大的为维提岛（n=309），其次为坎达武岛（n=109）（详见附表S1）。采样覆盖海拔3至1328米的区域，采用扫网法采集本土与外来植物上的蜂类，同时采集路边筑巢集群的个体。每个采样点均通过GPS设备（主要使用Garmin 550型GPS）记录纬度、经度与海拔信息。采集后的蜂类立即转移至装有98%乙醇的单独离心管中，将离心管置于约5℃环境保存，并在采样后一周内更换乙醇以降低DNA降解程度。

### 地理信息系统分析
为探究种群动态历史模式是否与历史时期的露海陆地面积相关，研究团队借助水深图分析末次盛冰期以来露海陆地面积与连通性的变化。水深数据通过R语言3.6.2版本中的`marmap`包（版本1.0.4）下载（Pante & Simon-Bouhet, 2013），并使用该包绘制地图并计算当前及末次盛冰期的露海陆地面积。

### COI数据生成
本研究针对前期研究（Dorey et al., 2020b）中的部分COI数据开展不同分析，以验证针对单一物种种群动态历史的全新假说，而非针对多物种间演化关系的分析。DNA提取的组织样本取自每头标本的单条后足。2015年之前采集的样本在安大略省生物多样性研究所的加拿大DNA条形码中心（CCDB）完成测序（Groom et al., 2013），此类样本的DNA扩增使用通用引物对LepF1与LepR2（Groom et al., 2013; Hebert et al., 2004）。其余样本的DNA提取与PCR扩增均在南澳大利亚分子生态与进化区域设施（SARFMEE）完成，DNA测序与纯化则由韩国Macrogen公司完成。SARFMEE的DNA提取使用Gentra Puregene® DNA纯化试剂盒（Gentra Systems Inc.），严格遵循制造商的操作流程。PCR扩增使用引物LCO1490（正向）与HCO2198（反向），扩增线粒体DNA（mtDNA）中COI基因的710 bp片段。25 μL的PCR反应体系包含以下试剂：无菌去离子水（15.9 μL）、MRT缓冲液（5 μL）、1 μL（5 μM）LCO1490引物、1 μL（5 μM）HCO2198引物、Immolase Taq酶（0.1 μL）以及标本的mtDNA模板（2 μL）。热循环程序设置为：94℃预变性10分钟；随后5个循环：94℃变性60秒、46℃退火90秒、72℃延伸75秒；再进行35个循环：94℃变性60秒、51℃退火90秒、72℃延伸75秒；最后72℃终延伸10分钟，25℃保温2分钟。

### COI序列质控与比对
研究团队使用BLAST（blastn与blastx）将所有序列与NCBI数据库比对，以筛选非Homalictus属的非目标序列。每头斐济熊蜂（H. fijiensis）标本的正向与反向序列经比对后，使用Geneious 10.2.2版本（Kearse et al., 2012）查看色谱图，以检查终止密码子及核苷酸赋值错误。所有存在无法可靠判定的单个或多个碱基对的序列均被排除出数据集。为去除引物序列并避免错配分布分析与贝叶斯天际线溯祖分析中因缺失数据导致的伪影结果，将斐济熊蜂的序列比对结果裁剪至630 bp（Grant, 2015; Joly et al., 2007）。最终，研究团队对覆盖整个斐济群岛的474条斐济熊蜂序列开展分析，其中包括来自维提岛的309条序列。

### SNP数据生成
本研究对前期研究（Dorey et al., 2020b）中的原始单核苷酸多态性（SNP）数据开展更为严格的过滤与分析，最终得到与当前研究问题更相关的、经大幅优化的初始数据集子集。研究采集了5个物种各19头来自维提岛的雌蜂，分别为斐济熊蜂（H. fijiensis）、H. tuiwawae、H. ostridorsum、H. groomi以及H. sp. S（Dorey et al., 2020a; Dorey et al., 2019），取样组织为胸部与前足。研究使用澳大利亚堪培拉的多样性芯片技术（DArTseq™）开展限制性酶切位点相关DNA测序（Jaccoud et al., 2001）。DArTseq技术通过降低基因组复杂度结合二代测序平台，其原理与双酶切RADseq技术相似（Shams et al., 2019），本次分析使用的限制性内切酶为PstI与HpaII。

本次测序在5个物种中共检出62426个SNP位点。研究使用R语言的`DArTR`包（版本1.3.4，Gruber et al., 2018）对数据进行过滤：首先将原始SNP数据集限定为仅包含斐济熊蜂；随后移除单态（非变异）位点，最终保留7719个位点；接着依次过滤缺失数据（保留4046个位点）、筛选技术重复数据集的可重复性得分（保留3928个位点），并移除次级位点（即单个片段上存在多个连锁SNP的位点，最终保留3768个位点）。全基因组SNP数据常存在大量连锁位点，此类连锁关系会破坏诸多分析的独立性假设并导致结果偏倚（O'Leary et al., 2018）。因此，研究团队针对连锁不平衡（LD）筛选了多种标准：按r²阈值分别为0.9、0.7与0.2过滤连锁位点，依次保留1811、1646与171个位点。移除连锁位点时，按照连锁程度从高到低的顺序依次剔除，以保留可能被误删的位点（详见脚本S1）。最终分析中使用了后5种过滤水平的数据集。

### COI数据单倍型分析
研究团队使用PopART 1.7版本（Leigh & Bryant, 2015）构建完整COI数据集的最小生成网络（Bandelt et al., 1999），并使用Geneious 10.2.2版本（Kearse et al., 2012）分析独特单倍型与氨基酸序列。

### COI数据错配分布分析
错配分布分析与扩展贝叶斯天际线图（EBSP）分析均对输入数据存在若干假设，包括：（1）样本为种群的随机抽样；（2）种群为泛交种群；（3）所使用的标记基本符合中性进化假设（Grant, 2015）。为验证泛交假设，研究使用Arlequin 3.11版本（Excoffier et al., 2005）计算所有岛屿及合并岛屿数据集的COI单倍型成对FST值，筛选出无显著差异的分组（p>0.05）。随后开展错配分布分析，以探究种群动态历史变化是否符合种群向现代的扩张模式：将观测到的成对核苷酸差异与近期种群扩张模型下的预期差异进行绘图，并通过2000次模拟生成预期核苷酸差异分布。错配分布图中的单峰分布可对应种群的突然扩张，而多峰分布则提示历史上存在种群瓶颈或种群结构分化（Rogers & Harpending, 1992）。

### COI数据扩展贝叶斯天际线图（EBSP）分析
研究使用PartitionFinder 2.1版本（Lanfear et al., 2016）为所有COI数据集确定最适的分子进化模型。由于种群结构会违反泛交假设（Heller et al., 2013），因此全岛屿数据集未纳入最终的贝叶斯种群动态分析。针对样本量大于50的岛屿类群（详见附录）——维提岛与坎达武岛，研究使用BEAST 2.6.3版本（Bouckaert et al., 2019a; Drummond et al., 2012）中的COI数据集扩展贝叶斯天际线图分析，以推断其种群动态历史变化。分析限定在第三密码子位点开展，因为该位点是同义突变的主要发生区域（Ermolaeva, 2001）。我们应用严格分子钟与PartitionFinder 2确定的最适模型：维提岛种群使用HKY替换模型，坎达武岛种群使用HKY+Γ替换模型。每项分析均开展4次独立运行以确认平稳性：维提岛种群的每次运行包含4条加热链，共进行3亿次迭代，每30000次迭代抽样一次，使用BEAST的`CoupledMCMC`包（版本1.0.2，Altekar et al., 2004）完成，多链设置可确保在系统发育空间的多个最优解间充分抽样。坎达武岛种群则使用单链，进行5亿次迭代，每100000次迭代抽样一次。所有运行的日志文件使用Tracer 1.7.1版本（Rambaut et al., 2018）检查，均在达到平稳性后舍弃10%的预烧期样本（有效样本量均大于200）。使用LogCombiner 2.5.0版本（Bouckaert et al., 2019a）合并日志与树文件，并使用R 3.5版本中的`plotEBSP`脚本（Bouckaert et al., 2019a）分析EBSP日志文件。

本研究以秀丽隐杆线虫（Caenorhabditis elegans）全线粒体基因组第三密码子位点的突变率（1.09×10^-7 每个位点每世代，Denver et al., 2000）作为所有EBSP分析的突变率，以推代为单位构建近似时间尺度。该直接测定的突变率适用于近期种群动态变化的推断，且与其他经验测定值基本一致（Gratton et al., 2008; Papadopoulou et al., 2010）。秀丽隐杆线虫的AT碱基偏好性（基于21条秀丽隐杆线虫COI序列，为70.3%；BOLD, 2020）与本研究中斐济熊蜂COI片段的AT偏好性（基于474条序列，为74%）相似。参考Groom et al.（2014）的方法，研究假设Homalictus属蜂类每年有4个世代，将EBSP分析得到的世代时间尺度转换为年代时间尺度。同时，研究也探索了假设每年3或5个世代的影响，这在分析上等价于假设更高或更低的每世代突变率。

### 岛屿祖先状态重建
为重建岛屿祖先状态，研究使用各岛屿内的独特COI序列，结合BEAST 2的`CoupledMCMC`包开展分析，设置8条链与EBSP树先验。共开展4次独立运行以确认平稳性，外类群选用H. groomi与H. sp. O。分子数据划分为单个分区，应用HKY+I替换模型与无关联松弛分子钟模型。将“岛屿”作为离散性状纳入分析，采用对称变换模型与严格分子钟（更复杂的模型无法使该分区收敛）。每次运行包含1亿次迭代，每20000次迭代抽样一次。日志文件使用Tracer检查，舍弃10%的预烧期样本。使用TreeAnnotator 2.6.3版本（Bouckaert et al., 2019b）以节点高度中位数构建最大类群置信度树并计算后验支持率。使用FigTree 1.4.4版本（Drummond, 2016）可视化系统发育树。

### 基于SNP数据的扩展贝叶斯天际线图与有效种群大小分析
研究使用PartitionFinder 2为不同连锁过滤水平的SNP数据集分别确定最适分子进化模型。针对本研究的SNP数据集，在BEAST 2中开展EBSP分析时采用如下模型：r²阈值为0.2的数据集使用K80模型，r²阈值为0.7与0.9的数据集使用HKY+Γ模型，其余数据集（包含或不包含次级位点）使用GTR+Γ模型。所有SNP的EBSP分析均使用松弛对数正态分子钟（Lemey et al., 2010）。除保留次级位点的数据集外，其余所有运行均进行1亿次迭代，每50000次迭代抽样一次，且重复4次以确认EBSP结果收敛。保留次级位点的运行则使用BEAST的`CoupledMCMC`包开展4条加热链运行，进行1亿次迭代，每10000次迭代抽样一次。所有运行的日志文件使用Tracer检查，随后使用LogCombiner合并，舍弃20%的预烧期样本（所有有效样本量均大于200，且均已达到平稳性）。最终使用`plotEBSP`脚本（Bouckaert et al., 2019a）分析EBSP日志文件。

### 补充方法概述
研究使用EstimateS 9.1.0版本（Colwell, 2013）中的稀疏分析法，探究可靠种群动态推断所需的单倍型样本量（详见附录）。同时，研究开展嵌套抽样（BEAST 2中的NS包），并使用R语言中的DIYABC-RF 1.0.12版本（Collin et al., 已投稿; Durif & Collin, 2021）对多种潜在种群动态模式开展显式比较（详见附录）。