Taphonomic deformation, the distortion of fossils as a result of geological processes, poses problems for the use of geometric morphometrics in addressing paleobiological questions. Signal from biological variation, such as ontogenetic trends and sexual dimorphism, may be lost if variation from deformation is too high. Here, we investigate the effects of taphonomic deformation on geometric morphometric analyses of the abundant, well known Permian therapsid Diictodon feliceps. Distorted Diictodon crania can be categorized into seven typical styles of deformation: lateral compression, dorsoventral compression, anteroposterior compression, ‘saddle-shape’ deformation (localized collapse at cranial mid-length), anterodorsal shear, anteroventral shear, and right/left shear. In simulated morphometric datasets incorporating known ‘biological’ signals and subjected to uniform shear, deformation was typically the main source of variance but accurate ‘biological’ information could be recovered in most cases. However, in empirical datasets, not only was deformation the dominant source of variance, but little structure associated with allometry and sexual dimorphism was apparent, suggesting that the more varied deformation styles suffered by actual fossils overprint biological variation. In a principal component analysis of all anomodont therapsids, deformed Diictodon specimens exhibit significant dispersion around the ‘true’ position of this taxon in morphospace based on undistorted specimens. The overall variance associated with deformation for Anomodontia as a whole is minor, and the major axes of variation in the study sample show a strong phylogenetic signal instead. Although extremely problematic for studying variation in fossil taxa at lower taxonomic levels, the cumulative effects of deformation in this study are shown to be random, and inclusion of deformed specimens in higher-level analyses of morphological disparity are warranted. Mean morphologies of distorted specimens are found to approximate the morphology of undistorted specimens, so we recommend use of species-level means in higher-level analyses when possible.

Empirical data collection and analyses

522 crania of Diictodon feliceps (composed of 518 specimens from South Africa, three specimens from Zambia, and one specimen from China) were examined for this study. Of these specimens, 485 were complete enough to have landmarks digitized on at least one side of the skull in dorsal view, with 387 of those complete enough for a bilaterally-symmetric configuration of dorsal landmarks, and 464 specimens were complete enough to have landmarks digitized in lateral view. Specimen images were digitized using ImageJ and TPSDig (Abràmoff, Magalhães & Ram, 2004; Rohlf, 2010).

In dorsal view, we digitized a configuration of 16 landmarks consisting of four midline landmarks and six pairs of bilaterally-symmetric lateral landmarks: (1) anterior edge of premaxilla; (2 & 11) Prefrontal-lacrimal sutural border at orbital margin; (3 & 12) Anteroventral edge of postorbital bar; (4 & 13) Posteroventral edge of postorbital bar; (5 & 14) Posterior extent of squamosal; (6 & 15) Anteromedial edge of temporal fenestra; (7 & 16) Prefrontal-frontal sutural border at orbital margin; (8) Mid-frontal sutural border with preparietal; (9) Anterior edge of pineal foramen; (10) Mid-parietal sutural border with postparietal.

We analyzed two permutations of the dorsal view dataset. In the first dataset (‘dorsal’), we reflected bilaterally symmetric landmarks across the midline and averaged the positions of the resulting pairs of landmarks to create a series of ‘half specimens.’ In cases where a symmetric landmark was not preserved on one side of the skull, the coordinates of the single preserved landmark were used. This approach follows common practice for dealing with incomplete specimens in paleontological datasets, but the reflecting and averaging process can affect both the biological and taphonomic signals in a dataset including distorted specimens (e.g., elimination of natural asymmetry, or creation of misleading mean values when one side of the skull is highly sheared relative to the other) (Angielczyk & Sheets, 2007). Therefore, we also utilized a second dataset (‘dorsal complete’) consisting of only those specimens in which all sixteen landmarks could be digitized. The resulting bilaterally symmetric landmark configurations were then utilized in subsequent statistical analyses without reflecting and averaging symmetric landmarks.

