Materials

Specimens of Archaeopteryx in this study are designated through a commonly used numerical sequence that roughly corresponds to their succession of discovery2. The fifth specimen of Archaeopteryx2 (JM 2257) is a nearly complete and largely articulated skeleton of the smallest Archaeopteryx specimen known to date. It is also known as the Eichstätt Specimen and housed at the Jura Museum in Eichstätt, Germany (JM). The seventh specimen of Archaeopteryx2 (BSP 1999 I 50) is represented by a comparably complete skeleton exhibiting a substantial degree of articulation. It is formally named Solnhofen-Aktien-Verein Specimen but generally referred to as Munich Specimen, and is kept at the Paläontologisches Museum München in Munich, Germany (PMM). Skeletal elements of both the fifth and seventh specimen of Archaeopteryx have experienced brittle deformation during post-depositional compaction that resulted in splintering of the bone cortex. The ninth specimen of Archaeopteryx2 (BMMS-BK1a) preserves a partially disarticulated right wing skeleton of comparably large size that is presently housed at the Bürgermeister–Müller–Museum in Solnhofen, Germany (BMM). It is officially named “Exemplar der Familien Ottman & Steil”, also known as Bürgermeister–Müller Specimen, and colloquially referred to as “Chicken Wing”. Although a certain degree of post-depositional compaction is evidenced by the presence of several fractures that propagate through the long bone cortex, cortical splintering has not occurred. Its elements have therefore largely preserved their original three-dimensional geometry2.

Comparative material (Supplementary Data 1) was sourced from the collections of the European Synchrotron Radiation Facility, Grenoble, France (ESRF), the Musée des Confluences, Lyon, France (MdC), the Muséum national d’Histoire naturelle, Paris, France (MNHM), the Museum of Evolution, Uppsala, Sweden (MoE), and the University of Manchester, Manchester, England (TUoM).

Data acquisition

The humeral and ulnar cross-sectional geometry of the three specimens of Archaeopteryx, 28 species of neornithine birds, the small coelurosaur Compsognathus longipes, the rhamphorhynchid pterosaur Rhamphorhynchus sp., the anhanguerid pterosaur Brasileodactylus araripensis, and the crocodile Crocodylus niloticus were visualised through propagation phase-contrast synchrotron X-ray microtomography (PPC-SRµCT) at beamlines BM05 and ID19 of the European Synchrotron Radiation Facility. An ulnar cross section of aff. Deinonychus antirrhopus was imaged and subsequently paired with a humeral section of D. antirrhopus from literature32 through morphological and dimensional comparison.

Synchrotron X-ray tomography was conducted by utilising an optimised polychromatic beam with sufficient coherence to permit the application of PPC-SRµCT. Propagation phase-contrast imaging relies on a certain propagation distance between the sample and the detector that allows for the exploitation of the phase-contrast effect towards emphasising low-contrast features33. The fifth and seventh specimen of Archaeopteryx were imaged in accumulation mode, a novel acquisition protocol developed for imaging fossils encased in lithic slabs. The motivation for and implementation of the accumulation mode are explained in Supplementary Note 1 Further details of the adopted data acquisition parameters for each sample are provided in Supplementary Data 3.

Data processing

Three-dimensional volume reconstruction was conducted through filtered back projection following a phase retrieval protocol that relies on a homogeneity assumption by using a modified33 version of the algorithm developed by Paganin et al.34. Virtual two-dimensional cross-sectional slides were extracted directly from the reconstructed volumes at the developmental mid-diaphyseal plane oriented perpendicular to the local bone long axis in VGStudio MAX 2.2 (Volume Graphics, Heidelberg, Germany). Avian ulnae that carry a quill knob (ulnar papilla) at mid-diaphysis were virtually sampled at the level nearest to mid-diaphysis where no quill knob was present. The data set was supplemented with avian samples used in earlier studies8,10,35. Furthermore, scaled figures depicting complete perpendicular humeral and ulnar cross sections of avian and non-avian archosaurs were sourced from literature17,36,37,38,39 and processed in tandem with data obtained through PPC-SRµCT.

Cross-sectional geometry

