Mosaic anatomy in an early fossil squamate – Nature

Preservation and taphonomy of NMS G.2023.7.1

Disarticulated bones of NMS G.2023.7.1 are spread over a diameter of approximately 19 cm on an undulating bedding surface. This preservation is similar to that of specimens of P. estesi (NHMUK (Natural History Museum, London, UK) PV OR 41388) and cf. P. estesi (NHMUK PV R8851), from the Early Cretaceous Purbeck Limestone Group⁸. Various elements are visible at the surface of the block, including a right mandible, braincase, various other skull bones, vertebrae, ribs, a partial right coracoid, humeri, a right ilium and right femur (Fig. 1b). Numerous additional bones are present in the matrix (Fig. 1b). Of the skull, NMS G.2023.7.1 includes the left jugal, left postfrontal, left and right parietals, left and right squamosal, braincase, left vomer, left palatine, right pterygoid, right dentary, right angular, and right compound bone incorporating the surangular, prearticular and articular. Of the axial skeleton NMS G.2023.7.1 includes 32 vertebrae or partial vertebrae, many dorsal ribs, and a cervical intercentrum (Extended Data Fig. 5 and Supplementary Data 2). Of the appendicular skeleton, NMS G.2023.7.1 includes the right coracoid, right and left humeri missing distal ends, right ilium, fragment of right pubis, right femur missing distal end, left femur missing epiphyses, right tibia missing epiphyses, one metapodial, and three phalanges, including one ungual phalanx. Details of all these elements, including links to 3D digital models, are included in Supplementary Data 1.

The morphologies and lack of duplication of all squamate elements preserved in NMS G.2023.7.1, are consistent with assignment to a single relatively large-bodied individual (see main text). Although the bones are disarticulated, they are spread in a rostral to caudal pattern with skull elements and anterior vertebrae at one side of the block, and hindlimb and caudal vertebrae at the opposite edge, with dorsal vertebrae, ribs, and forelimb elements between (Fig. 1b and Extended Data Fig. 5). NMS G.2023.7.1 is isolated from other skeletons recovered from the same bedding surface by at least two metres distance. Those skeletons belong to mammals, amphibians and fish, but not squamates. Thus, it is unlikely that any squamate elements reported here represent distinct squamate individuals or species, transported from adjacent areas of the lagoon floor.

The presence of large numbers of easily transported Voorhies Group 1 elements (for example, vertebrae and ribs⁵¹) indicates that the skeleton was subjected to currents that were neither strong enough to remove these elements, nor strong enough to transport allochthonous Group 2 and 3 elements (that is, limb bones and skull elements⁵¹) from other individuals into the bone scatter before it was buried. Other non-squamate vertebrate elements such as bone crumbs and a tritylodontid tooth are present in the region of the skeleton, but these Group 3 lag elements are most likely to have been present on the lagoon floor when the skeleton was deposited, and similar material is very common throughout vertebrate-bearing levels of the sequence and around other skeletons.

These observations, the rarity of squamate remains compared to other groups in the Kilmaluag Formation assemblage¹⁴, and the rarity of large squamate remains specifically (of almost 200 specimens collected from the Elgol Coast SSSI from 2014–2024, NMS G.2023.7.1 is the only specimen to preserve any large-bodied squamate remains) provide strong support for the view that NMS G.2023.7.1 represents a single squamate individual.

Reconstruction of skull and body proportions

The skull reconstruction shown in Fig. 3 was conducted in Blender 3.5.0 by arranging elements of NMS G.2023.7.1 in 3D space, with reference to extant squamate anatomy. The outline of the maxilla was based on the maxilla of NHMUK PV OR 48388, the holotype of P. estesi. The skull length is estimated at 41.4 mm and primarily uses information from the preserved portions of the dentary and compound bone. The reconstructed skull has long, low proportions. Evidence for this comes from the braincase dimensions, relative to the lengths of the combined palatal elements and mandibular elements. We allowed additional vertical height at the back of the skull to accommodate slight crushing of the braincase.

The life reconstruction shown in Fig. 1a was produced by Mick Ellison at the American Museum of Natural History, in consultation with R.B.J.B. and S.E.E., using measurements derived from the specimen. Measurements were made from 3D digital models, using Meshlab 2023.12⁵², and are reported in Supplementary Data 2. The length of the presacral vertebral column was estimated based on the summed lengths of the 23 definite presacral vertebrae with measurable lengths (84.5 mm, excluding the condyles), giving an estimated presacral length of 99.2–110.2 mm if 27–30 presacral vertebrae were originally present (allowing for the missing atlas and the possibility of missing cervicals or dorsals). Of this, the summed cervical lengths give an estimated neck length of 25.9 mm or more. The straight-line length of the longest complete dorsal rib (19 mm) and ilium length (15 mm) informed reconstruction of body depth, and the partial humerus length (10.1 mm), femur length (19 mm) and tibia length (11.5) informed reconstruction of limb lengths. Other aspects, such as hand and foot morphology, are not informed by evidence and should be considered as generalized.

