Hoofed mammal premolars show a range of occlusal crown morphology from molariform to caniform – the placement of taxa on this spectrum can be described as the relative molarization of the premolars. Molarized premolars function together with the molars in grinding mastication in most of the lineages in which these unique premolars appear. The degree of molarization varies across dietary ecologies, which has led to cheek tooth morphology being designated as an important contributor to dietary predictions in extant and extinct taxa. Recent research into mammalian occlusal cheek tooth patterning have found independent patterning mechanisms of the premolars and molars. A research gap exists in understanding how molarization of the premolars has occurred so frequently in hoofed mammals if these dental regions are independent in their patterning. In this study, we tested the application of the patterning cascade model to the lower premolar-molar boundary in a geometric morphometrics framework in hoofed mammals. We used 2D geometric morphometrics to study occlusal cuspid covariation at the lower p4-m1 boundaries of 16 artiodactyl and 18 perissodactyl species. Phylogenetically informed modularity analyses were used to test alternate a priori hypotheses originating from evolutionary, developmental, and functional considerations of cheek tooth morphogenesis. Our results showed artiodactyls and perissodactyls differ significantly in their p4-m1 boundary covariation patterns, which we hypothesize could be caused by heterochronic shifts between premolar and molar development. To our knowledge, our study is the first to contribute a comprehensive yet accesible 2D geometric morphometric method to further investigate the evolution of molarized premolars.

Specimen selection

We took occlusal images of the lower p4-m1 boundary of extinct and extant hoofed mammals to collect shape data. Right p4-m1 boundaries were placed on an imaging surface with the occlusal surface perpendicular to the line of sight of the lens to minimize parallax. If a specimen lacked the right p4-m1 boundary but had an intact left p4-m1 boundary, specimens were imaged and then reflected for landmarking. Images were taken using a Nikon D200 camera with either an AF-S Nikkor 17-55 mm 1:2.8 G lens (large specimens) or an AF-S Micro Nikkor 105 mm 1:2.8 G lens (small specimens). Specimens were sampled from the University of Calgary Museum of Zoology, the University of Nebraska State Museum Zoology and Vertebrate Palaeontology collections, and the Canadian Museum of Nature Zoology and Vertebrate Palaeontology collections. Additional specimens were imaged by colleagues (Redpath Museum; McGill University, New Brunswick Museum, University of Alaska Museum; University of Alaska, and Washington State Museum of Natural History and Culture; University of Washington).

2D geometric morphometrics is subject to statistical artifacts if the relationship between landmark placement and biological variation in the landmark scheme is not carefully evaluated (Buser et al. 2018). Specimens were chosen as representatives for their corresponding species based on minimizing wear and damage in the occlusal surface – a single specimen represented each species in our analysis (Supplement 1). Specimens included had the lower fourth premolar and first molar fully erupted in-situ with no visible damage to the cuspids in occlusal view (See 2D GMM section). Tooth wear can also change the occlusal surface of the tooth in most mammals making it difficult to landmark important cuspids (Galbany et al. 2011; Famoso et al. 2013; Schultz et al. 2021). To account for this, we developed a method to quantify cuspid wear (cuspid wear coefficient) derived from a previous study on cusp wear in Amboseli baboons (Galbany et al. 2011). Cuspid wear coefficients were used to exclude specimens that were too worn to landmark effectively (Supplement 1). An additional criterion was that species selected must have at least two trigonid cusps and two talonid cusps to be landmarked. Species that had less than four cusps on the p4 or m1 were not included in our study.

Specimens that passed inclusion criteria and cuspid wear quantification were cross-referenced with conservative composite meta-analysis phylogenetic trees based on a variety of work that has contributed to our understanding of artiodactyl and perissodactyl relationships (Janis et al. 1998; Price et al. 2005; Maguire and Stigall 2009; Fraser et al. 2015). While we understand that these phylogenetic trees may not be accurately tip dated at this point, they preserve  the key relationships among ungulate taxa studied in relation to node placement. Our focus in this study is understanding the changes in molarization that occur through the tree (node focused) rather than understanding the evolution of particular modern molarized premolars (tip-focused). We encourage future workers to adapt this framework with new and updated trees as they are published. Phylogenetic information was used to calculate phylogenetic signal in shape data collected in addition to permitting the use of phylogenetic comparative methods in our study (See modularity analyses). Species that were not present in this phylogeny were excluded from our study. The Nexus files and phylogenies used in our study can be found in the supplementary materials (Supplement 2). Overall, the specimens presented in the following results consist of 34 species of extinct and extant hoofed mammals including 16 terrestrial artiodactyls and 18 perissodactyls (listed in Supplement 3)

2D Geometric Morphometrics

The landmark scheme we developed captures the cuspids and their manner of growth that are comparable between the p4 and m1. The four principal cuspids of  lower cheek teeth that have the longest time to develop within the premolar and molar specific PCM are our primary cuspids of interest. The landmark scheme includes 58 landmarks to capture cuspid covariation by landmarking the placement of the cuspids and their respective contours for the p4 and m1 together (Fig. 4). Thirty-two of these landmarks were assigned as sliding semi-landmarks (SSLM) to capture the contour of each major cuspid. A detailed description of our landmark scheme and the classification given for each landmark can be found in the supplementary materials (Supplement 4).

Once landmarked, shape data were extracted from occlusal images. These shape data were separated into two data sets – one for artiodactyls and one for perissodactyls. For each of these data sets the landmark configurations were scaled, translated, and rotated to the same position and relative size using Procrustes superimposition to minimize the impact of size on shape differences (Zelditch et al. 2012). SSLM were slid to minimize Procrustes distance. Recent work has shown prioritizing bending energy as the sliding criterion can generate spurious integration between modules that are not truly biologically related (Zelditch and Swiderski 2023). Because our work is focused on differentiating modules, we chose Procrustes distance as our sliding criterion.