Based on the characteristically unfolded nature of the Solnhofen Plattenkalk40, the geometry of Archaeopteryx wing elements was assumed to have experienced only brittle deformation during unidirectional compaction with insignificant movement of bone fragments perpendicular to the visualised cross sections. Two-dimensional restoration was conducted with image editing software by virtual extraction of the bone fragments and visually applying optimal fit of local fracture geometry, periosteal and endosteal curvature across adjacent fragments, and internal structures (e.g., canalisation). For the ninth specimen of Archaeopteryx, humeral and ulnar parameters were obtained by averaging the values found for two reconstructed circa mid-diaphyseal cross sections each. As the humeral and ulnar geometry of the fifth and seventh specimen are distorted to a markedly larger degree than those of the ninth specimen, they are represented by the single best-preserved cross section present in the circa mid-diaphyseal domain.

The elements of Compsognathus and Rhamphorhynchus used in this study were recovered from the Solnhofen Plattenkalk as well, and were reconstructed following the same protocol as the Archaeopteryx material (Supplementary Fig. 10). One fragment of cortical bone is conspicuously absent at the optimal sample location for the Compsognathus ulna in the upper right quadrant of the bone in the extracted slide, as also evidenced by an ulnar cross section extracted 3.58 mm proximal to the used sample location (Supplementary Fig. 10). The geometry of this cortical fragment at the sample location was reconstructed through close comparison with the bone and fracture geometry visible in the referred more proximal cross section.

