Dating the origin and diversification of Pan-Chelidae (Testudines, Pleurodira) under multiple molecular clock approaches

Pan-Chelidae (Testudines, Pleurodira) is a group of side-necked turtles with a currently disjointed distribution in South America and Australasia and characterized by two morphotypes: the long-necked and the short-necked chelids. Both geographic groups include both morphotypes, but different phylogenetic signals are obtained from morphological and molecular data, suggesting the monophyly of the long-necked chelids or the independent evolution of this trait in both groups. In this paper, we addressed this conflict by compiling and editing available molecular and morphological data for Pan-Chelidae, and performing phylogenetic and dating analyses over the individual and the combined datasets. Our total-evidence phylogenetic analysis recovered the clade Chelidae as monophyletic and as sister group of a clade of South American extinct chelids; furthermore Chelidae retained inside the classical molecular structure with the addition of extinct taxa in both the Australasian and the South American clades. Our dating results suggest a Middle Jurassic origin for the total clade Pan-Chelidae, an Early Cretaceous origin for Chelidae, a Late Cretaceous basal diversification of both geographic clades with the emergence of long-necked lineages, and an Eocene diversification at genera level, with the emergence of some species before the final breakup of Southern Gondwana and the remaining species after this event.


Introduction
Pan-Chelidae is one of the two main lineages of crown Pleurodira (e.g., turtles that retract their neck inside the shell in a horizontal plane). Pan-Chelidae contains nowadays around 25 extinct species (Maniel & de la Fuente, 2016;de la Fuente et al., 2017a, b;Oriozabala et al., 2019) and approximately 58 extant species (Turtle Taxonomy Working Group, 2017) and shows a disjointed distribution in South America and Australasia. The oldest fossils that could be attributed to Pan-Chelidae are from the early Cretaceous (110 million years ago [mya]) from Argentina and Australia Lapparent de Broin & Molnar, 2001;Smith, 2010;de la Fuente et al., 2011), suggesting that they originated in the south region of the Gondwanan supercontinent (Broin & de la Fuente, 1993).
Until now, the phylogenetic relationships among extant and extinct chelids are still under debate. Since the XIX Century (e.g., Boulenger, 1889) morphological studies divided chelids into two groups regarding the length of the neck. One group is formed by the long-necked chelids, where the length of the neck is longer than the length of thoracic vertebrae. The other group is formed by the short-necked chelids where the length of the neck is shorter than the length of the thoracic vertebrae. Following Gaffney (1977) and subsequent morphological studies (Bona & de la Fuente, 2005;de la Fuente et al., 2015de la Fuente et al., , 2017aManiel et al., 2018), long-necked chelids form a monophyletic group (except in the study of de la Fuente et al., 2017a) represented by the extant South American species Hydromedusa tectifera Cope, 1870, Hydromedusa maximiliani Mikan, 1825 and Chelus fimbriata Schneider, 1783, the extant Australasian species belonging to the genus Chelodina, and the South American extinct taxa H. casamayorensis de la Fuente & Bona, 2002, Yaminuechelys gasparinii de la Fuente et al., 2001, and Yaminuechelys maior Staesche, 1929, and probably (because they have not being included in a phylogenetic analysis up to now), Chelodina alanrixi Lapparent de Broin & Molnar, 2001 and Ch. murrayi Yates, 2013. In this sense, the morphological hypothesis suggests that the origin of the long neck occurred only once in the evolutionary history of chelids. In addition, the biogeographic scenario under this hypothesis accounts for the diversification of each clade before the final breakup of Southern Gondwana (de la Fuente et al., 2015). On the other hand, the molecular phylogenies (Seddon et al., 1997;Shaffer et al., 1997;Georges et al., 1998;Guillon et al., 2012;Rodrigues & Diniz-Filho, 2016;Pereira et al., 2017) suggest that the Australasian chelids form a monophyletic group, including short and long-necked species, and the South American chelids form another monophyletic group also including short and long-necked species. Following this hypothesis, the presence of the long neck would have evolved independently in the South American and Australasian clades.
In spite of the significant amount of incompleteness surrounding the fossil record (Wilkinson, 1995;Wiens, 2003aWiens, , b, 2005Escapa & Pol, 2011;Wiens & Morrill, 2011), extinct taxa have become a key component of evolutionary studies because they contribute helpful information to elucidate the phylogenetic relationships of stem and crown groups and have intrinsic highly valuable information about the time of evolutionary events. However, this information is not always easy to understand. For example, the phylogenetic uncertainty, usually due to the loss of morphological information (i.e., non preserved parts), and the fossil age uncertainty may increase the lack of fit between phylogeny and stratigraphy. Consequently, this might lead to the existence of long ghost lineages as part of the best hypothesis for the evolutionary history of a group (e.g., Norell & Novacek, 1992;Benton & Storrs, 1994;Pol et al., 2004;Pol & Norell, 2006;Sterli et al., 2013). Regardless of these drawbacks, researchers involved in the study of the evolution of different groups of organisms, which contribute to the final goal of dating the "tree of life" (e.g., Philippe & Forterre, 1999;Benton & Ayala, 2003), have increasingly assimilated the temporal information from fossils in their research. Consequently, and in parallel with a wide range of advances in the field of molecular clock methods (e.g., Sanderson, 1997;Huelsenbeck & Ronquist, 2001;Sanderson, 2003;Drummond et al., 2006;Drummond & Rambaut, 2007;Yang, 2007;Drummond et al., 2012), extinct taxa have been widely used as node age calibrations (i.e., the "node-dating" method), and discussions about the shortcomings associated with the time estimates based on these calibrations became frequent (Ho & Phillips, 2009;Inoue et al., 2010;Ksepka et al., 2015;Parham et al., 2012). The main flaws of this method seem to be related to arbitrary decisions, including phylogenetic position of the extinct taxa used as calibrations (Near & Sanderson, 2004;Near et al., 2005;Marshall, 2008;Lee et al., 2009;Pyron, 2010), shapes and limits of prior age distributions (Marshall, 2008;Dornburg et al., 2011;Parham et al., 2012;Joyce et al., 2013;Ksepka et al., 2015;Warnock et al., 2015), number of calibrations and the inclusion of more than the oldest known fossil belonging to a lineage as calibration (Near et al., 2005;Marshall, 2008;Pyron, 2010;Dornburg et al., 2011;Duchêne et al., 2014;Gavryushkina et al., 2014).
The Mk model allowed Lewis (2001) to estimate phylogenetic trees from morphological data under Maximum Likelihood criterion for the first time, and later, Nylander et al. (2004) implemented this model to infer phylogenies from morphological data in Bayesian context. During the following years, the estimation of divergence times using Bayesian inference has become widely adopted, leading to the development of many algorithms, models and software (e.g., BEAST, Drummond & Rambaut, 2007;MCMCTree, Yang, 2007;DPPDiv, Heath et al., 2012;MrBayes, Ronquist et al., 2012a). Nevertheless, the direct integration of the information from extinct taxa by estimating their phylogenetic position and divergence dates at the same time remained virtually undeveloped until Pyron (2011) combined the Mk model and Bayesian relaxed molecular clock model to estimate divergence times in Lissamphibia from a "total evidence" matrix. This approach, later improved by Ronquist et al. (2012b) and referred to as "tip-dating", was applied by other researchers in subsequent studies (e.g., Near et al., 2014;Arcila et al., 2015;Cannatella, 2015;Bapst et al., 2016;Puttick et al., 2016).
One important aspect of the molecular (and morphological) clock development was the constant advance in the models describing branching processes (i.e., tree models). Some of them were developed to describe microevolutionary processes (at intraspecific or population level) like the Coalescent model (Kingman, 1982) and others were developed to describe macroevolutionary processes. These latter models include the pure Birth model (Yule, 1925;Gernhard, 2008), the Birth-Death model (BD, Gernhard, 2008)  This work aims to increase knowledge about the evolution of Pan-Chelidae. To carry out this goal we built new morphological and molecular datasets to analyze them individually and in combination, and produce new phylogenetic results on this group, as well as estimate the age of their main nodes using multiple methodological tools. Furthermore, we evaluated our results by comparing the estimated ages in each dating analysis with those presented in previous comparable studies.

