# Bayesian Total-Evidence Dating Reshapes the Age and Historical Biogeography of Parnassiinae

Noémie M.-C. Hévin, Fabien L. Condamine

► To cite this version:

Noémie M.-C. Hévin, Fabien L. Condamine. Bayesian Total-Evidence Dating Reshapes the Age and Historical Biogeography of Parnassiinae. Bulletin of the Society of Systematic Biologists, 2024, 3 (2), 10.18061/bssb.v3i2.9882. hal-04896879

HAL Id: hal-04896879

https://hal.science/hal-04896879v1

Submitted on 20 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bayesian Total-Evidence Dating Reshapes the Age and Historical Biogeography of Parnassiinae

Noémie M-C. Hévin¹, Fabien L. Condamine¹

¹ CNRS, Institut des Sciences de l'Evolution de Montpellier, Université de Montpellier, Place Eugène Bataillon, 34095 Montpellier, France
https://doi.org/10.18061/bssb.v3i2.9882

Bulletin of the Society of Systematic Biologists

## Abstract

Although incomplete, the fossil record offers direct evidence for the existence of a lineage, providing insights into its age and geographic location. Reconstructing time-calibrated phylogenies, including both extant and fossil taxa as lineages (total-evidence dating) under the fossilized birth-death process, can provide new information about the phylogenetic position of fossil specimens, the time at which they diverged from closely related species, and can impact inferences of clades' historical biogeography. Here, we focus on the origin and radiation of Apollo butterflies (Papilionidae: Parnassiinae), whose origin has been estimated in the late Eocene in Central Asia. The two fossil taxa recovered in the subfamily Parnassiinae, dated to the late Oligocene (†*Thaites ruminiana* Scudder) and late Miocene (†*Doritites bosniaskii* Rebel) in the Western Palearctic, thus challenge the Central Asian origin. We performed a Bayesian total-evidence dating approach to explore the impact of dating analyses and past fossil distributions on the estimation of the evolutionary history of the group. We obtained a more credible dating and historical biogeography for the group, placing its origin in the late Paleocene (ca. 57 Ma) in the Western Palearctic + Western Asia and Caucasus regions, followed by dispersals to Central Asia and the Himalayas. This study also highlights the importance of investigating fossil position in addition to clock partitioning and models for molecular dating. Furthermore, we confirm that a few fossils are sufficient to cast doubt on the origin and biogeographic history of a group, especially when those fossils were dispersed outside of the current center of diversity.

## INTRODUCTION

Dated phylogenies are widely used since they provide the basis for testing many hypotheses in ecology and evolutionary biology (Zhang et al., 2016). With the increasing amount of molecular data and methods available, the question of whether fossil data are even needed has been asked (e.g. Monroe & Bokma, 2010; Venditti et al., 2011). While studying the body size evolution of mammals, Monroe & Bokma (2010) excluded fossils to avoid artefacts, and Venditti et al. (2011) argued that the inclusion of fossils is unlikely to alter trends in extant species, as their method accounts for different rates of change in extinct groups. However, given that molecular divergence can only provide a relative timeline, calibration with data from an external source remains necessary (Heath et al., 2014; Ho & Phillips, 2009; Rieux & Balloux, 2016). It is widely recognized that fossils offer unique direct evidence for the existence of a lineage (Mongiardino Koch et al., 2021; Mongiardino Koch & Parry, 2020) and information on the absolute timing of the lineage's diversification, providing insights into its age,

preserved traits, and geographic location (Pett & Heath, 2020).

Two powerful and complementary methods represent important conceptual advances for estimating divergence times by directly integrating fossils: the total-evidence dating (Pyron, 2011; Ronquist et al., 2012a), and the fossilized birth-death process (Heath et al., 2014). Total-evidence dating (TED) consists of taking into account molecular and morphological characters to put together in a single dataset both extant and extinct species. Contrary to node dating, TED directly integrates fossils in the phylogeny and uses all fossils and not only the oldest one (Ronquist et al., 2012a; Matschiner, 2019). However, O'Reilly and Donoghue (2016) showed that TED and node dating should be used in a complementary way. Early attempts of total-evidence dated trees typically show exaggerated divergence time estimates (O'Reilly et al., 2015). In fact, a uniform tree prior was used, which assigned equal weight to trees even if they were less likely and treated all possible distributions of branching periods as equally probable (Hunt & Slater, 2016). Instead of the uniform prior, the fossilized birth-death (FBD) process

<figure><img src="image_2.png" /></figure>

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (CC BY-4.0). View this license's legal deed at http://creativecommons.org/licenses/by/4.0 and legal code at http://creativecommons.org/licenses/by/4.0/legalcode for more information.
is a mechanistic model that allows fossil tips and sampled
ancestors while providing a model of lineage diversification
by taking into account speciation, extinction, fossilization,
and taxon sampling (Gavryushkina et al., 2017; Heath et al.,
2014). Total-evidence dating under the FBD process thus
makes it possible to use the morphology but also the strati-
graphic age of the fossils as evidence to infer both phy-
logenetic relationships and divergence times (Donoghue
& Yang, 2016; Luo et al., 2020; Mongiardino Koch et al.,
2021).

Total-evidence dating provides a major advance in the field of dating, but it comes with new challenges ranging from the incompleteness of the fossil record to the morphological matrix and the clock coming with it as morphological characters rarely evolve in a clock-like way (Barba-Montoya et al., 2021; Luo et al., 2020; Mongiardino Koch et al., 2021). The most concerning is that older time estimates are often inferred with a TED approach (Arcila et al., 2015; Barba-Montoya et al., 2021; Casali et al., 2020; Hunt & Slater, 2016; Renner et al., 2016; but see Ronquist et al., 2012a; Gavryushkina et al., 2017), a phenomenon qualified as ‘deep root attraction’ (Luo et al., 2020; Ronquist et al., 2016). Several hypotheses have been proposed about what can influence this attraction: for example the failure to account for diversified sampling (Ronquist et al., 2016; Zhang et al., 2016), inadequate morphological models (Ronquist et al., 2016), the morphological matrices often coming with large amount of missing data (Ronquist et al., 2012; Casali et al., 2020; Hunt & Slater, 2016; Matschiner, 2019; Mongiardino Koch et al., 2021; O’Reilly & Donoghue, 2021; Ronquist et al., 2016), or the use of fixed node aged or hard bounds instead of soft bounds (Barido-Sottani et al., 2019, 2020; Ho & Phillips, 2009; Yang & Rannala, 2006). The precision of time can also be influenced by the proximity of the fossil taxa with the dated nodes, as well as by the root priors (Arcila et al., 2015). The effects of partitioning (Angelis et al., 2018), clock models (Ronquist et al., 2012a; Gavryushkina et al., 2017; Zhang et al., 2016) and choice of outgroups (Spasojevic et al., 2021) can also impact the divergence time estimates, and have not been investigated by many TED studies yet.

Simulation and empirical studies have shown that fossil
data can improve inferences about trait evolution (Hunt &
Slater, 2016). Total-evidence dating approaches are poten-
tially able to detect homoplasy and thus improve inferences
about phylogenetic relationships, character evolution and
divergence times (Lee & Yates, 2018; Tejada et al., 2024).
More and more studies have relied on a TED approach for
empirical studies to understand the evolutionary history of
an array of groups (Darlim et al., 2022; Heritage & Seif-
fert, 2022; Jiangzuo & Flynn, 2020; Lee & Yates, 2018;
Matschiner et al., 2017). Groups with only a few extant rep-
resentatives but a lot of fossils like Sirenians (Heritage &
Seiffert, 2022), Xenarthrans (Tejada et al., 2024), or groups
with a rich fossil record like fishes (Near & Kim, 2021),
penguins (Gavryushkina et al., 2017), and turtles (Selvatti
et al., 2023) have been investigated. However, more recent
studies have shown that the use of a few fossils can have
significant effect on topologies, but also on the historical

biogeography scenarios, especially when fossils are recov-
ered in a different distribution than the one known for the
extant species (Bacon et al., 2022; Coiro et al., 2023; Parks
et al., 2022; Rule et al., 2020; Santos et al., 2022; Wood et
al., 2013). From a macroevolutionary perspective, direct in-
corporation of fossils through TED under the FBD process
is a valuable method to better understand the evolutionary
history of the group that can improve the timing of diversi-
fication (Coiro et al., 2023; Mongiardino Koch et al., 2021;
Near & Kim, 2021; Wisniewski et al., 2022), and should be
preferred over node dating when sufficient morphological
data are available for both fossil and extant taxa (e.g. 100
characters at least; Pyron, 2011; Ronquist et al., 2012a).
Globally, it is important to keep in mind that each case is
potentially unique, which needs to be considered to draw
some conclusions (Ronquist et al., 2016).

Biological model system

Another biological model to investigate the impact of a
few fossils on the dating and historical biogeography of a
group is the Apollo butterflies (Papilionidae: Parnassiinae),
which consists of about >85 recognized species distributed
in the Holarctic. Two fossil taxa are unambiguously recov-
ered in the subfamily Parnassiinae (Nazari et al., 2007):
†Doritites bosniaskii Rebel and †Thaites ruminiana Scudder.
For 35 extant species and the two fossils, a morphological
matrix containing 236 characters has been made by Nazari
et al. (2007), later used by Condamine et al. (2018) to in-
vestigate the phylogenetic position of these fossils with
the extant species by reconstructing a phylogeny using a
total-evidence approach (without dating). †Doritites bosni-
askii is recovered as sister to the genus Archon Hübner,
while †Thaites ruminiana is found to be sister to the tribe
Parnassiini. However, it is well-known that fossil ages are
also used as a source of evidence in TED to infer the rela-
tion between extinct and extant taxa (Barido-Sottani et al.,
2020). Furthermore, the origin of the subfamily has been
estimated to be in the late Eocene (ca. 38 Ma) in Cen-
tral Asia (Condamine et al., 2018) but it highlighted a high
uncertainty in the ancestral estimation of intergeneric re-
lationships which represent the backbone and have long
branches that are difficult to estimate. Dated to the late
Oligocene (†Thaites ruminiana) and late Miocene (†Doritites
bosniaskii) in the western Palearctic, the two fossils thus
challenge the Central Asian origin and can help clarify the
ancestral states of the backbone.