All transverse cross sections were converted to binary cortical bone profiles by tracing the periosteal and endosteal surfaces and subsequently filling the area of the original cortical bone white41. Occasionally occurring spongy bone and obvious irregularities, such as cracks or protruding splints, were digitally removed. The area of the few small splints in fossil material that could not be accurately repositioned in their exact original orientation was taken into account during restoration of typically the periosteal margin. Cross-sectional geometric parameters were calculated with MomentMacro 1.4 (http://www.hopkinsmedicine.org/fae/mmacro.html) in the public-domain image analysis software ImageJ (https://imagej.nih.gov/ij/).

Most species in our data set are represented by the humerus and ulna of a single adult individual, although some were included as a composite of elements sourced from two individuals, or as average values derived from elements from two, three or four individuals (see Supplementary Data 1). Individuals sampled in this study are believed to represent adults based on element size and bone structure. Archaeopteryx specimens are often considered to be juveniles, which has been specifically concluded for the specimens of Archaeopteryx included in this study through relative size and bone surface texture2,13. The sampled Compsognathus specimen was also reported to represent a juvenile individual42. The studied Rhamphorhynchus individual is of comparably small size, suggesting juvenility as well. Gender composition across the data set is generally unknown and was therefore not considered.

Locomotor modes and body mass

Avian flight mode categorisation notoriously suffers from the qualitative, non-discrete nature of faunal flight strategies20. To overcome classification-specific effects in discriminant analysis, we considered the classifications suggested by Viscor et al.21 and Close et al.20 independently. Both avian flight mode divisions were expanded with one group that encompasses volant wing-propelled diving auks, and supplemented with alternative archosaurian locomotory strategies represented by exemplary taxa (Supplementary Data 1). The avian flight mode categories sensu Viscor et al.21 encompass (1) short flight, (2) forward flapping/bounding flight, (3) high-frequency flapping flight, (4) undulating flight and (5) gliding–soaring flight, which were assigned following the proposed taxonomical designations21. Taxa not included in their work were assigned flight modes according to the provided description21. Geococcyx californianus was classified as ‘short flight’ rather than as ‘forward flapping/bounding flight’ proposed for Cuculidae. A second, more recent avian flight mode division by Close et al.20 separates (1) burst flight, (2) intermittent bounding flight, (3) continuous flapping flight, (4) flap-gliding flight and (5) soaring flight, and was applied through description. We chose to score volant wing-propelled divers separately in both subdivisions as their aquatic locomotory strategy is known to profoundly influence wing bone morphology17,43 and, consequently, the expression of flight-related adaptations recorded therein43. Both referred avian flight classifications were complemented with the following locomotor categories: (6) long-tailed pterosaurian flight, (7) short-tailed pterosaurian flight, (8) (avian) non-volant wing-propelled diving, (9) ratite bipedal, (10) (non-avian) dinosaurian bipedal, (11) (non-avian) dinosaurian omnipedal and (12) crocodilian quadrupedal.

Body mass values for extant taxa were either directly available for the referred individuals or sourced from online databases44,45 and literature46,47 as species averages (see Supplementary Data 1). For extinct forms, either specimen-specific body mass estimates13,39,48 or average specific body mass estimates were available49,50. The Malagasy shelduck Alopochen sirabensis, reported to have been “slightly larger” than the Egyptian goose Alopochen aegyptiaca51 (average body mass of 1900 g44), was assigned a reconstructed body mass of 2000 g. Body mass for the Rhamphorhynchus sp. MdC 20269891 was reconstructed through the relation between body mass and wing span for basal pterosaurs proposed by Witton52. Total wing length for MdC 20269891 was measured as the cumulative length of the humerus (19 mm), radius (34 mm), wing metacarpal (14 mm), phalange I (47 mm), phalange II (40 mm), phalange III (35 mm) and phalange IV (44 mm) taken from photographic and scan data, and amounts to 233 mm. In the Dark Wing specimen of Rhamphorhynchus muensteri (JME SOS 4785; Jura Museum Eichstätt), the distance between left and right glenoid measures 1.56 × humeral length, which proposes an original interglenoid distance of 30 mm for for McD 20269891. Its corresponding wingspan, calculated as twice the wing length plus the interglenoid distance, amounts to 0.496 m. From the relation of Witton52 follows a reconstructed body mass of 95 g. The body mass for the Brasileodactylus araripensis individual in our study was inferred through close morphological and dimensional agreement between its humerus (length 168 mm, maximum distal width of 47) and the humerus of AMNH 22552 (length 170 mm, maximum distal width of 46 mm)53,54, for which a reconstructed wingspan of 3270 mm was reported55. From the described relation between wingspan and body mass in pterodactyloids56 follows a reconstructed body mass of 6540 g.

Body mass values for the studied specimens of Alligator mississippiensis (141 cm38) and the domestically bred Crocodylus niloticus (200 cm; personal observation PT) were reconstructed through specific allometric scaling relations between body length and body mass offered in literature57,58.

Cross-sectional parameters

Relative cortical thickness12 (CA/TA) and mass-normalised resistance against torsional forces12 (J/M) were quantified for archosaurian humeri and ulnae (Supplementary Data 1, Supplementary Figs. 4 and 5). CA/TA describes element hollowness as the ratio of cortical bone area to the total area delimited by the external surface of the bone in cross section (Supplementary Fig. 2). As such, CA/TA is proportionate to the corticodiaphysary index (CDI)59 and inversely related to the K-parameter6,10. Polar moment of inertia of an area J quantifies the mechanical resistance against torsion around the longitudinal axis of the considered element. J mathematically equals the sum of the maximum second moment of area (Imax) and minimum second moment of area (Imin) that quantify resistance against deflection along the respective orthogonal major and minor principal axes (Supplementary Fig. 2) through the relative distribution of matter12. Values for J obtained from cross sections with an Imax/Imin >1.50 are typically overestimated12,60,61, but remain informative when considered proportionally rather than quantitatively (as is its derivative Zp41,61). J was normalised over body mass to permit comparison in a highly body mass-diversified comparative framework that spans well over five orders of magnitude (Supplementary Data 1).

Cortical vascular density, expressed as the amount of canals per mm2 of bone area in cross-section62,63,64,65, was considered qualitatively for a modest selection of archosaurs for which high-resolution data was available, but not challenged statistically (Supplementary Data 1). Bone area in section was calculated as CA (see ‘Cross-sectional geometry’ above) with MomentMacro 1.4 in ImageJ 1.49. In the fifth and ninth specimen of Archaeopteryx, cortical canals were counted visually. Absolute canal abundance in the cross-sectional cortex of archosaurs other than Archaeopteryx was obtained through selection of the darkest grey levels in the greyscale histogram (including canals) by thresholding the cortical domain and subsequently counting the amount of elements within canal size range using the Analyse Particles function of ImageJ 1.49.

The ratio of Imax over Imin provides a reliable measure for the ellipticity of the transverse bone shaft and has been considered as such in biomechanical explorations12,60,66,67. These approaches traditionally assume that the degree and orientation of ellipticity reflect an adaptation that offers optimised resistance against bending, with the direction of Imax corresponding to the orientation of the maximum bending moment. However, an opposite functional interpretation of cross-sectional element ellipticity in which a preferred bending direction is achieved through orientation of Imin has also been proposed specifically for avian wing bones68. Such conflicting explanations of the same parameter illustrate the complexity of interpreting cross-sectional bone ellipticity in a functional context and thereby obscure the information offered by other characters when assessed in a multivariate context. We therefore chose not to include quantified bone ellipticity measures in our comparative study.

Tree inference and divergence chronogram

Mesozoic topology and timing used in this study (Supplementary Data 3 and Supplementary Fig. 1) were derived from the Paleobiology Database (PaleoDB; https://paleobiodb.org) on 29 January 2016. Divergence nodes were adopted as the older bound date for the oldest report of a taxon nested beyond the respective split. Mesozoic terminal nodes and the Tertiary terminal node for Mancalla cedrosensis were placed at the younger bound date of their occurrence. Alopochen sirabensis is placed at 656 AD, which represents the median calibrated radiocarbon age for the last-occurrence date for the species69. The 19th century terminal node for Pinguinus impennis was dated through its well-documented last observation in 1844 AD70. Topology and timing within the extant avian subset were largely adopted from the well-resolved phylogeny by Jarvis et al.71 (Supplementary Data 4 and Supplementary Fig. 1), since the more recent neoavian phylogeny proposed by Prum et al.72 was found to conflict PaleoDB on numerous crucial accounts. Several higher-order divergence times in Aves were obtained from PaleoDB following the procedure described above. Two specific inconsistencies in PaleoDB were negotiated through literature (see Supplementary Note 3). Insufficiently resolvable divergence nodes were placed at a standard + 4 MY with respect to their closest established crownward node. The three Archaeopteryx specimens were included as a polytomy at + 4 MY with respect to the older bound date for the genus in recognition of taxonomic and ontogenetic uncertainty2. The phylogenetic tree was constructed in Mesquite 3.0473.

Statistical analyses

The relations between individual geometric parameters and locomotor divisions in the training taxa set were statistically assessed through phylogenetic analysis of covariance using the PDAP module of Garland et al.74. For each parameter, 10000 unbound simulations were performed along the constructed tree (Supplementary Fig. 1) under a Brownian motion regime in PDSIMUL. ANCOVA was performed with a grouping of the training taxa according to their locomotor classes as response variable, parameter values as predictor variable, and body mass as covariate (Supplementary Data 5).

Phylogenetic PCA75 scores for the studied taxa, founded on humeral and ulnar CA/TA and J/M (Supplementary Data 2), were obtained with the phyl.pca function (method: BM; mode: cor) of the phytools package76 in the R-environment77 through RStudio 0.99.48478. The phylogenetic PCA scores were subsequently subjected to Partitioning Around Medoids specified to two clusters with the pam function of the cluster package79 in RStudio 0.99.484.

The archosaurian groups outside Archaeopteryx serve as training taxa that represent known locomotor modes and thus form a morphological reference environment for discriminant analysis. The optimum value of Pagel’s λ, the scaling factor of autocorrelation for a certain parameter on a given phylogenetic tree80 to be applied in phylogenetically informed discriminant analysis, was found using the approach described by Schmitz et al.81 in RStudio 0.99.484 (Supplementary Fig. 8). Linear discriminant analysis and classification of mystery taxa (Supplementary Data 2) was conducted in PAST 3.1082, as was one-way MANOVA (Supplementary Tables 1–3) among individual locomotor strategies. Additional motivation for and information on the statistical approach used here is available as Supplementary Note 4.

Data availability

All data underlying the study are available in Supplementary Data 1.