Dataset
Complete and partial sequences of five mitochondrial genes (12s rRNA, 16s rRNA, COI, CytB and ND4) and three nuclear genes (R35, Rag1 and cmos) from 51 species of pleurodiran turtles (44 chelid species as ingroup plus 7 pelomedusoid species as outgroup) were compiled from GenBank (supplementary material S1) and integrated in individual matrices. Additionally, the RAG-1 matrix was completed with a sequence from Hydromedusa maximiliani sequenced in the Lab of Genetic Identification (IdeGen), dependent of the Instituto de Diversidad y Evolución Austral (IDEAus-CONICET). Later, each matrix was aligned with ClustalX (Larkin et al., 2007) using default parameters. These matrices were concatenated with SequenceMatrix 1.7.8 (Vaidya et al., 2011) resulting in a unique alignment of 7,815 base pairs (supplementary material S2). Next, we applied PartitionFinder (Lanfear et al., 2012) under linked branch lengths and greedy search algorithm, using the Bayesian Information Criterion (BIC) as the selection method for the statistically best-fit partitioning scheme. In addition, the codon position for coding genes was considered in the search (see supplementary material S3).
The morphological matrix (supplementary material S4) was taken from the study of de la Fuente et al. (2017b) and modified after careful revision (see supplementary material S3).
In order to perform the combined phylogenetic analysis, the molecular and the morphological matrices were merged into a total-evidence matrix using TNT v1.1 (Goloboff et al., 2008a, b), while the dating analysis based on the total evidence was performed after reading the separated datasets (morphology and molecules) with BEAST v2.4.3 (Bouckaert et al., 2014) and saved as a totalevidence xml file.