Objectives

In this study, we investigate the impact of fossils on the
dating analyses and estimation of the historical biogeog-
raphy using the Apollo butterflies. To compare with Con-
damine et al. (2018), we performed a Bayesian total-ev-
idence dating approach under the fossilized birth-death
process based on the same dataset (four mitochondrial
genes, one nuclear gene, and 180 morphological characters)
and concentrated on the backbone (intergeneric relation-
ships). Specifically, we (i) assessed the effect of dating phy-
logenies on fossil placements and node supports; (ii) car-
<figure><img src="image_3.png" /><figcaption>Figure 1. Original illustrations of the two fossils belonging to Parnassiinae. Illustrations of †Doritites bosniaskii on the left were made by Rebel (1898), and illustrations of †Thaites ruminiana on the right were made by Scudder (1875).</figcaption></figure>

ried out total-evidence dating under the fossilized
birth-death process to examine the effect of clock parti-
tioning (from 1 to 12 partitions) and clock models (auto-
correlated, uncorrelated, and strict clocks) on the node age
estimates; and (iii) inferred parametric time-stratified bio-
geographic analyses with the resulting total-evidence dated
trees, thus considering the past western Palearctic distrib-
utions of the two fossils.

METHODS

Taxon sampling, molecular and
morphological datasets

Parnassiinae fossils

Within Parnassiinae, †*Thaites ruminiana* (Scudder, 1875) (Fig. 1) is the first fossil that undoubtedly belongs to the subfamily. It is a compression fossil from limestone in the “Niveau du gypse d’Aix Formation” of France (Bouches-du-Rhône, Aix-en-Provence) within the Chattian (23.03–28.1 Ma) of the late Oligocene (Rasnitsyn & Zherikhin, 2002; Sohn et al., 2012). †*Thaites ruminiana* was found to be most likely the sister of Parnassiini (Condamine et al., 2018), however it is also occasionally found to be the sister of Luehdorfini and Sericinini. The second fossil is †*Doritites bosniaskii* (Rebel, 1898) (Fig. 1), an exoskeleton and compression fossil from the Messinian stage (5.33–7.25 Ma, late Miocene) that was found in Tuscany, Italy (Rebel, 1898). In accordance with Carpenter’s (1992) findings, Condamine et al. (2018)’s study strongly supports †*Doritites bosniaskii* as sister to the genus *Archon*.

Molecular dataset

The dataset from Condamine et al. (2018) was used. It con-
tains the 70 species currently recognized by morphology
(Weiss, 1991–2005), 15 lineages that can be recognized as
valid species (Condamine et al., 2018), the two fossils un-
ambiguously belonging to Parnassiinae (†*Thaites ruminiana*
and †*Doritites bosniaskii*; Nazari et al., 2007), and nine out-
groups species (*Baronia brevicornis* Salvin, and eight species
belonging to the sister subfamily Papilioninae). The mol-
ecular dataset comprises 4,535 nucleotides with five genes
(four mitochondrial genes and one nuclear gene) for all
species (except the two fossils and *Bhutanitis ludlowi*
Gabriel): cytochrome oxidase I (COI), NADH dehydrogenase
1 (ND1), NADH dehydrogenase 5 (ND5), rRNA 16S (16S),
and elongation factor 1 alpha (EF1α). The summary statis-
tics from AMAS (Borowiec, 2016) recovered 12% of missing
data and 41% of variable sites in the molecular alignment.
Molecular partitioning was made using PartitionFinder
1.1.1 (Lanfear et al., 2012), which recovered eight partitions
for the molecular dataset.

Morphological dataset

The morphological dataset is a subsampling of the morpho-
logical dataset from Nazari et al. (2007), which now con-
tains only 180 characters for 35 extant species and the two
fossils. The subsampling was made to only retain the char-
acters for which at least one of the three species without
molecular data (the two fossils and *Bhutanitis ludlowi*) was
coded, and thus reduce issues coming with missing data
(Ronquist et al., 2016). In this matrix, morphological char-
acters were coded for all extant species, except in the genus
*Parnassius* Latreille, for which only one species per sub-genus was coded since there are no fossils related to this genus (Nazari et al., 2007). For one analysis (see below), the morphological dataset was divided into four partitions adapted from the body categories (body [1-53], wing basal structure and venation [54-106], wing pattern [107-158], and other [159-180]) determined by Nazari et al. (2007) (Appendix S1). Molecular and morphological datasets were then combined to construct a total-evidence dataset.

## Bayesian total-evidence dating

Bayesian TED approaches under the FBD process were carried out to estimate divergence times. Contrary to Condamine et al. (2018) who used BEAST, total-evidence phylogenetic and dating analyses were carried out using MrBayes 3.2.7a (Ronquist et al., 2012b) as it allows the use of multiple runs along with the reversible-jump Markov Chain Monte Carlo (rj-MCMC) option, as well as to test for autocorrelated and uncorrelated clocks. The rj-MCMC option was used to sample the entire substitution rate model space for each partition recovered by PartitionFinder (Huelsenbeck et al., 2004). The one-parameter Markov (Mk) model was used to model morphological data (Lewis, 2001). Three fossils were used for the calibrations: †*Doritites bosniaskii* and †*Thaites ruminiana* which are directly included in the matrix with uniform distributions (5.33–7.25 Ma and 23.03–28.1 Ma, respectively) and the genus †*Praepapilio* (Durden & Rose, 1978) from the early Lutetian (Eocene) of the Green River Formation (Colorado, U.S.A.), which is used to constrain the crown age of Papilionidae with a uniform distribution (47.8–100.5 Ma; de Jong, 2017; Smith et al., 2003). For each test (see below), two independent analyses with eight rj-MCMC, each was run for 50 million generations, sampled every 50,000 generations. A consensus tree was built after discarding 25% samples as burn-in. The nexus files with all parameters for all analyses following are available online.

## Effect of dating phylogenies on fossil placements and node supports

The placement of the two fossils belonging to Parnassiinae was tested by evaluating the impact of (i) dating, (ii) the choice of outgroups, and (iii) the level of missing data. To do that, a subsampling of the 53 extant and extinct species with morphological data was used to reduce missing data. For all fossil placement tests, one topological constraint was made bringing Papilioninae monophyletic, and a unique uncorrelated clock (IGR; Drummond et al., 2006) was used, as it has been identified as the best model fitting dataset including the fossil record (Ronquist et al., 2012a). To test for the effect of outgroups, *Baronia brevicornis* was removed because of its long branch, which can cause phylogenetic artifacts in Parnassiinae (Allio et al., 2020; Nazari et al., 2007). To test for the effect of missing data, subsampling matrices containing the characters of either †*Thaites ruminiana* (81 characters) or †*Doritites bosniaskii* (96 characters) are used. Finally, to test the effect of dating, all these tests were made by constructing the phylogeny only