For the lateral view, we digitized 11 landmarks: (1) Anteroventral tip of premaxilla; (2) Anterior edge of canine/caniniform process; (3) Septomaxillary-nasal sutural border at narial margin; (4) Posterior edge of canine/caniniform process; (5) Prefrontal-lacrimal sutural border at orbital margin; (6) Prefrontal-frontal sutural border at orbital margin; (7) Ventral edge of maxillary-jugal suture; (8) Postorbital-postfrontal sutural border at orbital margin; (9) Anteromedial edge of temporal fenestra; (10) Posterior extent of squamosal; (11) Posterior edge of parietal. Images of specimens photographed in left lateral view were reflected prior to digitization. If a specimen could be digitized for all lateral landmarks in left and right views, the mean landmark coordinates of the two images were used for subsequent analysis. Otherwise, only the single complete side was utilized. 

            Specimens were grouped by four variables: sex, size class, assemblage zone (AZ), and deformation style (see Supplementary Data). Sex was determined by the presence or absence of maxillary tusks, following Sullivan et al. (2003) in considering tusked individuals to be male and tuskless individuals to be female (although accurate identification of whether the tusked cohort represents males or females is not important for the purposes of our analyses, only the existence of a dimorphic pattern in the sample). In the smallest observed tusked specimen of Diictodon (BP/1/102, total skull length 4.87 cm), the tusks are just erupting—all known smaller skulls are tuskless. Because the sex of these presumed juveniles cannot be determined, they were excluded from analyses using this variable. The specimens excluded from analyses of sex make up the ‘small’ size class, consisting of all specimens with a total skull length less than 5 cm. The ‘medium’ size class consisted of specimens with total skull length ranging between 5 and 9 cm, and the ‘large’ size class consisted of specimens with total skull length in excess of 9 cm. Assemblage zone data were available for 419 specimens used in the dorsal analysis, 337 specimens used in the dorsal complete analysis and 391 used in the lateral analysis. Assemblage zone data were derived from specimen labels, Haughton & Brink (1954), Kitching (1977), Smith (1993), and Rubidge (1995). Specimens were considered of ‘unknown’ assemblage zone and excluded from zonal analyses if they lacked locality data altogether, had only vague locality information (e.g., “Cape Province”), or were collected without stratigraphic context at a locality known to span multiple assemblage zones. The largest subset of specimens is from the Tropidostoma AZ, with 216 dorsal, 178 dorsal complete, and 208 lateral landmarked specimens. The second largest sample is from the Cistecephalus AZ, with 112 dorsal, 89 dorsal complete, and 102 lateral. From the Tapinocephalus AZ there are 47 dorsal, 37 dorsal complete, and 40 lateral. From the Daptocephalus AZ there are 24 dorsal, 16 dorsal complete, and 22 lateral. Lastly, from the Pristerognathus AZ there are 20 dorsal, 17 dorsal complete, and 19 lateral specimens. Deformation style was determined qualitatively a priori. Many Diictodon skulls are subject to multiple forms of distortion, in which case the dominant style of deformation was given precedence for binning. 

            Procrustes superimposition and principal components analysis (PCA) of landmark data were performed using the program MorphoJ (Klingenberg, 2008). Meaningful PCA axes were determined using the broken-stick method described by Jackson (1993), which distinguishes between eigenvalues providing significant data structure and those that describe random noise. The digitized specimens of Diictodon listed above were also included in a broad-scale PCA covering all of Anomodontia, with the sample composed of 1876 specimens in dorsal view and 1921 specimens in lateral view (including Diictodon). Landmark protocol for this analysis was identical to that for the within-group Diictodon study, with ‘dorsal complete’ landmarking used for skulls in dorsal view. Anomodont specimens were binned into the higher-level taxa described by Kammerer & Angielczyk (2009) for calculation of within-group means when measuring disparity. Procrustes variance-based morphological disparity was calculated for the anomodont-wide datasets using the ‘morphol.disparity’ function of the R (R Core Team, 2018) package geomorph (Adams & Otárola-Castillo, 2013). The significance of pair-wise differences in disparity among groups was assessed via resampling.

Simulations