Phylogenetic analyses
The phylogenetic analyses were performed under MP criterion, using the software TNT v1.1 (Goloboff et al., 2008a, b), with 100 replicates of heuristic search under Random Addition Sequences (RAS) and Tree Bisection-Reconnection (TBR), keeping 100 trees per replicate. The extended Implied Weighting method (Goloboff, 1993(Goloboff, , 2014 was used during the analyses, computing the implied weights separately for each column (character), according to its homoplasy. This searching method was used because it has been recently demonstrated (Goloboff et al., 2018) that IW outperforms equally weighted parsimony and probabilistic model-based methods when it is applied to empirical (i.e., nonsimulated evolutionary rates) datasets.
Following Mirande (2019), five reference concavity values (k = 1, 5, 10, 15 and 20) were used in order to produce weighting strengths from 50 to 99% of the weight of a perfect (nonhomoplastic) character. After the searches, a stability criterion (based on the one described by Mirande, 2019) was used to select the most stable analysis. In this procedure, for each Most Parsimonious Tree (MPT), obtained under a specific k value, the topological distance for each of the other MPTs (obtained under the remaining k values) is calculated. Then, the average topological distance is calculated for each k value, and the analysis with the lowest average topological distance is selected. The topological distances were calculated with the R function "dist.topo" from the APE package (Paradis & Strimmer, 2004). This procedure was applied for the morphological, the molecular and the combined matrix, and when more than one MPT was obtained after the searching phase a strict consensus tree was calculated in order to obtain one tree per data matrix. Branch supports were calculated over each final tree by 1000 replicates of Bootstrapping, recording absolute node frequencies.

Dating analyses
In order to obtain a range of results based on different data, topologies and branching models (table 1), we performed three dating analyses: i) a total-evidence tip-dating (TE TD) analysis; ii) a morphology-based tipdating (M TD) analysis; and iii) a node-dating analysis (ND). The three analyses were run with partially constrained topologies, using the MP trees as a backbone (for references on constrained nodes see figs. 2 and 3) to avoid the inference of wrong clades as a result of incorrectly rooted phylogenies (see Springer et al., 2018). The TE TD and the M TD analyses were ran under the FBD model implemented in the "Sampled Ancestor" package (Gavryushkina et al., 2014(Gavryushkina et al., , 2016. The node-dating analysis was run under the BD model, with constraints for the Australasian and the South American clades, based on the oldest extinct taxa included in the TE TD analysis for each clade (Chelodina alanrixi and Yaminuechelys gasparinii, respectively). All dating analyses were carried out using BEAST v2.4.3 (Bouckaert et al., 2014), and the setting priors were defined in the BEAUti platform (Drummond et al., 2012). For all partitions, the clock model selected was the relaxed uncorrelated lognormal (Drummond et al., 2006).
The parameters included in the FBD model are: birth (i.e., speciation) rate (λ), death (i.e., extinction) rate (μ), probability of sampling extant species (ρ), and the fossil sampling rate (ψ). This model assumes that the fossilization is a constant-rate Poisson process. The priors for these parameters were estimated using available information from extinct and extant pleurodiran species using available methodological tools: fossil sampling rate was estimated using TRiPS (Starrfelt & Liow, 2016), birth and death rates were estimated using PyRate (Silvestro et al., 2014), and the probability of sampling extant species was estimated as the ratio: number of species present in the sample/number of known extant pleurodiran species. The calculated values for these parameters were used in the molecular clock analyses as starting points. As example, the TRiPS and PyRate output were used to define limits of uniform prior distribution on the relevant parameters, as well as the starting values (see supplementary material S3 and S5). The time of origin for the FBD process, in the TD analyses, and the root age in the ND analysis, were established as a uniform prior distribution, ranging from the lower limit of the age of Proterochersis robusta Fraas, 1913 (228.4 mya, Norian), to the upper limit of the age of Araripemys barretoi (100.5 mya, Albian). This period was selected because it covers from the oldest known stem turtle to the oldest known crown pleurodiran turtle.
The morphological character evolution was modeled under the Mk model (Lewis, 2001) with rate variation modeled using a discrete gamma distribution. The age calibration densities of clades defined by topological constraints during the ND analysis, were estimated by the CladeAge method (Matschiner et al., 2016), setting the first occurrence age according with the stratigraphic limits Downloaded from Brill.com12/17/2019 03:42:25AM via free access of the oldest fossil record assigned to each clade, and a sampling gap of 2 mys (as recommended in the program documentation). In the tip-dating analyses, the uncertainty associated with fossil ages was incorporated by editing the xml file, following the example in the BEAST2 user web site (http://beast2.org/ divergence-dating-with-sampled-ancestorsfbd-model/). The MCMC chains were run for 200 million generations, sampling every 1,000 or 5,000 generations. This number of generations produced ESS values of the parameters above 200. The convergence of estimated parameters on the stationary phase of the MCMC chains was checked using Tracer v1.5.0 (Rambaut & Drummond, 2009). The final values of divergence time estimations and its associated 95% HPD were summarized over the target phylogenies with TreeAnnotator v1.8.0, discarding the initial 10% of generations as burn-in, and selecting as node heights the median of the ages.

Phylogenetic results
Morphology.   1A). DNA. The phylogenetic analyses performed over the DNA matrix resulted in 928 MPTs (200 under k = 1 and k = 5; 295 under k = 10; and 112 under k = 15 and k = 20) ranging from 8,180 to 8,114 steps in length. The lowest average topological distance (19.08) was obtained with k = 20, so the strict consensus of the 112 MPTs was calculated and presented as the best hypothesis for the DNA dataset ( fig. 1B). This topology is similar to those recovered in previous molecular studies of the clade Chelidae. Each geographic clade was recovered as monophyletic, and both have groups of longnecked and short-necked chelids. In the case of the Australasian clade, the Chelodina group is represented by Chelodina oblonga as the sister taxon of a polytomy with the remaining species of Chelodina. In turn, this clade is the sister group of the Australasian short-necked species, a clade fully resolved with Pseudemydura umbrina at the base, followed by (Rheodytes leukops (Elusor macrurus, Falviemys purvisi)) + (Elseya, (Emydura, Myuchelys)).
Total-evidence. The phylogenetic analyses performed over the TE matrix resulted in 10,900 MPTs (3,200 under k = 1; 4,900 under k = 5; 100 under k = 10; 1,200 under k = 15; and 1,500 under k = 20) ranging from 8,418 to 8,348 steps in length. The lowest average topological distance (40.15) was obtained with k = 15, so the strict consensus of the 1,200 MPTs was calculated and presented as the best hypothesis for the combined dataset ( fig. 1C). The molecular signal mainly dominates this topology, however, the South American extinct chelids form a clade located in the stem of Chelidae as in the morphological tree. Furthermore, inside the clade Chelidae, both geographic clades are recovered as in the molecular tree, but the presence of extinct taxa produces some topological differences. The Australasian clade forms a polytomy including all Chelodina species and a branch leading to the clade of shortnecked species. The short-necked clade has almost identical structure to the same clade in the molecular tree, except for the inclusion of Birlimarr gaffneyi between Pseudemydura umbrina and the clade containing Rheodytes leukops, Elusor macrurus, Falviemys purvisi, and the genera Elseya, Emydura and Myuchelys. The South American clade also is almost identical to the one recovered in the molecular tree, differing from this only for the presence of a polytomy containing Yaminuechelys gasparinii, Y. maior and Hydromedusa casamayorensis as the sister group of extant species of Hydromedusa. As expected because of the combination of data sources, the branch support values were low for the stem-group and the node Chelidae (mainly, bootstrap supports from 1 to 25), and the Australasian clade, in general, was better supported than the South American clade ( fig. 1C).

Origin and diversification of Pan-Chelidae
The results of our node age estimations are presented and compared with estimations from four previous molecular clock studies, selected from the literature according to comparable conditions such as dataset, methods used to date the phylogenies and proportion of shared nodes. These studies are Dornburg Chelidae + South American extinct chelids (node 2). This node was recovered and dated only in the TE TD and the M TD analyses. The estimations produced by these analyses for this node were 146.82 mya (95% HPD = 121.67-177.12 mya) and 130.28 mya (95% HPD = 108.4-158.5 mya), respectively (figs. 2-4,  . 4, table 2). However, is interesting noting that this comparison could change if the position of Pseudemydura at the base of the clade is considered (i.e., the topology recovered by Kehlmaier et al., 2019).
Chelodina (node 5). This node represents the crown group of Australasian long-necked chelids. It was recovered and dated in all our analysis. The TE TD analysis (figs. 2-4,   Among previous studies, only Pereira et al. (2017) included the two extant species of Hydromedusa, but they did not recover the clade as monophyletic, so comparisons were not possible.
Despite differences in the dataset, type of molecular clock analysis, calibration schemes, branching models, statistical algorithms, software, taxonomic sampling, and others, our results were notably different from the previous analyzed studies only in the estimated age for Pleurodira, which our analyses dated as having a significantly more extensive ghost range (considering the age rages from the 95% HPDs), reaching ages much older than those proposed in these (and other) previous studies. However, taking into account the discrepancies in the analytical conditions mentioned above, is worth noting that: i) our analyses produced, in general, broader range ages, possibly as a consequence of the complexity of the branching models used here (BD and FBD), which defines the distribution of node age, instead of directly by the user; ii) as noted in previous tip-dating vs. node-dating studies (e.g., Arcila et al., 2015), nodes including extinct taxa (TD analyses), are generally older than nodes being calibrated (ND analyses) with the same extinct taxa (e.g., node stem Chelodina including C. alanrixi and C. murrayi, or stem Hydromedusa including H. casamayorensis; figs. 3 and 4); iii) some observed differences in age estimations between the M TD analysis and any other analysis, may be explained by the topological differences more than any other reason. The latter is the case for the Australasian and South American short-necked chelids. Although these nodes were recovered in the M TD analysis, the phylogenetic structure of the morphological tree is notably different to the rest of the analyses that recovered these nodes (figs. 3 and 4, nodes 6 and 10).

