The Arnstadt Formation of Saxony-Anhalt, Germany has yielded some of Germany’s most substantial finds of Late Triassic tetrapods, including the sauropodomorph Plateosaurus and the stem-turtle Proganochelys. Here, we describe an almost complete skull of a new sphenodontian taxon from this formation (Norian, 227-208 Ma), making it the oldest known articulated sphenodontian skull from Europe and one of the oldest in the world. The material is represented by the dermal skull roof and by the complete maxilla and temporal region, as well as parts of the palate, braincase, and lower jaw. A phylogenetic assessment recovers it as a basal sphenodontian closely related to Planocephalosaurus and to Eusphenodontia, making it the earliest-diverging sphenodontian known with an articulated skull. Its cranial anatomy is generally similar to the well-known Diphydontosaurus avonis from the Rhaetian of England, showing that this successful phenotype was already established in the clade around 10 myr earlier than long assumed. An analysis of evolutionary change rates recovers high rates of evolution in basal sphenodontians, with decreasing rates throughout the evolution of the group. However, contrary to previous studies, reversals in this trend were identified, indicating additional peaks of evolutionary change. These results improve our understanding of the early sphenodontian diversity in Europe, providing critical information on evolutionary trends throughout the history of the clade and sparking renewed interest in its evolution.

(a)   Computed Tomography

The holotype MB.R.4520.2 is housed in the fossil reptile collection of the Museum für Naturkunde in Berlin (Germany). µCT scans were performed in house using a Phoenix|x-ray Nanotom S tomography machine (GE Sensing and Inspection Technologies GmbH, Wunstorf, Germany) with a voltage of 110 kV and a current of 150 µA, capturing 2,000 images with an exposure time of 750 ms and a resolution of 10.23 µm. Slices were reconstructed using the datos|x v.2.3.0.844 - RTM reconstruction software (GE Sensing and Inspection Technologies GmbH, Phoenix—x-ray) and the resulting volume was manually segmented and analysed in VG Studio Max 3.3 (Volume Graphics GmbH, Heidelberg, Germany) at the 3D-Visualisation Laboratory at the Museum of Natural History Berlin using the region grower and pen tools.

(b)   Time-calibrated phylogenetic analysis

A sphenodontian phylogeny was constructed to explore the relationship of Parvosaurus harudensis within Rhynchocephalia using the morphological matrix by [23]. In the analysis performed by [24], incompletely resolved phylogenies are prepared for the change rates analyses by breaking all polytomies at random, inevitably losing phylogenetic information. We solved this problem by running preliminary analyses and identifying rogue taxa to be excluded a priori, thus resolving all polytomies based on highest likelihood. Sophineta cracoviensis and the Clevosaurus complex (except for C. hudsoni) were pruned from the parsimony phylogenetic and evolutionary rates analyses. This, with the addition of Parvosaurus harudensis, resulted in a matrix comprising of 30 taxa and 131 characters.

To provide a timescale for the phylogeny, first and last appearance dates (FAD and LAD, respectively) were compiled for each taxon from the Paleobiology Database (see Supplementary Data). The age of Parvosaurus harudensis is based on the Arnstadt Formation in which it occurs, and so the upper and lower boundaries of the Norian were used as LAD and FAD, respectively. Two different dating approaches were performed to facilitate assessment of the stability of the time-scaling and evolutionary rates methods. Following the methods and settings recommended by [72] and [73], branch lengths were computed in R v.4.1.0 (R Core Team, 2021) using the ‘equal’ and ‘minimum branch length’ [mbl] methods. For the ‘equal’ method, the root-length was fixed to a 2 myr duration. For the ‘mbl’ method, it was set to a minimum duration of 1 myr. This way of time-calibration was chosen over a Fossilised Birth Death model through Bayesian Inference to allow for more direct comparison to the evolutionary change rates analysis performed by [24]. In order to track the influence of Parvosaurus harudensis on the phylogeny and evolutionary change rates, all analyses were performed twice, once including, and once excluding Parvosaurus harudensis.

(1)   Bayesian analysis

The phylogenetic analysis was performed with the settings from [23] under the Mk model [74], applying a gamma distribution for variable rates. No outgroup was set. The analysis was run for 10⁷ generations with 2 runs for 4 chains each. Parameters were sampled for every 1,000 generations and relative burn-in was set to 50%. Posterior node probabilities for all Bayesian phylogenies were exported from FigTree v.1.4.4. All recovered phylogenies can be viewed in the Supplementary Data.

(2)   Maximum parsimony analysis

A phylogenetic maximum parsimony analysis was performed in TNT v.1.5 [75] under implied weights, with a concavity index (K) = 3 (all analyses performed with higher k -values resulted in the same singular tree). As for the Maximum Likelihood rates analyses, a single, fully resolved tree was beneficial, Parsimony was preferred as input for the evolutionary rates analyses over Bayesian Inference methods as seen in [23]. All branches were collapsed to zero length prior to the analysis. Trees were searched using 1,000 replications of Wagner trees with a single random seed and Tree Bisection Reconnection (TBR), saving 10 trees per each of the 1,000 replications. This resulted in the output of a single, fully resolved tree with a fit of 21.80595. Character state changes were counted for each branch using TNT v.1.5 and mapped on the phylogeny using the “Branch length” function (see Supplementary Fig. S10).

[c] Evolutionary rates analysis

Evolutionary rates were assessed using Maximum Likelihood methods following the protocols of [24, 73, 74, 76, 77]. The time-calibration and evolutionary rates analyses were performed in R v.4.1.0 (R Core Team, 2021) with the Claddis [77], paleotree [78] and tidyverse [79] packages with a modified version of the script by [24] and [73]. All analyses were repeated and averaged over 100 repetitions. To address uncertainty arising from sampling each terminal taxon’s age, node ages were randomised for each of the 100 replicates within the range of estimated ages. LRT significance testing was chosen over AIC for better comparability to the script by [24] with an alpha threshold of 0.01 used to evaluate significance, and with Benjamini-Hochberg false discovery rate correction. Partitions were used for time-bins, character partitions, clade partitions and branch partitions (see Supplementary Data). Time-bins were partitioned into geochronologic ages ranging from the Wuchiapingian (Upper Permian) to the Maastrichtian (Late Cretaceous). Characters were partitioned into cranial and postcranial characters. Evolutionary rates of the time-calibrated phylogeny were assessed over 100 runs with randomised dates and the resulting rates were illustrated by averaging over the 100 runs.