Phylogenomics and the position of turtles
Previous phylogenetic investigations of amniote phylogeny have failed to reach a clear consensus on the phylogenetic position of turtles, as the various studies have often produced ambiguous and sometimes conflicting results. The causes for this probably stem from the intrinsic difficulty of this phylogenetic problem, which involves ancient divergences. Most of the previous molecular studies addressing this question used either small datasets based on a few nuclear genes [19, 24] or genetically linked mitochondrial genes [22, 23, 37]. In general, phylogenetic analyses based on using mitochondrial data tended to recover a sister group relationship between turtles and Archosauria [20–22, 37, 38], whereas some of the nuclear data favoured a clade of turtles with crocodiles [23, 24, 26]. The only exception to this pattern is the study by Iwabe et al. , who reported statistical support for the turtles plus archosaurs clade, but this was based on only two nuclear genes and a minimal taxon sampling.
Resolving the branching patterns of the major lineages of amniotes requires gathering a considerable amount of informative sequence data from independent markers with adequate taxon sampling. Shen et al.  recently investigated this question using 23 (mostly nuclear) markers for 28 vertebrates, and estimated that with their taxon sampling, a total sequence length of more than 13,000 nucleotides from independent nuclear markers is necessary to provide significant statistical support for resolving the controversial relationships between turtles, birds, and crocodilians. Our phylogenomic results, based on 248 nuclear markers, corroborate their predictions about the challenge represented by resolving this phylogenetic question, and add support to the sister group relationship between turtles and archosaurs (birds plus crocodilians). Furthermore, our statistical analyses reject any alternative hypotheses to the sister group relationship of turtles to Archosauria (Table 2), thus advancing the resolution of this long-standing controversial issue of vertebrate evolutionary history.
However, as illustrated by the occurrence of conflicting signals in our phylogenetic analyses, the phylogenomic approach is not immune to statistical inconsistency , as highlighted here in the cases of ML and Bayesian analyses of nucleotide datasets under a single concatenated GTR + G model, and in species tree inference from nucleotide gene trees, which showed inconsistencies that are most likely due to saturation at third codon positions. The fact that turtles group with crocodilians in concatenation analyses is probably due to a long-branch attraction (LBA) artefact causing the faster evolving squamates to be attracted towards mammals and the outgroups (see Additional file 2). The same grouping of turtles + crocodiles was also retrieved with strong support by Tzika et al.  based on ML and Bayesian analyses of amino-acid datasets using site-homogeneous empirical models and a reduced taxon set. Those authors evoked the same kind of LBA artefact to explain what they considered as an artefactual result, as support for it disappeared in analyses including fewer sites but with fewer missing data . These observations confirm that phylogenomic reconstruction can lead to artefacts, especially when the taxon sampling is reduced and model assumptions are violated . When analysing large concatenations, use of best-fit models is required to account specifically for heterogeneities among genes and codons in the substitution process, and to alleviate the deleterious effects of substitution saturation. Similarly, we found that when using mixed models for analysing protein-coding gene concatenations, partitioning by codon position outperformed the widely used gene-partitioning scheme. The CAT-GTR + G mixture model nevertheless offers the most efficient solution to account explicitly for site heterogeneities in the substitution process as typically observed in phylogenomic datasets .
As illustrated by our results, statistical inconsistency is not restricted only to concatenation-based phylogenetic reconstruction methods. Although specifically designed to accommodate potential gene-tree discordances, species tree inference methods also seem to be sensitive to mis-specification of the substitution model used to infer gene trees. Indeed, the species tree obtained from the nucleotide dataset also strongly supported the artefactual grouping of crocodiles and turtles (Figure 3b). Therefore, in addition to their accounting for gene-tree heterogeneity, the use of the best-fitting substitution models seems to be equally important for these methods . These results also indicate a potential problem of stochastic error in reconstructing gene trees for which only a limited number of sites is available, and the consequent effect on species tree inference. Ultimately, species trees can only be as good as the gene trees from which they are inferred.
Finally, it is worth noting that a recent analysis of miRNA phylogenetic distribution  supported a branching order (turtles + squamates) that was strongly rejected by our data. This is not the first instance of a conflict between miRNA and sequence-based phylogenetic studies, as shown by the case of acoels for instance [43, 44]. Thus, our study suggests caution is needed when using miRNA markers in phylogenetic reconstructions, as they might not be as free from homoplasy as sometimes considered . For example, secondary loss of multiple families of miRNAs have already been reported in tunicates . Our results imply that the four miRNAs families exclusively shared by turtles and lizards  either have been lost secondarily in archosaurs, or have been independently recruited in turtle and lizard genomes. Upcoming reptile genomic data [47, 48] should allow verification of these predictions.
Consequences for interpreting character evolution in amniotes
The well-supported sister group relationship of the turtles to the archosaurs has important implications for the evolution of morphological, developmental, and ecological characters of amniotes. It implies, as previously proposed , that turtles experienced a secondary closure of the temporal fenestration, which therefore appears to be a derived character, rather than a reflection of the ancestral condition, as has long been assumed. In addition, because a basal position of turtles within reptiles is supported by the timing of events in organogenesis , our results suggest the occurrence of significant convergence in developmental timing characters, and advocate for the re-interpretation of sequence heterochrony data in the light of the phylogenetic position of turtles supported by our phylogenomic analyses. Finally, the assumption of an aquatic origin of turtles (the hypothesis that was brought forward due to the proposed sister group relationship between turtles and an extinct group of marine reptiles (Sauropterygia) and Lepidosauria ) also needs to be reconsidered. A recent study suggested, for example, that stem turtles could have occupied both terrestrial and aquatic habitats .
The proposed phylogeny of amniotes also provides a more solid background from which to investigate the evolution of the sex-determining systems and genomic characteristics of reptiles. Whereas mammals and birds have only genetic sex determination, non-avian reptiles have both genetic and temperature-dependent sex determination. Temperature-dependent sex determination also occurs in crocodilians, tuatara, and the majority of turtles, whereas it is less common in squamates [51, 52]. Studies on this subject have relied on a traditional view of the vertebrate phylogeny, with turtles being basal to the other reptiles (including birds) (compare, for example Janzen & Krenz  with Janes et al. ). Although the phylogenetic scenario supported by our data would not change the main conclusion that temperature-dependent sex determination evolved multiple times within amniotes, it does provide a basis from which to further investigate the possible adaptive evolutionary value of the temperature-dependent sex determination in amniotes and the evolution of sex chromosomes.
Recently, Matsuda et al.  reported a high degree of conservation between the chromosomes of a turtle (Pelodiscus sinensis) and the common chicken, in accordance with an archosaurian affinity of turtles. These authors also suggested that although no specific sex chromosomes could be identified in the studied turtle, which has genetically determined sex, the ancestral avian Z sex chromosome has been conserved in the turtle genome. However, other features of the genome, such as its average genome size, GC content, and distribution of transposable elements show a marked similarity between turtles and crocodiles, to the exclusion of birds [29, 55]. By rejecting the turtles plus crocodilians grouping, our analysis could possibly be interpreted as evidence for a parallel evolution of these genomic features in the two lineages, or, perhaps more plausibly, recent evolution of bird-specific features in the avian lineage.
Models of molecular clock relaxation
In our molecular dating analyses, we found discrepancies between the results obtained using standard substitution models (LG + G and GTR + G) and the CAT-GTR + G mixture model with both amino acids and nucleotides. These differences probably stem from underlying differences in branch-length estimates between the two types of models, indicating the need to apply the most appropriate models of sequence evolution currently available for conducting molecular dating analyses . Our results indicate that the CAT-GTR + G mixture model better accounts for the site-specific heterogeneities of our concatenated protein-coding gene datasets. Therefore, we consider that the divergence date estimates obtained under this model are the most reliable.
However, these small discrepancies between estimates obtained under different substitution models are almost negligible as compared with the large differences in estimates between the auto-correlated and uncorrelated models of rate change. In our case, the use of uncorrelated models generally led to unreasonably recent dating estimates for all nodes relative to the values reported in the literature . These results seriously question the adequacy of the uncorrelated models of molecular clock relaxation parameters for estimating divergence times, at least with our dataset. Based on Bayes factor comparisons, Lepage et al.  showed that auto-correlated models provide a significantly better fit than the uncorrelated gamma model, especially for large phylogenomic datasets. These results were recently confirmed in an empirical phylogenomic study focusing on the estimation of arthropod divergence times, for which the assumption of rate autocorrelation seemed to be the most realistic way of modelling evolutionary rate variations across the tree [56, 57]. For these reasons, we consider the results from the auto-correlated relaxed clock analyses under the CAT-GTR + G substitution model as our most reliable dating estimates (Table 3; Figure 4).
Our auto-correlated relaxed clock analyses based on the CAT-GTR + G model support a divergence between turtles and Archosauria around 255 Mya (274-233 Mya), which is in agreement with the estimates recently reported by Shen et al. . The dating obtained for other nodes also seems to be mostly compatible with current knowledge. For example, the Testudinoidea MRCA corresponding to the divergence between Emys and Chelonoidis is estimated at a mean of 91 Mya (range 131 to 54 Mya), which is comparable with that obtained by Lourenço et al. . Exceptions concern squamates and crocodilians, for which our estimates indicated a more recent time than generally reported (Table 3). We note that the confidence intervals are relatively large, however, as would be expected for such indirect estimates, in which dates are estimated jointly with the process of substitution-change variations over time .
The single major difference between our estimates and the previously published divergence dates concern the MRCA of living turtles. This is a controversial issue in the paleontological literature, with proposed ages of divergence between the two main turtle lineages (Pleurodira and Cryptodira) varying from the Upper Triassic to the Upper Jurassic. This debate is mostly due to the rarity and the need for better characterization of turtle fossils older than 160 Mya . Our analyses suggest that the MRCA of Chelonia (Pleurodira plus Cryptodira) is on average 157 Mya (range 207 to 104 Mya), taking the average mean and extremes of 95% credibility intervals for the CAT-GTR + G model amino-acid and nucleotide results. This means that an Upper Triassic origin (229 to 200 Mya) of extant turtle lineages is rejected, and that although a Lower Jurassic origin (200 to 176 Mya) is still possible, it seems unlikely. Remarkably, a similar conclusion was reached in a recent study using a fully independent molecular dataset, which only included turtle sequences and within-turtle fossil calibrations . Our 95% credibility intervals for the turtle ancestral node (185 to 104 Mya with amino acids, and 207 to 120 Mya with nucleotides) are nevertheless wider and probably less realistic in including the Lower Cretaceous (146 to 97 Mya). However, taken together, these two analyses begin to suggest that the origin of Chelonia may be in the Middle or Upper Jurassic (176 to 146 Mya) or later. If so, two controversial fossils, Proterochersis and Kayentachelys, attributed respectively to the Cryptodira and Pleurodira clades, would be currently misplaced on the turtle lineage. These placements, proposed by Gaffney , have been a subject of intense debate [62, 63] (references therein). This would also have important implications for molecular clock analyses, because these fossils are usually used for calibrating both turtle  and amniote  trees. Considering this, it may be prudent to consider these fossils Testudines incerte sedis until additional data can be obtained to confirm their phylogenetic placement.