Remarks on the phylogeny of Pan-Chelidae
Previous phylogenetic analyses for Pan-Chelidae (or Chelidae) were based on exclusively morphological or molecular data, and in the case of analyses under the MP criterion, equal weights were always used. Here, we compiled and modified molecular and morphological matrices in order to produce more complete datasets, and to use the technique of extended implied weights under a wide range of analytical conditions to analyze these data separately, and (for the first time) in a context of total evidence.
Morphology. The selected tree obtained from the morphological matrix (k = 10; fig. 1A) is comparable with those presented in previous studies (Gaffney, 1977;Bona & de la Fuente, 2005;de la Fuente et al., 2015de la Fuente et al., , 2017aManiel et al., 2018) in several aspects. This tree shows two main clades: Chelidae and the clade formed by South American extinct chelids. Chelidae includes extinct and extant species. A similar structure is shown in the tree obtained by de la Fuente et al. (2017b). However, our tree has a broader taxonomic sampling and better resolution of the phylogenetic relationships within the stem group. In most of the previous morphological phylogenies, there is a clade of long-necked chelids, with the typical structure: Chelus (Chelodina (Hydromedusa, Yaminuechelys)) (or similar) in derived position, nested within a hierarchical structure defined by short-necked chelids located from basal to more derived positions. The morphological phylogeny of the present work ( fig. 1A) does not have this structure, the long-necked chelids do not form a monophyletic group, Chelus fimbriata is located at the base of the clade Chelidae, and the clade with remaining long-necked species is the sister group of the clade of short-necked species. In de la Fuente et al. (2017a) the clade of long-necked chelids is not recovered as monophyletic either, Mendozachelys wichmanni (a short-necked chelid) is related to Chelus fimbriata, and these are the sister group of the remaining long-necked species. In our morphological phylogeny, Mendozachelys wichmanni is recovered as a part of the stem group, as in de la Fuente et al. (2017b). In its original description (de la Fuente et al., 2017a), Mendozachelys wichmanni was described as having cranial and post-cranial features similar to those present in South American and Australasian short-necked chelids, and even in panpelomedusoid species. This unusual combination of anatomical features can be translated as homoplastic characters that make the phylogenetic position of this taxon hard to be inferred. It has been exposed that this kind of conflict can be resolved using implied weights (Goloboff, 1993) since homoplasy is considered in the calculation of character weights.
DNA. Our phylogenetic hypothesis for the molecular data ( fig. 1B) is very similar to the previous molecular phylogenies performed over the clade Chelidae or higher taxonomic groups that include chelid representatives (Seddon et al., 1997;Shaffer et al., 1997;Georges et al., 1998;Guillon et al., 2012;Rodrigues & Diniz-Filho, 2016;Pereira et al., 2017). Nevertheless is worth noting that there is a notably difference with the results obtained in the study of Kehlmaier et al. (2019). In that study a complete mitogenome was obtained for Australasian chelids from alcohol preserved specimens and the phylogenetic Total-evidence. Our study presents for the first time a phylogenetic analysis over the total taxon Pan-Chelidae based on a matrix combining morphological and molecular data. The result of this analysis is a consensus tree ( fig. 1C) with a topological structure nicely fitted to the stratigraphic and biogeographic information. The main nodes recovered in this phylogeny (node 2, Chelidae, the Australasian and the South American clades), as well as the distribution of extinct taxa and their relationships with extant taxa, display a good complement between the morphological and the molecular data. The molecular signal dominates the structure of the tree; however, the morphological data are crucial in recovering the position of extinct taxa. Between these taxa, seven were recovered as stem group (Prochelidella cerrobarcinae, Rionegrochelys caldieroi, Bonapartemys bajobarrealis, Lomalatachelys neuquina, Palaeophrynops patagonicus, Mendozachelys wichmanni, and Phrynops paranaensis), and six were recovered within the crown group: Chelodina alanrixi, C. murrayi and Birrlimar gaffneyi, as part of the Australasian clade; and Yaminuechelys gasparinii, Y. maior and Hydromedusa casamayorensis, as part of the South American clade. The location of some of these taxa, in the TE tree, is shared by previous phylogenies. As an example, among the taxa included in the stem group, four of them (Rionegrochelys caldieroi, Bonapartemys bajobarrealis, Lomalatachelys neuquina, Mendozachelys wichmanni) were also recovered as stem group taxa by de la Fuente et al. (2017b); among the extinct taxa included in the crown group, we recovered Birrlimar gaffneyi as sister taxon of the Emydura group s.l. as in Megirian & Murray (1999), and Yaminuechelys as related to Hydromedusa as in Bona & de la Fuente (2005) and several subsequent analyses (de la Fuente et al., 2015(de la Fuente et al., , 2017aManiel et al., 2018).
Considering novel results or discordances between our analysis and previous studies, we can state that: i) Prochelidella cerrobarcinae was recovered as sister taxon of Yaminuechelys in Sterli et al. (2013), as sister taxon of Yaminuechelys maior and Chelidae in Joyce et al. (2016), and as sister group of Bonapartemys bajobarrealis + Yaminuechelys in Ferreira et al. (2018), while we recovered this taxon at the base of the stem group; ii) in our TE tree C. alanrixi and C. murrayi were not located within a monophyletic clade with the remaining species of this genus, instead, the entire group Chelodina forms a polytomy; iii) H. casamayorensis was recovered as sister taxon of Yaminuechelys instead of being more closely related to the extant species of Hydromedusa (Maniel et al., 2018); and iv) Phrynops paranaensis was recovered within the stem group, as sister taxon of the clade ((Palaeophrynops patagonicus, Bonapartemys bajobarrealis), (Mendozachelys wichmanni, Lomalatachelys neuquina)), instead of being related to extant species of Phrynops (Wieland, 1923;Maniel & de la Fuente, 2016). The position of P. paranaensis is the same as the position recovered from the morphological analysis; however, this is not the case for H. casamayorensis, which is recovered within the genus Hydromedusa in the morphological analysis ( fig. 1A). The presence of molecular data for the extant species of Hydromedusa could produce an effect of long-branch attraction (Felsenstein, 1978), resulting in the ejection of H. casamayorensis from the clade Hydromedusa; however this kind of behavior can also be explained by interactions between characters from different sources and/or the effect of missing entries (TE matrix = 45.1% of missing data) (Maddison, 1993).