A series of simulations was conducted using the program DefCat, with similar parameters to the simulation studies of Angielczyk & Sheets (2007). DefCat produces simulated deformed and non-deformed datasets in which both deformation and underlying biological signals are known. Our basic protocol consisted of: 1) using non-deformed empirical specimens near the ends of a biological continuum of variation to generate a series of simulated non-deformed specimens that fall along that continuum; 2) creating deformed datasets by using mathematical transformations to apply known types and amounts of deformation to the datasets of simulated non-deformed specimens; and 3) assessing the amount of variance added to the datasets by deformation and testing whether accurate biological signals could be recovered from the deformed datasets.    

            There were two simulated datasets with known biological signals that were subjected to deformation. The first included a known ontogenetic signal, and the second included a sexual dimorphism signal.

            To make the simulated dataset with an ontogenetic signal, two undistorted Diictodon specimens of different sizes were chosen, one a representative ‘small’ individual and one a representative ‘large’ specimen. SAM-PK-K7838 served as the ‘small’ specimen for all of the simulations because it is especially minute (skull length approximately 2.6 cm, near the lower end of known specimens) and because landmarks could be digitized on it in each view. For the dorsal complete dataset, the ‘large’ specimen was SAM-PK-K6041 (skull length approximately 10.1 cm). The ‘large’ specimen in both the dorsal and lateral datasets was USNM 25157 (skull length approximately 11.6 cm; no complete, undistorted specimens near the maximum size for Diictodon, ~15 cm, are known). DefCat was then used to generate as many simulated specimens as we had empirical specimens. This meant 485 dorsal and 464 lateral simulated specimens generated within each simulation. This evenly divides the size and shape differences between the endpoint specimens, creating a simulated ontogenetic series of undeformed specimens. A small amount of identical independent Gaussian noise was added to the data to simulate individual variation among the specimens.

            The second series of simulated datasets was created to test the effects of deformation on sexual dimorphism. Here, two datasets were created using sets of either likely male or likely female specimens, so that each dataset would be representative of either male or female Diictodon. Following Sullivan et al.’s (2003) hypothesis that tusked individuals of Diictodon are likely male and tuskless individuals are likely female, two undeformed specimens of either sex were used to create a series of simulated specimens that represent only one sex or the other. For the dorsal complete datasets, BP/1/293 and SAM-PK-K1650 were chosen to represent males, and BSPG 1934-VIII-48 and UCMP 42837 were chosen to represent females. For the dorsal view datasets, NHMUK PV OR 47052 (the holotype of Diictodon feliceps) and SAM-PK-K1650 were chosen to represent males, and BSPG 1934-VIII-48 and TM 299 were chosen to represent females. For the lateral view datasets, NHMUK PV OR 47052 and TM 373 were chosen to represent males, and BSPG 1934-VIII-48 and SAM-PK-K11484 were chosen to represent females. In each case, the amount of simulated specimens reflect those in each category for the empirical dataset, resulting in 251 simulated females and 220 males for the dorsal analysis, 194 females and 179 males for the dorsal complete analysis, and 238 females and 214 males for the lateral analysis. 

            For each analysis, 100 simulated specimens were generated that evenly divide shape variation among the starting specimens, with identical independent Gaussian noise added to the data to simulate individual variation.

            To generate our deformed datasets, we used the deformation model described in Angielczyk & Sheets (2007), which applies uniform shear and stretching to a series of landmark configurations. Although shear deformation is only one of several types of deformation observed in specimens of Diictodon (see below), it provides an intuitive starting point for investigating the potential effects of deformation on morphometric data for the species. The deformation model has two main parameters, θ and a. The θ parameter alters the ratio of the long axis to the short axis of the strain ellipse, whereas the a parameter varies the strain ellipse ratio and the orientation of the strain ellipse (details in Angielczyk & Sheets [2007]). In addition to these parameters, the orientation of a specimen’s landmark configuration relative to the strain ellipse will alter the details of its resulting deformation. For example, if the long and short axes of the specimen are aligned with the direction of applied stress, no shearing of the specimen will be apparent, whereas shear will occur when the axes are not aligned. Therefore, an additional parameter of the simulations is the range of angles specimens can take relative to the applied stress, which can be random (i.e., ranging from -180º to 180º) or constrained to a smaller range of angles to produce more stereotyped patterns of deformation in a given dataset. 

            A total of 18 groups of deformed datasets were generated, nine containing a simulated ontogenetic signal and nine with a simulated sexual dimorphism signal. For a given biological signal, each of the nine groups consisted of 16 individual datasets in which θ = 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, and 20.0, with additional parameters varying from group to group. In Group 1, a = 1.0, all specimens were deformed, and the specimens were oriented randomly relative to the direction of applied stress. Groups 2 and 3 had the same parameters as Group 1, except that variable numbers of specimens were left undeformed (50% deformed in Group 2, 94% deformed in Group 3). Groups 4 and 5 had the same parameters as Group 1 except that the orientation of the specimens relative to the strain ellipse was constrained (-45º to 135º in Group 4; 0º to 90º in Group 5). Groups 6 and 7 had the same parameters as Group 1, but deformation amplitude (a) was allowed to vary (0.95 to 1.0 in Group 6; 1.0 to 1.05 in Group 7). All parameters were allowed to vary in Groups 8 and 9.

            We calculated Procrustes variance-based morphological disparity for all of the simulated deformed datasets using the ‘morphol.disparity’ function of the R package geomorph (Adams & Otárola-Castillo, 2019). These values were then compared to disparity values for the corresponding empirical dataset (e.g., simulated deformed lateral view vs. empirical lateral view) to determine which combinations of deformation parameters produced datasets with levels of disparity comparable to our empirical sample. The significance of differences in disparity was assessed using re-sampling. To determine whether an accurate ontogenetic signal was present in our deformed datasets, we compared the size-correlated shape variation in each deformed dataset to its corresponding undeformed dataset using a homogeneity of slopes test and an ontogenetic trajectory analysis (Adams & Collyer, 2009). An accurate ontogenetic trajectory was considered recoverable if the slope and trajectory parameters for the deformed dataset did not differ significantly from the parameters of the corresponding undeformed dataset. Analyses were also carried out the geomorph package in R.

            To test the effects of deformation on the sexual dimorphism signal, we conducted Procrustes ANOVAs in MorphoJ (Klingenberg, 2008) to determine whether simulated male specimens differed significantly in mean shape from simulated females for each set of deformation parameters. Of particular interest in the tests of ontogenetic and sexual dimorphism signals was whether an accurate biological signal could be recovered from the simulated deformed datasets with disparity levels closest to the observed value for the corresponding empirical dataset.