Discovery, preparation and imaging of NMS G.2023.7.1

NMS G.2023.7.1 was discovered by S.A.W. in March 2015, at the Elgol Coast SSSI, during fieldwork led by R.B.J.B. and S.A.W., assisted by A. Wolniewicz, with permission of the landowner, the John Muir Trust, under permit from NatureScot (then, Scottish Natural Heritage). It was extracted as a block of micritic limestone approximately 220 mm long, 180 mm wide and 150 mm deep. This block was embedded in silicone and prepared from behind using acetic acid by S. Moore-Fay of Wavecut Platforms, to its current thickness of 15–30 mm.

We scanned the full slab, using the Nikon Metrology XT H 225 ST X-ray μCT scanner at the School of Earth Sciences X-ray Tomography Facility, University of Bristol, UK, providing a pilot scan of the whole specimen with all bones in their original positions. Segmentation of this and other scans in the current work was conducted using the software Mimics 19.0 (Materialise) primarily by E.F.G. The scan and parameters for this and all other μCT scans described in the current work are available on MorphoSource (links provided in Supplementary Data 1).

The pilot CT scan is the basis of the digital map of the specimen shown in Fig. 1b and Extended Data Fig. 5 and was also used to identify regions of matrix free of preserved bones that were then removed with a table-mounted disc cutter. The resulting reduction of slab diameter allowed a better signal-to-noise ratio in subsequent episodes of CT scanning. We also used the pilot scan to separate some portions of the specimen in smaller blocks for high-resolution μCT scanning at the University of Bristol facility. These blocks, and the remaining portion of the specimen, were given subpart numbers that follow from NMS G.2023.7.1: (1) NMS G.2023.7.1.1, the main part of the slab, excluding the following sections; (2) NMS G.2023.7.1.2, a small portion including the braincase, left jugal, right parietal, left postfrontal and a cervical vertebra; (3) NMS G.2023.7.1.3, small portion including the left palatine and a dorsal rib; (4) NMS G.2023.7.1.4, small portion including the left vomer; and (5) NMS G.2023.7.1.5, small portion including a tritylodontid tooth and unidentified bone fragments (grey elements Fig. 1b, bottom right).

We also targeted regions of the remaining slab (NMS G.2023.7.1.1) for phase-contrast synchrotron X-ray tomography on beamline ID19 of the European Synchrotron Radiation Facility (ESRF), Grenoble, France (described below): (1) the right dentary, angular and compound bone, right pterygoid, right squamosal, left parietal and left humerus; (2) the right ilium, partial pubis and femur, phalanges including an ungual phalanx, dorsal ribs, and eight vertebrae, including dorsals, caudals and a cervical intercentrum; and (3) the right humerus and scapulocoracoid, plus a phalanx and a cervical and two dorsal vertebrae.

Finally, the main slab (NMS G.2023.7.1.1) was split into four portions that were scanned separately at School of Earth Sciences, University of Bristol. This allowed higher quality 3D models of a few elements for which models from our other scans were not of sufficient quality, including the tibia, atlas and some other vertebrae.

We also completed CT scans of the holotype (NHMUK PV OR 48388) and referred (NHMUK PV R8851) specimens of P. estesi using a Nikon XTEK H 225 ST MicroCT scanner at Cambridge Biotomography Centre, University of Cambridge, UK, also available via MorphoSource (Supplementary Data 1).

Phase-contrast synchrotron X-ray tomography