Origin and diversification of Pan-Chelidae
Given the fact that our TE TD analysis collects all the information from different data sources, extant and extinct taxa are present in the matrix and it was ran under the FBD tree model, we selected it as the best representative to discuss our most relevant results.
Pleurodira. Although the aim of this study is not focused on dating the origin of Pleurodira and the taxonomic sampling between Pan-Chelidae and Pan-Pelomedusoides is unbalanced, this node represents the root of our trees, and consequently we produced age estimations for it. The median estimated age for this node is 172.63 mya (Middle Jurassic). Even when this age is older than the proposed in the previous studies here analyzed, is worth noting that two of them (Joyce et al., 2013 andPereira et al., 2017) also dated this node as Middle Jurassic, and the remaining estimations (with exception for the 95% HPD lower limit from the analysis of Rodrigues & Diniz-Filho [2016]) are included in the 95% HPD range covered by our three analyses ( fig. 4,  table 2).
The estimations for the origin of crown Pleurodira is also the estimation for the split between Pan-Pelomedusoides and Pan-Chelidae. Joyce et al. (2016) gave a specific date for the origin of Pan-Chelidae based on an MP calibrated morphological phylogeny. This study proposed that Pan-Chelidae would have risen around the Early Cretaceous in the Southern part of Gondwana. This age is much younger than previously proposed (e.g., Joyce et al., 2013) and those obtained here.
Chelidae + South American extinct chelids. The TE TD median date calculated in our study was 146.82 mya (Late Jurassic). During the Late Jurassic, the Southern continents (South America, Africa, Antarctica, India, and Australia) were forming a unique continental mass called Gondwana. Following our results, this clade originated well before the breakup of Gondwana. However, the oldest member of this clade registered in the fossil record is Prochelidella cerrobarcinae from the Early Cretaceous of Patagonia (de la Fuente et al., 2011), which implies a ghost lineage of almost 40 My.
Chelidae. This node represents the origin of the crown group Chelidae, but when the molecular (or the TE) phylogenetic topology is considered, it also represents the split between the Australasian and the South American clades. This node is dated in the TE TD analysis in 112.08 mya (Early Cretaceous). By the end of the Early Cretaceous (around 100 Mya), Africa and South America completed their breakup (Nürnberg & Müller, 1991;Granot & Dyment, 2015), but South America, Antarctica, and Australia were still united. Following our TE phylogeny and the fossil record, the oldest member of this clade, is not registered in the fossil record until the Late Cretaceous with Yaminuechelys gasparinii from Patagonia .
Several studies have produced dates for this node as well, and between those, the studies selected for comparison in this study yielded similar ages to those estimated by our three analyses. The range of our analyses covers the estimations of all the previous analysis ( fig. 4) except for the 95% HPD lower limits of the studies of Dornburg et al. (2011), andRodrigues &Diniz-Filho (2016). Consequently, our results (considering the 95% HPD), in agreement with most previous studies, suggest that the basal diversification of the crown group Chelidae would have begun in the Early Cretaceous, and would have taken place until the Late Cretaceous, so this was an exclusively Cretaceous process.
Australasian clade. The median date for this clade was dated to 85.26 mya (Late Cretaceous) for our TE TD analysis. Australasian chelids diversified after the beginning of the opening of the Tasman Sea (90 mya) that was completed in the early Eocene (50 mya following Woodburne & Case, 1996) or in the late Eocene (35 mya following Lawver et al., 2011). Although this clade would have originated in the Late Cretaceous (as suggested by our study), it is not registered in the fossil record until the Eocene with the fossil Chelodina alanrixi (Lapparent de Broin & Molnar, 2001).
Chelodina. Our estimations for this node (mainly in agree with previous analyses) suggest that the extant species of Chelodina (crown group Chelodina) begun to diversify in the Oligocene (21.93 mya). However, considering the basal position of extinct species related to crown Chelodina (e.g., Chelodina alanrixi, Ch. murrayi) retrieved in our TD analyses (figs. 2-3), the Australasian longnecked chelids would have originated in the Early Eocene (around 53.61 mya). These estimations imply that this group would have radiated at the time of the Early Eocene Climatic Optimum (EECO).
Australasian short-necked chelids. Our TE TD analysis calculated that this clade originated near the K/Pg boundary. However, this clade is not registered in the fossil record until the Oligocene-Miocene with Emydura s.l. (Warren, 1969;Gaffney, 1981).
South American clade. The structure of this clade is given by the South American longnecked taxa Hydromedusa and Yaminuechelys as the sister group of Chelus plus South American short-necked taxa (Phrynops, Mesoclemmys, Rhinemys, Platemys, and Acanthochelys). The age of this node dated by our TE TD analysis is 101.86 mya (Early Cretaceous). At this point the Southern Gondwana was still united. This clade is not registered in the fossil record until the Late Cretaceous with Yaminuechelys gasparinii from Patagonia . The study of Pereira et al. (2017) gave very similar results for this node (mean = 99.3 mya; 95% HPD = 81.51-117.62).
South American short-necked chelids. The origin of this clade was calculated in our TE TD analysis in 50.98 mya (Early Eocene). The basal diversification of this clade occurred during the early Eocene, close to the EECO. The clade formed by Acanthochelys + Platemys + Phrynops + Mesoclemmys hogei have, in general, older diversification dates (Late Eocene and Oligocene) than the clade formed by the remaining spp. of Mesoclemmys + Rhinemys, that mainly diversified during the Neogene. The short-necked SA clade is not registered in the fossil record until the Late Oligocene/ Early Miocene with Phrynops indet. from Brazil (Kischlat, 1993).