or by simultaneously inferring and dating the phylogeny. Altogether, this study resulted in eight analyses as follows: (i) total-evidence phylogeny (TEP) with 180 characters and all outgroups, (ii) TEP with 180 characters and Papilioninae as outgroups only (*Baronia* removed), (iii) TEP with 81 characters (no missing data for †*Thaites ruminiana*) and all outgroups, (iv) TEP with 96 characters (no missing data for †*Doritites bosniaskii*) and all outgroups, (v) TED with 180 characters and all outgroups, (vi) TED with 180 characters and Papilioninae as outgroups only (*Baronia* removed), (vii) TED with 81 characters (no missing data for †*Thaites ruminiana*) and all outgroups, and (viii) TED with 96 characters (no missing data for †*Doritites bosniaskii*) and all outgroups. For each analysis, the R package RoguePlots (Klopfstein & Spasojevic, 2019; available at: https://github.com/seraklop/RoguePlots) is used to estimate the probabilities and robustness of the placement of each fossil.

## Effect of clock partitioning and clock models on the node age estimates

The impact of the data partitioning was tested on the dataset containing all 94 species. Two topological constraints were made to enforce the monophyly of Papilioninae and Parnassiinae, and the IGR clock model is used. Five total-evidence dating were made as follows: (i) a unique clock, (ii) one clock for the molecular dataset and one clock for the morphological dataset, (iii) one clock for the mitochondrial genes, one clock for the nuclear gene, and one clock for the morphological dataset, (iv) one clock for each partition from PartitionFinder, and one clock for the morphological dataset, and (v) one clock for each partition from PartitionFinder and one clock for each partition of the morphological dataset. To compare the analyses, estimates of the marginal likelihoods obtained using the stepping-stone sampling approach (Xie et al., 2011) were tried. However, as pointed out by Ronquist et al. (2012a), the inclusion of fossils resulted in the non-convergence of runs. To select the best model, the convergence of the clock was checked and the analyses were ranked based on harmonic mean estimators.

The impact of clock models on the total-evidence dated phylogenies was tested on the dataset containing all species, and two topological constraints are enforced on Papilioninae and Parnassiinae. Without accounting for the results of previous tests, it was chosen to use three clock models (IGR, autocorrelated [TK02], and strict model) for three partitions (one for the mitochondrial genes, one for the nuclear gene, and one for the morphological dataset), resulting in 12 analyses, considering that the strict clock model was tested only for the morphological dataset. As above, harmonic mean estimators were used to rank the analyses, with the first one being considered as the best.

## Final total-evidence dated phylogeny

The final total-evidence dated phylogeny is constructed with the best partition clock and the corresponding best clock models. Additionally, a node dating analysis is made
using the same minimum calibration times as the TED analysis and a conservative maximum age of 100.5 Ma for the three constraints, considering †*Thaites ruminiana* as crown of Parnassiinae, †*Doritites bosniaskii* as crown of Luehdorfiini, and †*Praepapilio* as constraint crown age of Parnilionidae, as it has already been done by Condamine et al. (2018) and Allio et al. (2020). For both analyses, two independent analyses with eight rj-MCMC, each was run for 50 million generations, sampled every 50,000 generations, resulting in 1,000 sampled trees and log files. For each analysis, a consensus tree was built after discarding 25% of the generations as burn-in to compute each clade's posterior probability (PP), median age, and 95% highest posterior density (HPD). Convergence of runs was checked with the potential scale reduction factors (PSRF), which should be close to 1.0, and then evaluated graphically by computing the effective sample size (ESS) of relevant parameters under Tracer 1.7.2 (Rambaut et al., 2018), using the recommended threshold of 200.

## Historical biogeography analyses

Historical biogeography analyses were carried out using the Dispersal-Extinction-Cladogenesis (DEC) model (Ree & Smith, 2008) as implemented in the R package BioGeoBEARS 1.1.2 (Matzke, 2014). As an input tree, the final total-evidence dated phylogeny estimated with MrBayes was used after pruning all species not belonging to Parnassiinae using Mesquite (Maddison & Maddison, 2021). Six areas were defined based on Brummit et al. (2001) and Condamine et al. (2018) as follows: [WP] Western Palearctic (Europe and Northern Africa), [WAC] Western Asia and Caucasus, [CA] Central Asia, [HTP] Himalaya, Tibetan Plateau, and Indian foothills, [EP] Eastern Palearctic (gathering of [SI] Eastern Russia, Siberia, [MO] Mongolian steppes, Altai Mountains, [CJ] Northern China, Korea and Japan in Condamine et al., 2018), and [WN] Western Nearctic. A maximum number of three areas for ancestral ranges was allowed to reflect the largest distribution range of the species (maximum three areas recorded for *Parnassius simo* Gray). To better consider paleogeographic changes through time, the five time slices and the adjacency matrix from Condamine et al. (2018) were used.

Two alternative models were tested in BioGeoBEARS: [M0] a null model with no constraints or time slices, and [M1] a reference model using the area adjacency matrix and the time slices from Condamine et al. (2018). For each of these models, four analyses were performed resulting in a total of eight analyses: [M0-1/M1-1] with the dated phylogeny from Condamine et al. (2018) based on node dating (birth-death model using molecular and morphological datasets, see Condamine et al., 2018 for details), [M0-2/M1-2] with the dated phylogeny from TED without the fossils, [M0-3/M1-3] with the dated phylogeny from TED assuming that fossil ranges are true (no-fossil areas coded as “0”), [M0-4/M1-4] with the dated phylogeny from TED using the positive constraint strategy (no-fossil areas coded as “?” instead of true absence, which means that the range of the fossils can be extended by the model beyond the

known area, see Coiro et al., 2023). A model ranking was done relying on the log likelihoods (LnL) of each model.

## RESULTS

The results and therefore the discussion are mostly focused on the backbone tree (i.e. intergeneric relationships), however the results on the intrageneric relationships are available in the Figshare repository. The main analyses and results are summarized on Figure 2.

### Bayesian total-evidence dating

### Effect of dating phylogenies on fossil placements and node supports

All analyses testing for the placement of the fossils resulted in well supported phylogenies with 72% to 77% of the nodes having PP ≥ 0.95 (Appendices S2, S3). †*Doritites bosniaskii* is recovered sister to *Archon* with high support (PP ≥ 0.95 in six of the eight analyses, Table 1), with a minimum support of PP = 0.88 (TEP with †*Doritites bosniaskii* subsampling). †*Thaites ruminiana* is either recovered sister to Luehdorfiini + Sericinini, †*Doritites bosniaskii* + *Archon*, or to Parnassiini with PP ranging from 0.456 to 0.84 (see Table 1). The two analyses containing the subsampled characters for †*Doritites bosniaskii* recovered †*Thaites ruminiana* sister to †*Doritites bosniaskii* + *Archon*. For the analyses that recovered †*Thaites ruminiana* as sister to Luehdorfiini + Sericinini, when looking at the RoguePlots †*Thaites ruminiana* is however recovered sister to Parnassiini with a support between 0.40 and 0.50 (see Appendix S3; Table 1). Overall, the long branch of *Baronia brevicornis* did not have any effect on the topology and branch support (from 72 to 73% of the nodes supported with or without *Baronia brevicornis*). Simultaneously inferring and dating the phylogeny led to higher supports for the two fossils compared to only inferring the tree (from an average of PP = 0.93 to 0.98 for †*Doritites bosniaskii*, and from an average of PP = 0.57 to 0.64 for †*Thaites ruminiana*; see Table 1). The subsampling of †*Doritites bosniaskii* characters did not increase the support of the fossil (see Table 1). However, in the case of †*Thaites ruminiana*, the subsampling of its characters led to higher support (from PP = 0.456 to 0.51 for phylogeny only, and from PP = 0.599 to 0.76 for simultaneous dating and phylogenetic inferences; see Table 1).

### Effect of clock partitioning and clock models on the node age estimates

The median crown age estimates for Parnassiinae, Parnassiini, Luehdorfiini and Sericinini did not vary widely across the different clock partition analyses (Fig. 3, 4; Appendices S2, S4). With the increase in the number of partitions, harmonic means substantially improved from -47,485.29 to -47,034.23. However, when many partitions are used (nine and twelve), runs did not converge for the clock models (ESS = 17 for *clockrate{all}* parameter with nine partitions and two runs; ESS = 4 for *clockrate{all}* parameter with 12 partitions and two runs; Appendix S4). Thus, three parti-
Effect of dating phylogenies on fossil placements and node supports

<table><caption>Figure 2. Workflow of the total-evidence phylogenetic analyses performed in this study.</caption><thead><tr><th colspan="4">Total-evidence phylogeny</th><th colspan="4">Total-evidence dating</th></tr><tr><th colspan="2">Outgroups</th><th colspan="2">Subsampling</th><th colspan="2">Outgroups</th><th colspan="2">Subsampling</th></tr></thead><tbody><tr><td>Papilioninae<br/>+ Baronia</td><td>Papilioninae</td><td>†Thaites<br/>ruminiana</td><td>†Doritites<br/>bosniaskii</td><td>Papilioninae<br/>+ Baronia</td><td>Papilioninae</td><td>†Thaites<br/>ruminiana</td><td>†Doritites<br/>bosniaskii</td></tr><tr><td colspan="8" style="text-align: center; font-size: 1.5em; line-height: 1.2em;">↓</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Effect of clock partitioning on node ages</td></tr><tr><td rowspan="2">1 clock<br/>All</td><td rowspan="2">2 clocks<br/>Molecular;<br/>Morphology</td><td rowspan="2">3 clocks<br/>Mitochondrial;<br/>Nuclear;<br/>Morphology</td><td rowspan="2">9 clocks<br/>1 by gene;<br/>1 for Morphology</td><td rowspan="2">12 clocks<br/>1 by gene;<br/>3 for Morphology</td></tr><tr><td></td></tr><tr><td colspan="8" style="text-align: center; font-size: 1.5em; line-height: 1.2em;">↓</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Effect of clock models on node ages</td></tr><tr><td rowspan="2">Mitochondrial</td><td rowspan="2"></td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Nuclear</td><td rowspan="2"></td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Morphology</td><td rowspan="2"></td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td><td>IGR</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td colspan="8" style="text-align: center; font-size: 1.5em; line-height: 1.2em;">↓</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Final total-evidence dating</td></tr><tr><td colspan="8" style="text-align: center; font-size: 1.5em; line-height: 1.2em;">↓</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Node dating</td></tr><tr><td colspan="8" style="text-align: center; font-size: 1.5em; line-height: 1.2em;">↓</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Historical biogeography analyses</td></tr><tr><td colspan="2">M0</td><td colspan="2">M1</td><td colspan="2">M0</td><td colspan="2">M1</td></tr><tr><td colspan="2">Condamine et al. 2018</td><td colspan="2">Condamine et al. 2018</td><td colspan="2">Condamine et al. 2018</td><td colspan="2">Condamine et al. 2018</td></tr><tr><td colspan="2">TED without fossils</td><td colspan="2">TED without fossils</td><td colspan="2">TED without fossils</td><td colspan="2">TED without fossils</td></tr><tr><td colspan="2">TED true absence</td><td colspan="2">TED true absence</td><td colspan="2">TED true absence</td><td colspan="2">TED true absence</td></tr><tr><td colspan="2">TED missing</td><td colspan="2">TED missing</td><td colspan="2">TED missing</td><td colspan="2">TED missing</td></tr></tbody></table>

<table><caption>Table 1. Fossil placements and associated posterior probabilities (PP). The eight models tested are represented, with 'TEP' corresponding to total-evidence phylogeny, 'TED' to total evidence dating, 'Papilioninae' when *Baronia brevicornis* was removed, '†Thaites ruminiana' when the subsampling matrix containing the characters of †Thaites ruminiana was used (81 morphological characters), and '†Doritites bosniaskii' when the subsampling matrix containing the characters of †Doritites bosniaskii was used (96 morphological characters). The RoguePlot placements and supports are represented for †Thaites ruminiana.</caption><thead><tr><th rowspan="2">Analysis</th><th colspan="2">†Doritites bosniaskii</th><th colspan="4">†Thaites ruminiana</th></tr><tr><th>Placement</th><th>Support (PP value)</th><th>Placement</th><th>Support (PP value)</th><th>RoguePlot placement</th><th>RoguePlot support</th></tr></thead><tbody><tr><td>(i) TEP all outgroups</td><td>sister Archon</td><td>0.95</td><td>sister Luehdorfini + Sericinini</td><td>0.456</td><td>sister Parnassiini</td><td>0.40-0.50</td></tr><tr><td>(ii) TEP Papilioninae</td><td>sister Archon</td><td>0.926</td><td>sister Luehdorfini + Sericinini</td><td>0.497</td><td>sister Parnassiini</td><td>0.40-0.50</td></tr><tr><td>(iii) TEP †Thaites ruminiana</td><td>sister Archon</td><td>0.97</td><td>sister Luehdorfini + Sericinini</td><td>0.51</td><td>sister Parnassiini</td><td>0.40-0.50</td></tr><tr><td>(iv) TEP †Doritites bosniaskii</td><td>sister Archon</td><td>0.88</td><td>sister Archon + †Doritites bosniaskii</td><td>0.84</td><td>sister Archon</td><td>0.80-0.90</td></tr><tr><td>(v) TED all outgroups</td><td>sister Archon</td><td>0.98</td><td>sister Parnassiini</td><td>0.599</td><td>sister Parnassiini</td><td>0.60-0.70</td></tr><tr><td>(vi) TED Papilioninae</td><td>sister Archon</td><td>0.978</td><td>sister Parnassiini</td><td>0.687</td><td>sister Parnassiini</td><td>0.70-0.80</td></tr><tr><td>(vii) TED †Thaites ruminiana</td><td>sister Archon</td><td>0.999</td><td>sister Parnassiini</td><td>0.76</td><td>sister Parnassiini</td><td>0.70-0.80</td></tr><tr><td>(viii) TED †Doritites bosniaskii</td><td>sister Archon</td><td>0.976</td><td>sister Archon + †Doritites bosniaskii</td><td>0.53</td><td>sister Archon</td><td>0.90-1.00</td></tr></tbody></table>

tions (mitochondrial, nuclear and morphology; ESS = 86 for *clockrate*[all] with two runs, and ESS = 379 with four runs) are used for the final TED analysis.

However, the median age estimates varied widely across the different clock model analyses, ranging from 51 to 74

Ma for Parnassiinae, with the youngest boundary of the 95% HPD from 33 Ma to 47 Ma, and the oldest boundary of the 95% HPD from 71 to 94 Ma (Fig. 3, 4; Appendices S2, S5). The TED analyses with the three best harmonic means were the one using IGR clock model for mitochondrial and
nuclear partitions, and either IGR, TK02 or Strict clock models for the morphological partition (harmonic means of -47,220.31, -47,236.49, and -47,253.25, respectively). The TED analysis with the best harmonic mean is the one using the IGR clock model for all three partitions, which is therefore used for the final TED analysis.

## Final total-evidence dated phylogeny

The runs of the final TED analysis with three partitions - each set with the IGR clock model - converged well (ESS = 379 for *clockrate*[all]; see Appendix S6) and resulted in a globally well supported dated phylogeny (78% of PP ≥ 0.95). In this analysis, †*Doritites bosniaskii* was recovered sister to *Archon* (PP = 0.99), and †*Thaites ruminiana* was recovered sister to *Parnassiini* (PP = 0.58). The final TED analysis estimated the age of *Parnassiinae* ca. 57 Ma in the late Paleocene (95% HPD = 37.19-75.67 Ma), the age of *Parnassiini* ca. 31 Ma in the early Oligocene (95% HPD = 19.80-43.27 Ma), the age of *Luehdorfiini* ca. 32 Ma in the early Oligocene (95% HPD = 20.51-45.21 Ma), and the age of *Sericinini* ca. 38 Ma in the late Eocene (95% HPD = 25.02-52.09 Ma) (Fig. 3, 4, red bar). The origins of the subfamily and tribes inferred with the final TED analysis were estimated at older ages than the node dating analyses, which recovered *Parnassiinae* ca. 38 Ma in the late Eocene (95% HPD = 26.17-46.25 Ma), *Parnassiini* ca. 21 Ma in the early Miocene (95% HPD = 14.40-27.11 Ma), *Luehdorfiini* ca. 22 Ma in the early Miocene (95% HPD = 14.58-28.88 Ma), and *Sericinini* ca. 27 Ma in the Oligocene (95% HPD = 17.53-32.91 Ma) (Fig. 4, dark brown bar). The final TED analysis recovered the origin of *Parnassius* during the early Miocene (ca. 19 Ma, 95% HPD = 12.84-27.05 Ma), while the origins of the other genera are inferred during the late Miocene-Pliocene (see Appendix S6). Overall, the divergence times between sister species are estimated during the late Pliocene.

## Historical biogeography analyses

Overall, the biogeographic analyses made with the M1 constraints were more likely than the analyses with the M0 constraints (see Table 2; Appendix S8). The four M1 analyses (whether or not including fossils) estimated a most likely ancestral origin for *Parnassiinae* in Western Palearctic + Western Asia and Caucasus (Table 2; Fig. 5; relative probability from 39% to 55%). The second-best inferred origin included only Western Palearctic (relative probability from 14% to 45%). Without fossils, there was a higher uncertainty for the ancestral range at the root (53% when summing the percentages of the two best ancestral ranges) compared to inferences made when fossils are included in the dating (maximum of 94%, when summing the percentages of the two best ancestral ranges). Overall, the ancestral geographic ranges for all tribes and genera of *Parnassiinae* were congruent between the four analyses (Table 2; Fig. 5) and were the same for the two analyses including both extinct and extant taxa. The analysis coding fossil geographic ranges with missing data instead of true absence considered the geographic ranges of the two fossils as Western

Palearctic + Western Asia and Caucasus instead of Western Palearctic only (Table 2; Fig. 5; Appendix S7). The former was kept as the best analysis to discuss the biogeographic history of *Parnassiinae* since it has a likelihood slightly better compared to the analysis coding fossil geographic ranges with true absence (LnL = -186.88 versus LnL = -189.39; Table 2).

When considering deep node geographic ranges inferred with the DEC model including the two fossils, the origin in Western Palearctic + Western Asia and Caucasus ca. 57 Ma was followed by *in situ* speciation within Western Palearctic, then by colonization of Central Asia ca. 31 Ma for the ancestor of *Parnassiini*, followed by the colonization of the Himalaya and Tibetan Plateau ca. 19 Ma by the ancestor of *Parnassius*. There were six independent colonization events of the Eastern Palearctic during the middle-late Miocene: three from Western Palearctic by the ancestors of *Bhutanitis Atkinsoni*, *Luehdorfia Crügeri*, and *Sericinus Westwoodi*, and three from Central Asia and the Himalaya and Tibetan Plateau by the ancestors of *Parnassius (Driopa)*, *P. (Parnassius)*, and *P. (Sachaia)* (see Appendix S8).

## DISCUSSION

Overall, this study shows that integrating (only) two fossils as tips to estimate total-evidence dated phylogeny and historical biogeography of *Parnassiinae* can offer a different view of their macroevolutionary history. Hereafter, the dating and biogeographic history of *Parnassiinae* are discussed and used to provide recommendations for future empirical TED studies.

### Bayesian total-evidence dating

### Effect of dating phylogenies on fossil placements and node supports

The exploration of fossil placements led to similar results than Condamine et al. (2018), with †*Doritites bosniaskii* recovered as sister species to *Archon*, and †*Thaites ruminiana* more probably sister to *Parnassiini*, but sometimes sister to *Luehdorfiini* + *Sericinini*. The results of total-evidence dating compared to total-evidence phylogeny are in line with the importance of the ages of fossils as an additional source of information for their placement in the phylogeny (Donoghue & Yang, 2016): node support for the two fossils globally increased with the dating, and †*Thaites ruminiana* is inferred sister to *Luehdorfiini* + *Sericinini* in total-evidence phylogeny but is sister to *Parnassiini* in total-evidence dating.

Concerning the fossils which are more difficult to place (†*Thaites ruminiana* in this case), this study demonstrates the importance of exploring the placement through Rogue-Plots and character subsampling. On the one hand, even if sometimes †*Thaites ruminiana* is inferred sister to *Luehdorfiini* + *Sericinini*, it is also inferred sister to *Parnassiini* with a pretty high probability when looking at RoguePlots (Table 1, Appendix S3). On the other hand, the use of subsampling allows to limit missing data for the fossil of interest and is
<figure><img src="image_4.png" /><figcaption>Figure 3. Dated phylogenies of Parnassiinae for the best analyses of (a) clock partitioning, (b) clock models, and (c) the final total-evidence phylogeny.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4. Highest posterior density (HPD) and median ages of Parnassiinae, Parnassiini, Sericinini, and Luehdorfiini most recent common ancestors. Age estimates depend on clock partitioning and model tests, final total evidence dating (TED), node dating and previous studies. The tests for clock partitioning are represented in dark green and ranked from less to more complex; the tests for clock model are represented in light green and ranked from the best to the worst harmonic mean; the final TED analysis is represented in red; the node dating in dark brown; and the results of previous studies in light brown.</figcaption></figure>

a good way to investigate the placement of fossils (Tejada
et al., 2024): for †*Thaites ruminiana*, the use of subsampling
reinforced its sister relationship to Parnassiini. These dif-
ferent lines of evidence converge towards †*Thaites rumini-
ana* being more likely sister to Parnassiini, confirming what
has been shown previously (Condamine et al., 2018). How-
ever, it is important to note that the subsampling should be
used to explore only the placement of the fossil of interest.
Indeed, when using character subsampling to investigate

the placement of †Doritites bosniaskii, the analyses recov-
ered †Doritites bosniaskii sister to Archon but with a weak
support and †Thaites ruminiana was also sister to Archon
with a very high support (Table 1). This could be due to a
high proportion of missing data and fossil attraction, espe-
cially as Archon was previously classified within Parnassiini
(see Nazari et al., 2007 for details), due to the similar mor-
phology of Archon, Hypermnestra Ménétriés and Parnassius
(Munroe, 1961; Omoto et al., 2004).
<table><caption>Table 2. Results of the historical biogeography analyses. The eight models tested are represented, with 'M0' corresponding to the null model, 'M1' to the reference model, 'Condamine et al.' when using the dated phylogeny from Condamine et al. (2018), 'without fossils' when fossils were removed after the total evidence dating (TED) analysis, 'true absence' when the fossil distribution is restricted to the area where it was found, and 'missing' the range of the fossils can be extended by the model beyond the known area. The regions used correspond to: [WP] Western Palearctic, [WAC] Western Asia and Caucasus, [CA] Central Asia, [HTP] Himalaya, Tibetan Plateau, and Indian foothills, and [EP] Eastern Palearctic.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Analysis</th><th rowspan="2">Log likelihood</th><th colspan="10">Ancestral range reconstructions of most common recent ancestors</th></tr><tr><th>Parnassiinae</th><th>Luehdorfiini</th><th>Zerynthiini</th><th>Parnassiini</th><th>Luehdorfia</th><th>Archon</th><th>Bhutanitis</th><th>Zerynthia</th><th>Allancastria</th><th>Hypermnestra</th><th>Parnassius</th></tr></thead><tbody><tr><td rowspan="4">null model</td><td>[M0-1] Condamine et al.</td><td>-189.69</td><td>CA + HTP + EP</td><td>WAC + EP</td><td>EP</td><td>CA</td><td>EP</td><td>WAC</td><td>EP</td><td>WP</td><td>WP + WAC</td><td>WAC + CA</td><td>CA + HTP + EP</td></tr><tr><td>[M0-2] TED without fossils</td><td>-187.72</td><td>CA + HTP + EP</td><td>WAC + EP</td><td>EP</td><td>CA + HTP + EP</td><td>EP</td><td>WAC</td><td>EP</td><td>WP</td><td>WP + WAC</td><td>WAC + CA</td><td>CA + HTP + EP</td></tr><tr><td>[M0-3] TED true absence</td><td>-194.63</td><td>WP + CA + EP</td><td>WP + WAC + EP</td><td>EP</td><td>CA + HTP + EP</td><td>EP</td><td>WAC</td><td>EP</td><td>WP</td><td>WP + WAC</td><td>WAC + CA</td><td>CA + HTP + EP</td></tr><tr><td>[M0-4] TED missing</td><td>-188.99</td><td>WP + CA + EP</td><td>EP</td><td>EP</td><td>CA</td><td>EP</td><td>WAC</td><td>EP</td><td>WP</td><td>WP + WAC</td><td>WAC + CA</td><td>CA + HTP + EP</td></tr><tr><td rowspan="4">reference model</td><td>[M1-1] Condamine et al.</td><td>-184.51</td><td>WP + WAC</td><td>WP + EP</td><td>WP</td><td>WAC + CA</td><td>EP</td><td>WAC</td><td>EP</td><td>WP</td><td>WP + WAC</td><td>WAC + CA</td><td>CA + HTP + EP</td></tr><tr><td>[M1-2] TED without fossils</td><td>-185.36</td><td>WP + WAC</td><td>WP + WAC</td><td>WP</td><td>WAC + CA</td><td>EP</td><td>WAC</td><td>EP</td><td>WP</td><td>WP + WAC</td><td>WAC + CA</td><td>CA + HTP</td></tr><tr><td>[M1-3] TED true absence</td><td>-189.39</td><td>WP + WAC</td><td>WP</td><td>WP</td><td>WAC + CA</td><td>EP</td><td>WAC</td><td>EP</td><td>WP</td><td>WP + WAC</td><td>WAC + CA</td><td>CA + HTP</td></tr><tr><td>[M1-4] TED missing</td><td>-186.88</td><td>WP + WAC</td><td>WP</td><td>WP</td><td>WAC + CA</td><td>EP</td><td>WAC</td><td>EP</td><td>WP</td><td>WP + WAC</td><td>WAC + CA</td><td>CA + HTP</td></tr></tbody></table>
## Effect of clock partitioning and clock models on the node age estimates

The question of what is the best way to partition the data is still open, and was even described as “more art than science” by Angelis et al. (2018). Here, increasing the number of partitions leads to better harmonic means but less and less run convergence between molecular clocks. The use of three partitions (mitochondrial, nuclear, and morphology) was the more complex partitioning strategy where clear differences between molecular clocks could be observed (see Appendix S4), probably explaining the statistical convergence of clock rates. However, we did not observe substantial variation in the ages between different analyses with different numbers of partitions. This is consistent with the observations by Angelis et al. (2018), who found that partitioning schemes have little effect on the posterior time estimates, as long as the prior assumptions are correct and the clock models not seriously violated.

The variations observed in terms of age depending on the clock model show the importance of testing for clock models. The IGR model is recovered as better fitting the data compared to TK02 and a strict clock which is consistent with what has been observed previously (Ronquist et al., 2012a and Zhang et al., 2016 for comparison with TK02; Gavryushkina et al., 2017 for comparison with strict clock). Using a relaxed clock compared to a strict clock, Gavryushkina et al. (2017) observed that age estimates slightly shifted toward the past with increased length of 95% HPD, however it is important to keep in mind that they found a younger age using TED for the crown of the penguin radiation. When comparing the results of the different analyses, a trend toward older age estimates with TK02 is observed (Fig. 1). It has already been noted by Zhang et al. (2016) that TK02 had difficulties to fit well the data since autocorrelated evolutionary rates are not expected to shift rapidly, and thus do not fit well datasets with rapid shifts in evolutionary rates, as it has been identified in Parnassiinae (Condamine et al., 2018). Furthermore, Ronquist et al. (2012a) showed on Hymenoptera that the IGR model fits the fossil record better than the TK02 model due to the extension of the deep branches resulting in the autocorrelated models. Therefore, our study is in line with previous research stating that the IGR model seems to be a better model when using the TED due to the integration of the fossil record. However, this study also points out that the estimation of the clock model for the morphological partition is the most complicated (the three best models recovered either IGR, TK02 and strict clock models for the morphological partition) and should be particularly explored, highlighting the need for robust models capable of comparing analyses including fossils and morphological data.

## Final total-evidence dated phylogeny

The origin of Parnassiinae inferred by TED analyses is older than previously thought by Condamine et al. (2018) with a larger length of 95% HPD: *ca.* 57 Ma (95% HPD = 37.19-75.67 Ma) compared to *ca.* 38.6 Ma (95% HPD =

31.4-46.6 Ma). The broader 95% HPD with TED analyses than with node dating analyses is consistent with what was previously observed (Arcila et al., 2015; O'Reilly et al., 2015; but see Ronquist et al., 2012a). This is probably due to the incorporation of morphological characters, which have been identified as leading to higher length of 95% HPD (Barba-Montoya et al., 2021). The use of large morphological matrices could be a solution to improve the precision of divergence time estimates while taking into account the quality of morphological characters (Barido-Sottani et al., 2020; Luo et al., 2020), as missing data decreases topological accuracy and precision (*sensu* Mongiardino Koch et al., 2021). The older age is expected as it has been observed in an array of groups while using TED (e.g. Arcila et al., 2015; Casali et al., 2020; O'Reilly et al., 2015; Wood et al., 2013). However, each example is different and has to be examined again before drawing any conclusions (Ronquist et al., 2016). Indeed, compared to node dating, TED analyses recovered younger ages (Gavryushkina et al., 2017), convergent ages (Parks et al., 2022), older ages (e.g. Arcila et al., 2015; Casali et al., 2020; Heritage & Seiffert, 2022; Matschiner et al., 2017; Near & Kim, 2021; Wood et al., 2013), as well as a mix of older and younger ages (Ronquist et al., 2012a).

The older ages inferred for Parnassiinae are probably not artificial as a moderate size morphological matrix of high-quality characters (*ca.* 180 characters, with at least 96 characters for the fossils) and a diversified sampling were used, which are considered as two parameters that can lead to increased deep root attraction if they are not taken into account (Ronquist et al., 2016). Furthermore, the use of the FBD process usually produces estimates of divergence time that are realistic (Truman et al., 2024), even when there is a biased sampling of fossils and extant species (Heath et al., 2014). In addition, a single node calibration at the root with tip-dates as it was done here has already proved to be an effective method for obtaining non-artificial ages (Casali et al., 2020) and mixing node and tip dating is highly encouraged (O'Reilly & Donoghue, 2016). Finally, when comparing the age of Parnassiinae inferred by TED with other recent node dating analyses, the ages inferred by TED are not that older and the length of 95% HPD are similar or even smaller: *ca.* 53.39 Ma (95% HPD = 30.51-91.87 Ma) by Allio et al. (2021), and *ca.* 45.97 Ma (95% HPD = 27.60-82.64 Ma) by Tian et al. (2023) (see Fig. 1).

## Historical biogeography analyses

### A historical biogeographic scenario for Parnassiinae

The historical biogeographic analyses including the two fossils using positive constraint strategy recovered an origin in Western Palearctic + Western Asia and Caucasus *ca.* 57 Ma (Fig. 4). On the one hand, ancestors of Luehdorfiini and Sericinini experienced *in situ* speciation within Western Palearctic region followed by three independent colonizations of Eastern Palearctic (genera *Bhutanitis*, *Luehdorfia* and *Sericinus*) that occurred during or after the Oligocene (Fig. 5), corresponding to the drying up of the Turgai Strait
<figure><img src="image_6.png" /></figure>

Figure 5. Comparison of time-calibrated phylogenies and historical biogeography for Parnassiinae without and with fossils. The time-calibrated phylogeny based on node dating associated with historical biogeography by Condamine et al. (2018) is shown at the top, while the time-calibrated phylogeny based on TED (M1-4) associated with historical biogeography including the two fossils is shown at the bottom. The colored squares on nodes indicate the most likely ancestral distribution. Unsupported nodes (P<0.95) are represented by '*' for both-time calibrated phylogenies. Images of Parnassiinae shown on the right side have been made by Ariane Chotard.
(ca. 30 Ma; Sanmartín et al., 2001). On the other hand, our analysis inferred the colonization of Central Asia ca. 31 Ma for the ancestor of Parnassiini. Then, the colonization of the Himalaya and Tibetan Plateau ca. 19 Ma by the ancestor of *Parnassius* follows the onset of the Qinghai-Tibetan plateau during the Eocene (Favre et al., 2015), and corresponds to the intensification of the Himalaya and Tibetan plateau orogeny during the early Miocene (Favre et al., 2015; Tapponnier et al., 2001; Valdiya, 1999; Wang et al., 2008). This is congruent with what have been called the 'HTP origin' hypothesis as found by previous studies (Allio et al., 2021; Condamine et al., 2018; Nazari et al., 2007; Tian et al., 2023; Zhao et al., 2022).

The fact that our biogeographic analysis inferred that *Parnassius* is the only lineage that colonized the Himalaya and Tibetan Plateau could explain their great diversity compared to the other genera (*Bhutanitis* being in Northern India but not in the high mountains). Indeed, the dispersal into a new region with many reliefs could explain their great diversity as mountains can be at the same time 'museum' and 'cradle' of diversity (Rahbek et al., 2019). A hypothesis could be that the adaptation of *Parnassius* to altitudinal habitats probably spurred their diversification as shown by diversification analyses indicating a higher speciation rate for mountain species compared to lowland species (Condamine et al., 2018) due to allopatric speciation by facing divergent selection pressures ('cradle', Rahbek et al., 2019). Adaptation to mountains could also have limited the extinction of Parnassiinae during strong climatic changes as the Middle Miocene Climatic Optimum (17-15 Ma; Zachos et al., 2008) as mountains are known to act as buffers or altitudinal refuges ('museum', Rahbek et al., 2019). Furthermore, Tian et al. (2023) identified *Corydalis* as the most likely ancestral host plant of *Parnassius*. The diversity of *Corydalis* is concentrated in China and Himalaya and is estimated to have originated in the Eocene (ca. 37 Ma by Pérez-Gutiérrez et al., 2015; and ca. 49 Ma by Xu et al., 2022). According to Tian et al. (2023), the host-plant shift of the ancestor of *Parnassius* on *Corydalis* could have permitted the colonization of the Himalaya and Tibetan Plateau, supporting the hypothesis that only the genus *Parnassius* colonized the HTP.

## Effect of two fossils on biogeographic history of Parnassiinae

With node dating, an origin of Parnassiinae ca. 39 Ma in Central Asia followed by colonization of the other biogeographic regions was inferred (Condamine et al., 2018). By comparing the historical biogeography analysis including the two fossils with what have been found by Condamine et al. (2018), many differences are found in the backbone of the phylogeny. It is important to note that there is an effect of the number of areas in the model. Indeed, even when the fossils are not incorporated, the historical biogeography analyses integrating the two fossils found an origin in Western Palearctic + Western Asia and Caucasus just by reducing the number of areas. However, the probability for ancestral ranges in Western Palearctic + Western Asia and Caucasus increases when the fossils are included, reinforc-

ing this hypothesis. The scenario inferred by Condamine et al. (2018) appears as less likely that the one inferred by the historical biography analysis including the two fossils (see above) mostly because it recovered that Luehdorfiini and Sericinini diversified in broad geographic ranges (India + Turkey-Caucasus + Central Asia, and India + Turkey-Caucasus + Central Asia + HTP, respectively; see Fig. 5), which do not correspond to the observed ecology of Parnassiinae, as most of the species are highly endemic. These findings thus join previous studies indicating that a few fossils are important to consider when inferring the biogeographic history of a group, especially when the past distribution is different from today's distribution (Parks et al., 2022; Rule et al., 2020; Santos et al., 2022).

## CONCLUSIONS

This study confirms the previous placement of the two fossils in the dated phylogeny with †*Doritites bosniaskii* recovered as sister species to *Archon*, and †*Thaites ruminiana* more probably sister to Parnassiini. Including the two fossils allowed obtaining a more plausible dating and historical biogeography of the group, with an origin ca. 57 Ma in the late Paleocene in Western Palearctic + Western Asia and Caucasus, followed by colonization of Central Asia and Himalaya. The genus *Parnassius* is now inferred as the only genus that would have colonized the Himalaya and Tibetan Plateau after undergoing a change of host plant, which may explain its great extant diversity. For future studies on the macroevolution of Parnassiinae, we recommend the use of TED under FBD to perform dating analyses. If not possible, as for example while using genomic data, this dated phylogeny could serve as a basis for secondary calibrations, or large uniform prior with fossil calibrations should be used. This study adds to the growing body of TED studies showing the importance of exploring the placement of fossils for molecular dating when a high-quality morphological matrix (low amount of missing data) is available (e.g. Coiro et al., 2023; Spasojevic et al., 2021; Tejada et al., 2024). The placement of the fossils should be investigated using different clock models, RoguePlots, and character subsampling, particularly when fossils are difficult to place. This study highlights the importance of a trade-off between the number of partitions and the effect of clock models, as well as the need for robust models capable of comparing analyses including fossils and morphological data. Finally, we argue that a few fossils are enough to assess and question the origin and biogeographic history of a group, especially when they were distributed in different regions than the present center of diversity. In those cases, authors should use TED if possible, but also other ways to integrate fossil information in historical biogeography (e.g. constrain deep nodes with fossil distributions), or discuss the historical biogeography of the group in the light of what is known about fossil placement and distributions.

<figure><img src="image_7.png" /></figure>
**Funding**

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (project GAIA, agreement no. 851188). This project benefited from the Montpellier Bioinformatics Biodiversity platform supported by the LabEx CeMEB, an ANR "Investissements d'avenir" program (ANR-10-LABX-04-01).

**Acknowledgements**

We would like to thank the editor Bryan C. Carstens and the anonymous reviewer for constructive comments and suggestions. We thank Benoit Nabholz for comments on the manuscript and we are grateful to Ariane Chotard for providing pictures of Parnassiinae taken at the Lepidoptera collection from the National Museum of Natural History of Paris.

**Supporting Information**

Appendices containing morphological matrices, MrBayes nexus files for all analyses performed, and additional results for the dating and biogeography are presented in the online repository Figshare 10.6084/m9 figshare.25755504.

**Appendix S1.** Morphological characters.

**Appendix S2.** Table of dating analysis results.

**Appendix S3.** Fossil placement and node support analyses.

**Appendix S4.** Clock partitioning analyses.

**Appendix S5.** Clock model analyses.

**Appendix S6.** Final total-evidence dating analysis.

**Appendix S7.** Node dating analysis.

**Appendix S8.** Historical biogeography analysis results.

Submitted: May 13, 2024 EST, Accepted: September 30, 2024

EST
References

Allio, R., Nabholz, B., Wanke, S., Chomicki, G., Pérez-Escobar, O. A., Cotton, A. M., Clamens, A. L., Kergoat, G. J., Sperling, F. A. H., & Condamine, F. L. (2021). Genome-wide macroevolutionary signatures of key innovations in butterflies colonizing new host plants. *Nature Communications*, 12(1), 354. https://doi.org/10.1038/s41467-020-20507-3

Allio, R., Scornavacca, C., Nabholz, B., Clamens, A. L., Sperling, F. A. H., & Condamine, F. L. (2020). Whole genome shotgun phylogenomics resolves the pattern and timing of swallowtail butterfly evolution. *Systematic Biology*, 69(1), 38–60. https://doi.org/10.1093/sysbio/syz030

Angelis, K., Álvarez-Carretero, S., Dos Reis, M., & Yang, Z. (2018). An evaluation of different partitioning strategies for Bayesian estimation of species divergence times. *Systematic Biology*, 67(1), 61–77. https://doi.org/10.1093/sysbio/syx061

Arcila, D., Pyron, R. A., Tyler, J. C., Ortí, G., & Betancur-R, R. (2015). An evaluation of fossil tip-dating versus node-age calibrations in tetraodontiform fishes (Teleostei: Percomorphaceae). *Molecular Phylogenetics and Evolution*, 82, 131–145. https://doi.org/10.1016/j.ympev.2014.10.011

Bacon, C. D., Silvestro, D., Hoorn, C., Bogotá-Ángel, G., Antonelli, A., & Chazot, N. (2022). The origin of modern patterns of continental diversity in Mauritiinae palms: the Neotropical museum and the Afrotropical graveyard. *Biology Letters*, 18(11), 20220214. https://doi.org/10.1098/rsbl.2022.0214

Barba-Montoya, J., Tao, Q., & Kumar, S. (2021). Molecular and morphological clocks for estimating evolutionary divergence times. *BMC Ecology and Evolution*, 21(1), 83. https://doi.org/10.1186/s12862-021-01798-6

Barido-Sottani, J., Aguirre-Fernández, G., Hopkins, M. J., Stadler, T., & Warnock, R. C. (2019). Ignoring stratigraphic age uncertainty leads to erroneous estimates of species divergence times under the fossilized birth-death process. *Proceedings of the Royal Society B*, 286(1902), 20190685. https://doi.org/10.1098/rspb.2019.0685

Barido-Sottani, J., Van Tiel, N. M., Hopkins, M. J., Wright, D. F., Stadler, T., & Warnock, R. C. (2020). Ignoring fossil age uncertainty leads to inaccurate topology and divergence time estimates in time calibrated tree inference. *Frontiers in Ecology and Evolution*, 8, 183. https://doi.org/10.3389/fevo.2020.00183

Borowiec, M. L. (2016). AMAS: a fast tool for alignment manipulation and computing of summary statistics. PeerJ, 4, e1660. https://doi.org/10.7717/peerj.1660

Brummitt, R. K., Pando, F., Hollis, S., & Brummitt, N. A. (2001). *World geographical scheme for recording plant distributions* (Vol. 951, p. 952). International working group on taxonomic databases for plant sciences (TDWG).

Carpenter, F. M. (1992). *Treatise on invertebrate paleontology, Part R, Arthropoda 3-4*. Geological Society of America.

Casali, D. D. M., Dos Santos Júnior, J. E., Miranda, F. R., Santos, F. R., & Perini, F. A. (2020). Total-evidence phylogeny and divergence times of Vermilingua (Mammalia: Pilosa). *Systematics and Biodiversity*, 18(3), 216–227. https://doi.org/10.1080/14772000.2020.1729894

Coiro, M., Allio, R., Mazet, N., Seyfullah, L. J., & Condamine, F. L. (2023). Reconciling fossils with phylogenies reveals the origin and macroevolutionary processes explaining the global cycad biodiversity. *New Phytologist*, 240, 1616–1635. https://doi.org/10.1111/nph.19010

Condamine, F. L., Rolland, J., Höhna, S., Sperling, F. A. H., & Sanmartín, I. (2018). Testing the role of the Red Queen and Court Jester as drivers of the macroevolution of Apollo butterflies. *Systematic Biology*, 67(6), 940–964. https://doi.org/10.1093/sysbio/syy009

Darlim, G., Lee, M. S., Walter, J., & Rabi, M. (2022). The impact of molecular data on the phylogenetic position of the putative oldest crown crocodilian and the age of the clade. *Biology Letters*, 18(2), 20210603. https://doi.org/10.1098/rsbl.2021.0603

de Jong, R. (2017). Fossil butterflies, calibration points and the molecular clock (Lepidoptera: Papilionoidea). *Zootaxa*, 4270(1), 1–63. https://doi.org/10.11646/ZOOTAXA.4270.1.1

Donoghue, P. C., & Yang, Z. (2016). The evolution of methods for establishing evolutionary timescales. *Philosophical Transactions of the Royal Society B*, 371(1699), 20160020. https://doi.org/10.1098/rstb.2016.0020

Durden, C. J., & Rose, H. (1978). Butterflies from the middle Eocene: the earliest occurrence of fossil Papilionidae. *Pearce-Sellards Series of the Texas Memorial Museum*, 29, 1–25.
Favre, A., Päckert, M., Pauls, S. U., Jähnig, S. C., Uhl, D., Michalak, I., & Muellner-Riehl, A. N. (2015). The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. *Biological Reviews, 90*(1), 236–253. https://doi.org/10.1111/brv.12107

Gavryushkina, A., Heath, T. A., Ksepka, D. T., Stadler, T., Welch, D., & Drummond, A. J. (2017). Bayesian total-evidence dating reveals the recent crown radiation of penguins. *Systematic Biology*, *66*(1), 57–73. https://doi.org/10.1093/sysbio/syw060

Heath, T. A., Huelsenbeck, J. P., & Stadler, T. (2014). The fossilized birth–death process for coherent calibration of divergence-time estimates. *Proceedings of the National Academy of Sciences of the USA*, *111*(29), E2957–E2966. https://doi.org/10.1073/pnas.1319091111

Heritage, S., & Seiffert, E. R. (2022). Total evidence time-scaled phylogenetic and biogeographic models for the evolution of sea cows (Sirenia, Afrotheria). *PeerJ*, *10*, e13886. https://doi.org/10.7717/peerj.13886

Ho, S. Y., & Phillips, M. J. (2009). Accounting for calibration uncertainty in phylogenetic estimation of evolutionary divergence times. *Systematic Biology*, *58*(3), 367–380. https://doi.org/10.1093/sysbio/syp035

Huelsenbeck, J. P., Larget, B., & Alfaro, M. E. (2004). Bayesian phylogenetic model selection using reversible jump Markov chain Monte Carlo. *Molecular Biology and Evolution*, *21*(6), 1123–1133. https://doi.org/10.1093/molbev/msh123

Hunt, G., & Slater, G. (2016). Integrating paleontological and phylogenetic approaches to macroevolution. *Annual Review of Ecology, Evolution, and Systematics*, *47*, 189–213. https://doi.org/10.1146/annurev-ecolsys-112414-054207

Jiangzuo, Q., & Flynn, J. J. (2020). The earliest Ursine bear demonstrates the origin of plant-dominated omnivory in Carnivora. *iScience*, *23*(6), 101235. https://doi.org/10.1016/j.isci.2020.101235

Klopfstein, S., & Spasojevic, T. (2019). Illustrating phylogenetic placement of fossils using RoguePlots: an example from ichneumonid parasitoid wasps (Hymenoptera, Ichneumonidae) and an extensive morphological matrix. *PLoS One*, *14*(4), e0212942. https://doi.org/10.1371/journal.pone.0212942

Lanfear, R., Calcott, B., Ho, S. Y., & Guindon, S. (2012). PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. *Molecular Biology and Evolution*, *29*(6), 1695–1701. https://doi.org/10.1093/molbev/mss020

Lee, M. S., & Yates, A. M. (2018). Tip-dating and homoplasy: reconciling the shallow molecular divergences of modern gharials with their long fossil record. *Proceedings of the Royal Society B*, *285*(1881), 20181071. https://doi.org/10.1098/rspb.2018.1071

Lewis, P. O. (2001). A likelihood approach to estimating phylogeny from discrete morphological character data. *Systematic Biology*, *50*(6), 913–925. https://doi.org/10.1080/106351501753462876

Luo, A., Duchêne, D. A., Zhang, C., Zhu, C. D., & Ho, S. Y. (2020). A simulation-based evaluation of tip-dating under the fossilized birth–death process. *Systematic Biology*, *69*(2), 325–344. https://doi.org/10.1093/sysbio/syz038

Maddison, W., & Maddison, D. (2021). *Mesquite 3. A Modular System for Evolutionary Analysis*. http://www.mesquiteproject.org

Matschiner, M. (2019). Selective sampling of species and fossils influences age estimates under the fossilized birth–death model. *Frontiers in Genetics*, *10*, 477091. https://doi.org/10.3389/fgene.2019.01064

Matschiner, M., Musilová, Z., Barth, J. M., Starostová, Z., Salzburger, W., Steel, M., & Bouckaert, R. (2017). Bayesian phylogenetic estimation of clade ages supports trans-Atlantic dispersal of cichlid fishes. *Systematic Biology*, *66*(1), 3–22. https://doi.org/10.1093/sysbio/syw076

Matzke, N. J. (2014). Model selection in historical biogeography reveals that founder-event speciation is a crucial process in island clades. *Systematic Biology*, *63*(6), 951–970. https://doi.org/10.1093/sysbio/syu056

Mongiardino Koch, N., Garwood, R. J., & Parry, L. A. (2021). Fossils improve phylogenetic analyses of morphological characters. *Proceedings of the Royal Society B*, *288*(1950), 20210044. https://doi.org/10.1098/rspb.2021.0044

Mongiardino Koch, N., & Parry, L. A. (2020). Death is on our side: paleontological data drastically modify phylogenetic hypotheses. *Systematic Biology*, *69*(6), 1052–1067. https://doi.org/10.1093/sysbio/syaa023

Monroe, M. J., & Bokma, F. (2010). Little evidence for Cope's rule from Bayesian phylogenetic analysis of extant mammals. *Journal of Evolutionary Biology*, *23*(9), 2017–2021. https://doi.org/10.1111/j.1420-9101.2010.02051.x

Munroe, E. (1961). The classification of the Papilionidae (Lepidoptera). *The Canadian Entomologist*, *17*, 1–51. https://doi.org/10.4039/entm9217fv
Nazari, V., Zakharov, E. V., & Sperling, F. A. H. (2007). Phylogeny, historical biogeography, and taxonomic ranking of Parnassiinae (Lepidoptera, Papilionidae) based on morphology and seven genes. *Molecular Phylogenetics and Evolution, 42*(1), 131–156.
https://doi.org/10.1016/j.ympev.2006.06.022

Near, T. J., & Kim, D. (2021). Phylogeny and time scale of diversification in the fossil-rich sunfishes and black basses (Teleostei: Percomorpha: Centrarchidae). *Molecular Phylogenetics and Evolution, 161*, 107156. https://doi.org/10.1016/j.ympev.2021.107156

Omoto, K., Katoh, T., Chichvarkin, A., & Yagi, T. (2004). Molecular evolution of the 'Apollo' butterflies of the genus Parnassius (Lepidoptera: Papilionidae). *Gene, 326*, 141–147. https://doi.org/10.1016/j.gene.2003.10.020

O'Reilly, J. E., & Donoghue, P. C. (2016). Tips and nodes are complementary not competing approaches to the calibration of molecular clocks. *Biology Letters, 12*(4), 20150975. https://doi.org/10.1098/rsbl.2015.0975

O'Reilly, J. E., & Donoghue, P. C. (2021). Fossilization processes have little impact on tip-calibrated divergence time analyses. *Palaeontology, 64*(5), 687–697. https://doi.org/10.1111/pala.12564

O'Reilly, J. E., Dos Reis, M., & Donoghue, P. C. (2015). Dating tips for divergence-time estimation. *Trends in Genetics, 31*(11), 637–650. https://doi.org/10.1016/j.tig.2015.08.001

Parks, R., Harrington, S. M., & Thomson, R. C. (2022). Divergence dating and biogeography of Xenosauridae including fossils as terminal taxa. *Journal of Herpetology, 56*(3), 349–354.
https://doi.org/10.1670/21-068

Pérez-Gutiérrez, M. A., Romero-García, A. T., Fernández, M. C., Blanca, G., Salinas-Bonillo, M. J., & Suárez-Santiago, V. N. (2015). Evolutionary history of fumitories (subfamily Fumarioideae, Papaveraceae): An old story shaped by the main geological and climatic events in the Northern Hemisphere. *Molecular Phylogenetics and Evolution, 88*, 75–92.
https://doi.org/10.1016/j.ympev.2015.03.026

Pett, W., & Heath, T. A. (2020). Inferring the timescale of phylogenetic trees from fossil data. In C. Scornavaca, F. Delsuc, & N. Galtier (Eds.), *Phylogenetics in the Genomic Era*.

Pyron, R. A. (2011). Divergence time estimation using fossils as terminal taxa and the origins of Lissamphibia. *Systematic Biology, 60*(4), 466–481.
https://doi.org/10.1093/sysbio/syr047

Rahbek, C., Borregaard, M. K., Antonelli, A., Colwell, R. K., Holt, B. G., Nogues-Bravo, D., Rasmussen, C. M. Ø., Richardson, K., Rosing, M. T., Whittaker, R. J., & Fjeldså, J. (2019). Building mountain biodiversity: Geological and evolutionary processes. *Science, 365*(6458), 1114–1119. https://doi.org/10.1126/science.aax0151

Rambaut, A., Drummond, A. J., Xie, D., Baele, G., & Suchard, M. A. (2018). Posterior summarization in Bayesian phylogenetics using Tracer 1.7. *Systematic Biology, 67*(5), 901–904. https://doi.org/10.1093/sysbio/syy032

Rasnitsyn, A. P., & Zherikhin, V. V. (2002). Appendix: alphabetic list of selected insect fossil sites. Impression fossils. In A. P. Rasnitsyn & D. L. J. Quicke (Eds.), *History of insects* (pp. 437–444). Kluwer Academic Publishers.

Rebel, H. (1898). Fossile Lepidopteren aus der Miocänformation von Gabbro. *Sitzungsberichte der Kaiserlichen Akademie der Wissenschaften. Mathematisch-Naturwissenschaftliche Classe, 107*, 731–745.

Ree, R. H., & Smith, S. A. (2008). Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. *Systematic Biology, 57*(1), 4–14. https://doi.org/10.1080/10635150701883881

Renner, S. S., Grimm, G. W., Kapli, P., & Denk, T. (2016). Species relationships and divergence times in beeches: New insights from the inclusion of 53 young and old fossils in a birth–death clock model. *Philosophical Transactions of the Royal Society B, 371*(1699), 20150135. https://doi.org/10.1098/rstb.2015.0135

Rieux, A., & Balloux, F. (2016). Inferences from tip-calibrated phylogenies: a review and a practical guide. *Molecular Ecology, 25*(9), 1911–1924.
https://doi.org/10.1111/mec.13586

Ronquist, F., Lartillot, N., & Phillips, M. J. (2016). Closing the gap between rocks and clocks using total-evidence dating. *Philosophical Transactions of the Royal Society B, 371*(1699), 20150136. https://doi.org/10.1098/rstb.2015.0136

Rule, J. P., Adams, J. W., Marx, F. G., Evans, A. R., Tennyson, A. J., Scofield, R. P., & Fitzgerald, E. M. (2020). First monk seal from the Southern Hemisphere rewrites the evolutionary history of true seals. *Proceedings of the Royal Society B, 287*(1938), 20202318. https://doi.org/10.1098/rspb.2020.2318
Sanmartín, I., Enghoff, H., & Ronquist, F. (2001). Patterns of animal dispersal, vicariance and diversification in the Holarctic. *Biological Journal of the Linnean Society*, 73(4), 345–390. https://doi.org/10.1111/j.1095-8312.2001.tb01368.x

Santos, B. F., Sandoval, M., Spasojevic, T., Giannotta, M. M., & Brady, S. G. (2022). A parasitoid puzzle: Phylogenomics, Total-evidence dating, and the role of Gondwanan vicariance in the diversification of Labeninae (Hymenoptera, Ichneumonidae). *Insect Systematics and Diversity*, 6(5), 3. https://doi.org/10.1093/isd/ixac015

Scudder, S. H. (1875). Fossil butterflies. *Memoirs of the American Association for the Advancement of Science*, 1, 1–99. https://doi.org/10.5962/bhl.title.964

Selvatti, A. P., Moreira, F. R. R., Carvalho, D. C., Prosdocimi, F., de Moraes Russo, C. A., & Junqueira, A. C. M. (2023). Phylogenomics reconciles molecular data with the rich fossil record on the origin of living turtles. *Molecular Phylogenetics and Evolution*, 183, 107773. https://doi.org/10.1016/j.ympev.2023.107773

Smith, M. E., Singer, B., & Carroll, A. (2003). 40Ar/39Ar geochronology of the Eocene Green River Formation, Wyoming. *Geological Society of America Bulletin*, 115(5), 549–565. https://doi.org/10.1130/0016-7606(2003)115

Sohn, J. C., Labandeira, C. C., Davis, D. R., & Mitter, C. (2012). An annotated catalog of fossil and subfossil Lepidoptera (Insecta: Holometabola) of the world. *Zootaxa*, 3286, 1–132. https://doi.org/10.11646/zootaxa.3286.1.1

Spasojevic, T., Broad, G. R., Sääksjärvi, I. E., Schwarz, M., Ito, M., Korenko, S., & Klopfstein, S. (2021). Mind the outgroup and bare branches in total-evidence dating: a case study of pimpliform Darwin wasps (Hymenoptera, Ichneumonidae). *Systematic Biology*, 70(2), 322–339. https://doi.org/10.1093/sysbio/syaa079

Tapponnier, P., Zhiqin, X., Roger, F., Meyer, B., Arnaud, N., Wittlinger, G., & Jingsui, Y. (2001). Oblique stepwise rise and growth of the Tibet Plateau. *Science*, 294(5547), 1671–1677. https://doi.org/10.1126/science.105978

Tejada, J. V., Antoine, P. O., Münch, P., Billet, G., Hautier, L., Delsuc, F., & Condamine, F. L. (2024). Bayesian total-evidence dating revisits sloth phylogeny and biogeography: a cautionary tale on morphological clock analyses. *Systematic Biology*, syad069. https://doi.org/10.1093/sysbio/syad069

Tian, X., Mo, S., Liang, D., Wang, H., & Zhang, P. (2023). Amplicon capture phylogenomics provides new insights into the phylogeny and evolution of alpine Parnassius butterflies (Lepidoptera: Papilionidae). *Systematic Entomology*, 48(4), 571–584. https://doi.org/10.1111/syen.12591

Truman, K., Vaughan, T. G., Gavryushkin, A., & Gavryushkina, A. (2024). The fossilised birth-death model is identifiable. *bioRxiv*. https://doi.org/10.1101/2024.02.08.579547

Valdiya, K. S. (1999). Rising Himalaya: advent and intensification of monsoon. *Current Science*, 76(4), 514.

Venditti, C., Meade, A., & Pagel, M. (2011). Multiple routes to mammalian diversity. *Nature*, 479(7373), 393–396. https://doi.org/10.1038/nature10516

Wang, C., Zhao, X., Liu, Z., Lippert, P. C., Graham, S. A., Coe, R. S., Yi, H., Zhu, L., Liu, S., & Li, Y. (2008). Constraints on the early uplift history of the Tibetan Plateau. *Proceedings of the National Academy of Sciences of the USA*, 105(13), 4987–4992. https://doi.org/10.1073/pnas.0703595105

Weiss, J. C. (1991–2005). *The Parnassiinae of the World. Part I 1991, II 1992, III 1999, IV 2005*. *Sciences Nat.*

Wisniewski, A. L., Lloyd, G. T., & Slater, G. J. (2022). Extant species fail to estimate ancestral geographical ranges at older nodes in primate phylogeny. *Proceedings of the Royal Society B*, 289(1975), 20212535. https://doi.org/10.1098/rspb.2021.2535

Wood, H. M., Matzke, N. J., Gillespie, R. G., & Griswold, C. E. (2013). Treating fossils as terminal taxa in divergence time estimation reveals ancient vicariance patterns in the palpimanoid spiders. *Systematic Biology*, 62(2), 264–284. https://doi.org/10.1093/sysbio/sys092

Xie, W., Lewis, P. O., Fan, Y., Kuo, L., & Chen, M. H. (2011). Improving marginal likelihood estimation for Bayesian phylogenetic model selection. *Systematic Biology*, 60(2), 150–160. https://doi.org/10.1093/sysbio/syq085

Xu, X., Li, X., & Wang, D. (2022). New insights into the backbone phylogeny and character evolution of Corydalis (Papaveraceae) based on plastome data. *Frontiers in Plant Science*, 13, 926574. https://doi.org/10.3389/fpls.2022.926574

Yang, Z., & Rannala, B. (2006). Bayesian estimation of species divergence times under a molecular clock using multiple fossil calibrations with soft bounds. *Molecular Biology and Evolution*, 23(1), 212–226. https://doi.org/10.1093/molbev/msj024
Zachos, J. C., Dickens, G. R., & Zeebe, R. E. (2008). An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics. *Nature*, **451**(7176), 279–283. https://doi.org/10.1038/nature06588

Zhang, C., Stadler, T., Klopfstein, S., Heath, T. A., & Ronquist, F. (2016). Total-evidence dating under the fossilized birth–death process. *Systematic Biology*, **65**(2), 228–249. https://doi.org/10.1093/sysbio/syy080

Zhao, Y., He, B., Tao, R., Su, C., Ma, J., Hao, J., & Yang, Q. (2022). Phylogeny and biogeographic history of Parnassius butterflies (Papilionidae: Parnassiinae) reveal their origin and deep diversification in West China. *Insects*, **13**(5), 406. https://doi.org/10.3390/insects13050406