For synchrotron X-ray tomography, the beamline was set up for filtered white beam (W150-B wiggler gap 39 mm; filtered with 10 mm Cu and 0.5 mm W) resulting in a total integrated detected energy of approximately 170 keV. Images were recorded with an indirect detector comprising a 500 µm LuAG scintillator, a set of two Hasselblad lenses (100 and 150 mm; Victor Hasselblad) set for a 0.67× magnification, and a PCO.edge 4.2 sCMOS camera (PCO), resulting in a measured pixel size of 8.96 µm. The sample-detector distance was set to 13.2 m for propagation phase contrast. We used 5,000 projections over a 360° rotation, with an exposure time of 0.06 s per projection (taking the average of 3 frames of 0.02 s each), 41 flat-field images and 40 dark-field images were used as calibration. The recorded field of view in this configuration was 8.89 mm vertically and 18.35 mm horizontally (992 × 2,048 pixels). The area scanned during each rotation was increased by shifting the centre of rotation by around 8 mm horizontally (corresponding to 900 pixels on the detector), allowing us to reconstruct tomograms across a field if view spanning 34.59 × 34.59 mm (3,861 × 3,861 pixels). We then combined scans of multiple fields of view to image wider areas. Data were processed using PyHST2 with the single distance phase retrieval approach^53,54. Post-processing included a ring artefact correction⁵⁵, change of dynamic range from 32-bits to 16-bits using the maximum and minimum 0.001% histogram clipping values of all the datasets of the series, and weighted averaging of duplicated tomograms resulting from the overlap on the vertical axis.

Osteohistology

Three elements were thin sectioned from the associated partial skeleton NMS G.2023.7.1. The midshafts of a partial humerus and partial femur (Fig. 2c,f), and unknown region of the body of a rib were manually prepared from the specimen block with a carbide needle. The resulting pieces of bone were embedded in EpoThin 2 low viscosity epoxy resin. Embedded specimens were cut transversely at or nearest the mid-diaphysis using an Isomet 1000 precision saw. After being mounted onto frosted glass slides with clear Gorilla superglue gel, specimens were ground to optical clarity (~30 µm) on a lapidary wheel and photographed in plane-polarized light on a Nikon Eclipse LV100POL microscope fitted with a Prior ProScan III automated stage adaptor. Composite images were processed using Nikon NIS-Elements BR (version 5.24.03) imaging software using the Extended Depth Focus function to autofocus and z-stack photomicrographs to improve image quality of microstructural details. Two thin sections were made of the humerus and femur, and one thin section was made from the rib. High-resolution photomicrographs are available from MorphoSource (Supplementary Data 1).

Phylogenetic analyses

We included NMS G.2023.7.1 and two other parviraptorid specimens in three phylogenetic datasets. Parviraptorids were represented in these datasets by three operational taxonomic units (OTUs): (1) the holotype of Breugnathair elgolensis (NMS G.2023.7.1); (2) the holotype of P. estesi (NHMUK PV OR 48388); and (3) the referred specimen of cf. P. estesi (NHMUK PV R8551). Parviraptorid specimens from other localities were not included in the analyses as they are less complete (or incompletely described) but are discussed in the Supplementary Discussion. Whiteside et al.⁴⁸ recently reported the Late Triassic Cryptovaranoides microlanius as an anguimorph, and therefore deeply nested within the squamate crown group. However, Brownstein et al.⁵⁰ presented highly differing anatomical observations and suggested that Cryptovaranoides may instead be an archosauromorph. Given the ambiguity of published anatomical data^48,49,50, we did not test the affinities of Cryptovaranoides here and await more well-resolved anatomical data. We also did not include the proposed stem lepidosaur Taytalura⁵⁶, also because of phylogenetic uncertainties³. Note that we closely considered the phylogenetic scores proposed by Caldwell et al.⁷ for parviraptorids when scoring these matrices (explained in Supplementary Discussion).

Phylogenetic inference was carried out using Bayesian inference in MrBayes 3.2.7a⁵⁷, using a fossilized birth–death tree prior^58,59 with a proportion of extant species sampled of 0.038 under diversified sampling, and default priors for speciation rate, extinction rate and fossil sampling rate. We used an offset exponential tree age prior with a minimum age of 240 Ma (Middle Triassic) and mean age of 250 Ma (Early Triassic), and a relaxed clock transition rate model with independent gamma rates with a default variance increase prior of 10 and a log-normal clock rate with a mean of −2.56 log units and variance of 1.08. The ages of all operational taxonomic units were specified using a uniform distribution between their minimum and maximum possible stratigraphic ages, extended from ref. ³ to encompass the wider set of fossil taxa included in the analyses conducted here. We constrained the ingroup relationships of squamates to reflect current consensus based on molecular phylogenetic analyses², by specifying a series of partial constraints matching those shown by ref. ³. The resulting backbone constraint is consistent with the phylogeny of ref. ², with trichotomies representing areas of uncertainty among recent analyses (that is, enforcing no specific relationships between the three constituent groups). For example, Gekkota and Dibamidae in an unresolved trichotomy with the clade including all other squamates. The toxicoferan clade (Iguania, Anguimorpha and Serpentes) was left as an unresolved trichotomy in analyses of datasets 1 and 2, but constrained to represent three different phylogenetic hypotheses of toxicoferan ingroup relationships for three separate analyses of dataset 3 (explained below).