Our sample of ungulate species covers a wide breadth of cheek tooth morphologies including variants of selenodonty, lophodonty, and equid-like dentiiton. To ensure that our landmark scheme is capable of capturing shape data that accounts for this morphological diversity, we ran two principal components analyses (PCA) to visualize shape data collected; one PCA for all artiodactyls and another for all perissodactyls. Ordination results showed clear separation among the different cheek tooth morphotypes validating that our landmark scheme can accurately describe the wide variety of cheek tooth morphotypes seen in our specimen pool (Supplement 5)

We chose not to analyze allometry in these shape data considering our focus on phylogenetic representation rather than precision. Assuming that all species studied share a common impact of allometry on occlusal dental complexity would be inappropriate. Allometry has played a role in diversification in many lineages of hoofed mammals and would require replicates of individual species to evaluate this effect with biological significance (Illius and Gordon 1992; Ungar 2014; Klingenberg 2016). This study does not use the framework necessary to detect differences in allometric components among species studied (Mitterocker et al. 2013).

Landmarking was conducted using the TPSSuite including TPSDig (Version 2.31; Rohlf 2013a) and TPSUtil (Version 1.78; Rohlf 2013b).

Modularity analyses

Modularity was introduced as a systems-based approach to understand biological variation (Olson & Miller 1958), but has recently been applied to multivariate shape spaces (morphospaces) to test a-priori mechanistic hypotheses of morphological evolution (Polly 2008; Polly et al. 2016; Adams and Collyer 2019). Given it is impossible to directly study the development of extinct artiodactyls and perissodactyls and generate modularity hypotheses based on observation as is preferred by this method (Mitteroecker and Bookstein 2007; Klingenberg 2008; Zelditch and Goswami 2021), we chose to evaluate interspecific cuspid shape covariances as a proxy for developmental modularity of the PCM across the the p4-m1 boundary. We proposed six different modularity hypotheses. We used this approach to explore how these modularity hypotheses performed when compared to a null modularity hypothesis and to one another, to better understand how inter-regional shape data may elucidate the evolution of molarization in hoofed mammals (Fig. 5). Our first hypothesis (Validation hypothesis) was proposed to ensure the placement of the cuspid apex is associated with the maximum buccal or lingual contour of each cuspid (Fig. 5a). If we fail to reject this hypothesis, it would show that our landmark scheme can capture cuspid-specific modularity as patterned by the PCM. The remaining hypotheses are used to test whether modularity is present between the p4 and m1 and if so, which developmental processes might have contributed to organizing cuspid shape covariation across the boundary.

The Trigonid modularity hypothesis represents the functional importance of trigonid covariation between the p4 and m1 cuspids (Fig. 5b). The trigonid is strongly associated with shearing facets which are important in herbivorous diets that tend toward browsing or selective feeding (Crompton 1971; Crompton and Jenkins 1979; Davis 2011). The Talonid modularity hypothesis represents functional covariation of the talonid cuspids across the p4-m1 boundary (Fig. 5c). Talonid cuspids on the p4 are those most proximal to the m1 and occlude between the p4 and m1 with the upper dentition in herbivorous mammals (Crompton 1971; Crompton and Jenkins 1979; Davis 2011). This occlusion participates in grinding mastication against the corresponding upper teeth which is vital in grazing or bulk-feeding herbivore mastication (Schultz et al. 2021; Avedik et al. 2023).

The Buccal and Lingual modularity hypotheses emerge from both classic and recent work investigating the contribution of the buccal and lingual cusps to mammalian mastication (Butler 1952, 1956, 1972; Couzens et al. 2021). Simulations support the presence of a buccal module citing the importance of the buccal region in generating future cusps (Couzens et al. 2021). The upper buccal cusps are also suggested to have a large role in the diversification of tooth morphology leading to unilateral hypsodonty in early ungulates. To date there is limited discussion of the existence and/or function of a buccal or lingual module in the lower dentition. We proposed the Buccal and Lingualmodularity hypothesis to explore the strength of any cross-boundary signal on these function-driven hypotheses in the lower dentition. The Buccal modularity hypothesis represents the functional covariation between the molar cuspids and the premolar cuspids that are buccal to the midline of the tooth  (Fig. 5d) while the Lingual modularity hypothesis represents the functional covariation between the molar cuspids and the premolar cuspids that are lingual to the midline of the tooth (Fig. 5e).  

Lastly, we tested an Independence hypothesis that the p4 and m1 are patterned independently and therefore each have their own cuspid covariation structure (Fig. 5f), based on literature demonstrating that independent patterning mechanisms exist for the premolars and molars (Kavanagh et al. 2007; Sadier et al. 2023).

All modularity hypotheses were tested for artiodactyls and perissodactyls independently. If a significant phylogenetic signal was found in 2D morphometric data using the composite meta-analysis phylogenies mentioned previously (evaluated using Kmult; Adams 2014), modularity was evaluated with phylogenetic comparative methods (Klingenberg and Marugán-Lobón 2013). Modularity was quantified using the covariance ratio (CR) test statistic using the geomorph::phylo.modularity() function (Adams and Collyer 2019).

Calculated effect sizes for the various hypotheses were contrasted in pair-wise comparisons using the function geomorph::compare.pls() (Adams and Collyer 2019). Further detail on z-score calculation and comparison for modularity hypotheses can be found in Adams and Collyer (2019).

All code was developed using R version 4.1.2 with the package Geomorph version 4.0.3 (Adams and Otárola-Castillo 2013; R Core Team 2020). All shape files and R code used to complete these analyses are accessioned in this data dryad repository.