The external surface of non-avian dinosaur eggs is not usually smooth like those of their avian descendants. Unique ornamentation patterns sculpt the exterior of the eggs, a trait that is difficult to interpret because of its scarcity in modern taxa. One avian species that does homoplastically, present similar external eggshell ornamentation to that of non-avian dinosaurs is Dromaius novaehollandiae Latham, 1790, the emu. Here we use D. novaehollandiae eggs in conjunction with a clutch of oviraptorosaurian dinosaur eggs (NCSM 33576, Macroelongatoolithus carlylei) to test new methods of quantifying external eggshell ornamentation. Currently, the only scientific language for describing and comparing ornamentation styles in fossil ootaxa is restricted to qualitative categorization, which introduces issues of subjectivity and overly broad and overlapping typification. In this study, we derived and tested a new statistical quantitative approach to quantifying ornamentation that includes two existing functions of the molaR package in R previously applied to shape quantifications of fossil teeth, and ‘Orientation’, a novel function presented as a proxy for ‘direction’, needed to capture directionality. Results demonstrate that 1) the quantitative approach provides statistical backing to gross qualitative observations;  2) statistically significant differences exist between the ornamentation in D. novaehollandiae and M. carlylei, particularly in terms of relief; and 3) intranest variation of M. carlylei can be demonstrated from harmonic mean p-value differences between different pairs of eggs. This method offers a strong platform to consolidate quantitative measures with existing qualitative categories, improve the diagnoses of ootaxa, and answer broad ecological and evolutionary questions regarding dinosaur reproduction. Moreover, wider application of the technique is encouraged for a multi-proxy quantitative analysis of any paleontological or biological surfaces.

To compare eggs, we measured morphometrics (egg length and width) using digital calipers. We weighed D. novaehollandiae eggs using digital scales, and calculated egg volume using the equation from Hoyt [45]. We digitized the eight fossil and six extant eggs using the application Metascan on an Apple iPhone 12 Pro, following the methods of Avrahami and Herzog [46]. We coupled the iPhone camera with a 15x macro lens and a 10 cm diameter ring light. Meshes built this way achieved comparable resolutions to μCT scans (mean pixel size of c. 11 microns), and provided a methodology that could compare D. novaehollandiae eggs with M. carlylei that could not be CT scanned at high resolution in its current nest configuration. These high-fidelity scans also compare favorably to other portable surface scanners for the relatively low relief and high detail of eggshell surface texture versus other paleontological surfaces. Completed scans were exported as .gltf files from the application, preserving all raw data including texture and coloration. These files were imported to Blender (version 3.4). Here, we subsampled the eggshell surfaces so that we had uniformly-sized (10mm diameter) sections that could be analyzed in molaR, to avoid the effects of taphonomic cracks or non-representative breaks on the surface of the eggs that could alter their original surface topography. We extracted surface sections across the eggs by adding 10 mm diameter cylindrical meshes through the surface and using the Boolean modifier to intersect cylinders with the egg surface (Fig 2A). The sections were exported as .stl files using the batch export feature and were subsequently imported into MeshLab2022.02 where they were aligned on the X-Y axis (i.e., with a vertical Z axis; Fig 2B) and decimated to 5000 triangular faces using the Quadric Edge Collapse Decimation tool. To import these into R for molaR analysis, we exported aligned surfaces as .ply files. We binned the samples from each egg into one of five zones along the long axis of the egg, consistent with those described in [47] and matching eggs as they were laid in nidus. It is worth noting that this means Zone 1 in M. carlylei is the blunt end of the egg, and Zone 5 is the acute end, but in D. novaehollandiae Zone 1 is the acute end and Zone 5 the blunt end.

To quantitatively analyze the topography of each surface, we made use of the R package molaR (version 5.3) [48]. MolaR is an analysis tool with a suite of functions developed initially for the 3D analysis of dentition, and our study represents the first use of the analysis beyond this scope. Each .ply mesh was analyzed using two of the functions existing in molaR: Dirichlet Normal Energy (DNE) [49] and Slope [50]. These two metrics provide proxies for surface ‘complexity’ and ‘relief’ respectively. We also modified the Orientation Patch Count function [51] to write a new function called ‘Orientation’, which analyzes the relative surface area of faces occupying each of the eight 45 directional bins (Fig 2C) to calculate both the direction and strength of orientation (pole-to-pole along the egg long axis or around the egg short axis; full details below). We used pairwise, nonparametric, two-sample Kolmogrov-Smirnov tests to test the statistical significance of differences between distributions, and Welch’s two-sample T-test for mean differences, between the D. novaehollandiae eggs and M. carlylei and nests. For each metric, we used both full datasets and datasets averaged per egg zone to reduce the likelihood of sampling size bias. Mean values for DNE, Slope, and Orientation are reported for each zone to represent trends across the egg, as an average of sub-sampling within these zones that was intended only to give uniform, unaltered meshes. We calculated the inter-relationship of metrics using Pearson’s correlation coefficient. We calculated harmonic mean p-values [52,53] to meta-analyze the p-values found for the three metrics. We selected a harmonic mean over a meta-analysis like Fisher’s method [53] to account for the non-independence of DNE, Slope, and Orientation values.

‘Orientation’ is a new function in the molaR package [48] that draws on the existing Orientation Patch Count function (‘OPC’). The existing OPC function calculates the number of patches (where a patch is three or more triangles on a 3D mesh) that fall into each 45-degree bin about the X-Y plane (Figure 2D). This measure alone does not distinguish patch count size because a patch can incorporate three adjoining faces but also orders of magnitude more faces, and these are calculated by OPC as a single patch of equal value. To account for directionality, we use the surface area value from the ‘Patch Details’ output of this function, which calculates the total surface area occupied by faces in each bin. The function then calculates the ratio of surface area in the bins oriented along the long axis (bins 1, 4, 5, and 8) to those oriented around the short axis (bins 2, 3, 6, 7). An additional calculation is added to this ratio to produce a scale from -100 to +100 that indicates both the direction (positive along the long axis; negative around the short axis) and strength (higher numbers correlate to stronger directionality) of the orientation. This function will be available in the next update of the molaR package. We include functions both for individual meshes (‘Orientation)’ and a batch export to a .csv file with the DNE and Slope values, for the Statistical Complexity, Orientation, and Relief of Eggshell (SCORE).