The breakup of Southern Gondwana, the climatic changes, and the evolution of panchelids
The present and past distributions of panchelids are restricted to South America and Australasia. A similar disjointed distribution has also been observed in another clade of turtles, the Meiolaniformes (Sterli et al., 2015), and in numerous clades of extant and extinct taxa like the angiosperm Nothofagus, monotremes, and marsupials (e.g., Pascual et al., 1992;Broin & de la Fuente, 1993;Woodburne & Case, 1996;Sanmartín & Ronquist, 2004;Sigé et al., 2009;de la Fuente et al., 2011;Beck, 2012;Black et al., 2012). This particular distribution of organisms allowed Morrone (2002Morrone ( , 2009 to propose the Austral biogeographical kingdom, which is formed by southern South America, southeastern Australia, Tasmania, New Guinea, New Zealand, Antarctica, and southern South Africa. Furthermore, Sanmartín & Ronquist (2004) concluded that the distribution of the organisms in the Austral kingdom is linked to the breakup of Gondwana and that Antarctica was used as a path between Australasia and South America. The first landmass to separate from the rest of Gondwana was India. This separation started around 130 mya (Powell et al., 1988;Brown et al., 2006). Later, in the Early Cretaceous, Africa started separating from South America, losing the connection around 100 mya (Nürnberg & Müller, 1991;Granot & Dyment, 2015) entirely. In the Late Cretaceous (90 mya) the fragmentation of Southern Gondwana started with the origin of the Tasman Sea (between South America-Antarctica and Australia-New Guinea) and it was completed in the Early Eocene (50 mya following Woodburne & Case, 1996) or the Late Eocene (35 mya following Lawver et al., 2011). At the Eocene-Oligocene boundary (33.9 mya) the connection between South America and Antarctica was also lost, with the opening of the Drake Passage (Livermore et al., 2005). The isolation of Antarctica from South America and Australia, promoted the generation of the circumpolar current, causing a notable decrease in the global temperatures and precipitations (Zachos et al., 2001;Merico et al., 2008). Apparent changes in flora and fauna followed the climatic changes during the Cenozoic (Pascual, 1984;Ortiz-Jaureguizar & Cladera, 2006;Pascual & Ortíz-Jaureguizar, 2007;Woodburne et al., 2014). As previously suggested (e.g., de la Fuente et al., 2011) and based on our TE TD analysis, the origin and early evolution of pan-chelids predate the final breakup of Southern Gondwana near the Eocene-Oligocene boundary.
The Cretaceous was a crucial time for the origin and diversification of stem and crown chelids. Notably, the first main diversification events occurred during the Early Cretaceous in the stem chelids and the basal clades of crown chelids (split between South American and Australasian chelids and the split between Yaminuechelys + Hydromedusa clade, and Chelus + SA short-necked chelids). Joyce et al. (2016) suggested that the origin of South American and Australasian chelids is consistent with vicariant processes that took place during the Early Cretaceous. Later, during the Late Cretaceous a second diversification event occurred, including forms of stem chelids as well as the main basal clades of crown chelids (split between short and long-necked AU chelids, split between Yaminuechelys and Hydromedusa clades and split between Chelus and short-necked SA chelids). During the Late Cretaceous stem and crown chelids were present in Patagonia. All these diversification events could be correlated to the well-known Cretaceous Terrestrial Revolution (KTR) that is also recognized in other groups like mammals, crocodyliforms, squamates, dinosaurs, flowering plants, among others (Lloyd et al., 2008). Vlachos et al. (2018) recognized a "chelid turnover" in South America starting in the latest part of the Cretaceous (Campanian) that accelerated after the K-Pg extinction event. During the "chelid turnover" the basal lineages of chelids were replaced by derived lineages (Vlachos et al., 2018). In our TE phylogeny, the majority of stem chelids did not survive beyond the K-Pg boundary (except for "P." paranaensis), and only crown chelids survive until now, which could agree with this proposed "chelid turnover". However, "P." paranaensis does not follow the hypothesis in our analysis. Consequently, we consider that the position of "P." paranaensis among stem chelids should be taken by caution. Including more morphological characters in the matrix and a more detailed description of this extinct species could help to solve its phylogenetic position.
Our results suggest that few lineages originated during the Paleocene and Eocene, not only in South America but also, and maybe more evident, in Australasia. Furthermore, there are only three known species of chelids found in the Paleocene and Eocene of South America and Australasia (e.g., Chelodina alanrixi, Yaminuechelys maior, Hydromedusa casamayorensis). Vlachos et al. (2018) suggested that during the Paleocene-Eocene, there was a diversity crisis affecting South American turtles based on both the raw and phylogenetic diversities. During this diversity crisis, the number of species decreased. Our study seems to agree with the proposed "diversity crisis" recognized by Vlachos et al. (2018), and it would include not only the South American chelids but also the Australasian chelids.
After the Eocene and until the late Oligocene, there is a gap in the fossil record of around 20 mys in both continents. The next older fossil record of chelids for South America is from the late Oligocene/Early Miocene of Brazil (Phrynops spp., Kischlat, 1993) and for Australasia is from the Oligocene-Miocene (Emydura s.l., Warren, 1969;Gaffney, 1981).
Based on our total evidence phylogeny ( fig. 2), the Eocene-Oligocene time frame represents a critical time for the origin and early diversification of the main lineages of extant chelids.
The initiation of the circum-Antarctic current, promoted by the final isolation of Antarctica in the Eocene-Oligocene boundary, provoked a global decrease in temperature and precipitation (Zachos et al., 2001;Livermore et al., 2005;Merico et al., 2008). However, our TE TD analysis suggests that the origin and diversification of most extant genera occurred after this climatic deterioration. As turtles are ectothermic animals, the northern expansion of chelids in South America and the northern drift of the Australian landmass towards the north (towards hotter temperatures) could explain this diversification observed after the climatic deterioration in ectothermic animals, as we explain below.
Although nowadays chelids are not distributed in southern South America (Iverson, 1992), the fossil record of the clade shows Patagonia as an area where chelids have been present and diverse (de la Fuente et al., 2013;Vlachos et al., 2018). The pre-Eocene fossil record of chelids in South America is restricted (until now) to high paleolatitudes (Patagonia and the Mendoza province in Argentina; see Vlachos et al., 2018 for a summary), and the southernmost record of chelids in Patagonia seems to be the remains attributed to Testudines indet. by Otero et al. (2012) but here identified as Chelidae indet. (JS, pers. obs.). These remains come from the middle Eocene of Sierra Baguales in the Magallanes Region in Chile, more than 1600 km south of the extant distribution of chelids (Hydromedusa tectifera approximately 37°-38° South Latitude; Iverson, 1992). Interesting noting is that the remains from Sierra Baguales and Hydromedusa casamayorensis coming from the middle Eocene of Cañadón Hondo (Chubut province, Argentina) are the youngest records of chelids in Patagonia. The post middle Eocene record in South America is restricted to warmer and more humid areas outside Patagonia (de la Fuente et al., 2013). The fossil record is in accordance with the climatic deterioration registered after the Eocene-Oligocene boundary where the climate became drier and colder particularly in Patagonia (Zachos et al., 2001;Ortiz-Jaureguizar & Cladera, 2006). Previous authors have noticed this northern expansion in the fossil record of chelids (Broin & de la Fuente, 1993;Vlachos et al., 2018).
In Australia, the global cooling could have been buffered by the continental drifting towards the north (towards lower latitudes and higher temperatures) (McGowran et al., 2004). Nowadays, the southernmost distribution of chelids in Australasia corresponds to Chelodina longicollis Shaw, 1794 that reaches the southern tip of Victoria, Australia (39° South Latitude) (Iverson, 1992;Kennett, 2009). In the Australasian fossil record, the southernmost chelids were found in Oligocene-Miocene sediments in Tasmania, Australia (Warren, 1969), more than 400 km south to the extant distribution of chelids.
Early Cretaceous pleurodiran turtles from Australia: stem-chelids? or crown chelids? Smith (2010) described the oldest pan-chelid turtles from Australia coming from the middle Albian (Early Cretaceous) of Griman Creek Formation. Unfortunately, the fragmentary nature of the specimens and the absence of synapomorphic features do not allow them to be included in a cladistic analysis to test their phylogenetic position. Knowing the phylogenetic position of these pan-chelid turtles from Australia would have a significant impact on the evolution and paleobiogeography of pan-chelids ( fig. 2). We are going to discuss the leading two phylogenetic positions and their impact in the calibration of the pan-chelid tree.

1.
Early Cretaceous Australian pan-chelids as stem-chelids: This scenario proposes two hypothesis for the Early Cretaceous Australian pan-chelids: a, as part of the monophyletic clade today formed only by South American species, and b, the possibility that the Early Cretaceous Australian pan-chelids and the monophyletic stem-chelid turtles are successive paraphyletic groups with regards to crown chelids. Both hypotheses would have an impact on the calibration, affecting mainly the basal nodes such as Pan-Chelidae and Chelidae, pushing back the time of the origin of both clades. 2. Early Cretaceous Australian pan-chelids as crown chelids: This scenario proposes the hypotheses that the Early Cretaceous Australian pan-chelids are: c, located in the stem of the Australasian clade (more likely following the paleobiogeography patterns obtained in the total evidence analysis) or d, located in the stem of the South American clade. These hypotheses would have an impact on the time of origin mainly in the Australasian clade, South American clade, and Chelidae, pushing back their origin to at least the early Early Cretaceous.
Finding new, more complete pan-chelid specimens from Australia will undoubtedly provide valuable information to comprehend the origin and early evolution of pan-chelid turtles better.

Evolution of long necks in chelids
The fact that both current geographical clades contain long-necked species, and the absence of long-necked extinct species in the stem group, keeps alive the uncertainty about the origin of this trait in the group, even with the temporal framework provided by our study.
In light of current fossil evidence and according to our total-evidence time estimations, there are three possible and equally likely (or parsimonious in terms of cladistic optimization) hypotheses for the origin of the long neck: i) the long neck could have emerged at the base of the clade Chelidae during the Early Cretaceous, and later, reversed in the branch leading to the short-necked chelids within each geographical clade over the end of the Late Cretaceous; ii) the long neck could have emerged independently at the base of the South American clade in the Early Cretaceous and the base of the clade Chelodina in the Late Cretaceous or the Paleocene. Additionally there should have been a reversion to short-necks in the South American clade during the Late Cretaceous; iii) the long neck could have emerged at the base of the clade Chelodina in the Late Cretaceous or the Paleocene, and independently in each branch of the South American long-necked genera (clade Hydromedusa + Yaminuechelys and Chelus) during the Late Cretaceous.
Although the three hypotheses for the origin of the long neck are reasonable explanations and equally parsimonious, external evidence available up to now suggests that a single, common origin of the long neck in the group would seem unlikely. For example, until now there are no long-necked chelids in the stem group, the oldest known longnecked chelids (i.e., Yaminuechelys) are from the Late Cretaceous, and our age estimation for the origin of the clade containing these extinct long-necked taxa is also Late Cretaceous in age. However, we need to be cautious with these conclusions because new evidence could alter our hypothesis.