Dataset 1

The early reptile dataset of Talanda et al.³ (modified from refs. ^32,45,60), which includes an extensive sample of early members of the reptile crown group, as well as living and fossil rhynchocephalians and squamates, focussing on deep lepidosaur and deep squamate divergences, but includes relatively few snakes: three extant species plus the early fossil snakes Najash, Pachyrachis and Dinilysia. To this we added our three parviraptorid OTUs, the hypothesized early anguimorph Dorsetisaurus, and the Late Jurassic squamate Eoscincus⁵, resulting in a dataset of 125 taxa and 382 characters.

Dataset 2

An extended version of the squamate dataset of Meyer et al.⁶ (modified from Gauthier et al.²⁷). Meyer et al.⁶ extended the taxon and character sample of Gauthier et al.²⁷ by adding characters and taxa relevant to Jurassic squamate divergences, including new fossil taxa. The resulting matrix includes a strong sample of early squamates and relevant characters, as well as a strong sample of snakes. To this we added fifteen OTUs: the possible stem lepidosaurs Velbergia bartholomaei, Fraxinisaura rozynekae, Paliguana whitei, Sophineta cracoviensis and Marmoretta oxoniensis, the early rhynchocephalians Gephyrosaurus bridensis and Kallimodon pulchellus, the stem squamates Bellairsia gracilis, Oculudentavis naja, Oculudentavis khaungraae, the possible early anguimorph Dorsetisaurus purbeckensis, and our three parviraptorid OTUs. This resulted in a dataset of 637 characters and 168 tips, compared to the 155 tips included by Meyer et al.⁶.

Dataset 3

An extended version of the dataset of Zaher & Smith⁴⁷ (modified from Hsiang et al.⁶¹; after ref. ²⁷), to which we added various early squamate and non-squamate fossils. The matrix of Zaher & Smith⁴⁷ included a large sample of living and fossil snakes and relevant characters, as well as 11 iguanians, 21 anguimorphs and two rhynchocephalians. This sample was intended to focus on toxicoferans, with relevance to snake ingroup relationships.

We modified the matrix of Zaher & Smith⁴⁷ by adding many species of non-toxicoferan squamates, adopting the scores of ref. ⁶² (modified from ref. ²⁷ based on revisions proposed in more recent works such as refs. ^63,64) when, as in most cases, the identical characters were used in these studies. We also added fossil taxa intended to increase representation of the stem groups of Lepidosauria, Squamata, and squamate subgroups: the possible stem lepidosaurs V. bartholomaei, Taytalura alcoberi, F. rozynekae, P. whitei, S. cracoviensis and M. oxoniensis, the stem squamates B. gracilis, Huehuecuetzpalli mixtecus, O. naja and O. khaungraae, the early-diverging Jurassic and Cretaceous squamates Hongshanxi xiei, Yabeinosaurus tenuis, Dalinghosaurus longidigitus, Liushusaurus acanthocaudata, Scandensia ciervensis, Meyasaurus faurai, Jucaraseps grandipes, Chometokadmon fitzingeri, Ardeosaurus brevipes, Eichstaettisaurus schroederi, D. purbeckensis, and our three parviraptorid OTUs. We omitted the candidate stem snake Tetrapodopis amplectus from analyses due to substantial doubts about its anatomy and phylogenetic affinities^65,66.

In addition to this, we added two characters. Character 789: gastralia present (0); absent (1). Character 790: entepicondylar foramen present (0); absent (1). This resulted in a dataset of 790 characters and 189 tips, compared to 788 characters and 90 tips by Zaher & Smith⁴⁷. Because analyses of this matrix returned parviraptorids as toxicoferans we ran three separate analyses of this dataset, using topological constraints that represent different phylogenetic hypotheses for toxicoferan ingroup relationships: (1) snakes constrained as sister to anguimorphs; (2) snakes constrained as sister to iguanians; and (3) anguimorphs and iguanians constrained as sister taxa, to the exclusion of snakes.

Evaluating convergence

Analysis of dataset 1 was run for 133 million generations, dataset 2 was run for 160 million generations, and dataset 3 was run for up to 180 million generations prior to convergence (at an average standard deviation of split frequencies of 0.01 or less). Chains were sampled every 10,000th generation, with a burn-in of 50%. The effective sample size was greater than 200 for all tested parameters, and an average potential scale reduction factor was 1.01 or less on all parameters, indicating convergence.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.