Conclusions
The total taxon Pan-Chelidae is a compelling case study, in terms of divergence time estimations, because of its current disjointed distribution, and the implications of splits before or after the final separation of South America, Australia and Antarctica. In the present study, we approach this topic by building a dataset as complete as possible, compiling and modifying molecular and morphological matrices to analyze these data sources both individually and in combination, using several analytical techniques of phylogenetic search, evolutionary rate estimations, and molecular/ morphological clock analyses.
Our molecular phylogenetic analysis (fig. 1B) mainly agrees with the topology recovered in previous molecular analyses; however our morphological analysis ( fig. 1A) shows some differences with previous morphological analyses. We present for the first time a total-evidence analysis for Pan-Chelidae. Primarily, it is worth noting that our TE phylogeny recovered as monophyletic the clades: Chelidae and Australasian and South American chelids. Some extinct taxa were recovered as a monophyletic group in the stem group, and some extinct taxa are located inside the crown group. The position of those extinct taxa agrees, in general, with the state of knowledge, from descriptive studies, and previous phylogenies.
Regarding the dating analyses, there were nodes with similar dates obtained by the three analyses, and also nodes dated with significant differences (figs. 3 and 4). However, considering the dates produced by analyzing all the information available up to now (TE TD analysis, fig. 2), we can hypothesize that: i) the origin of the clade formed by Chelidae + South American extinct chelids could have occurred in the Late Jurassic (146.82 mya); ii) the origin and diversification of Chelidae occurred during the Cretaceous; iii) the South American clade would have begun to diversify in the Early Cretaceous, before the Australasian clade, which would have begun to diversify in the Late Cretaceous; iv) basal radiations in the South American and the Australasian clades, with emergence of longnecked lineages, would have occurred before the K/Pg boundary; v) main basal diversifications in the South American and the Australasian clades would have occurred before the final breakup of Southern Gondwana; vi) diversifications in the South American and the Australasian clades at genera and species level would have occurred once the Drake Passage opened and the land masses were separated. These hypotheses would be tested in future works when new extinct and extant pan-chelids are included.