Cheek teeth are filled to the cusp with information about mammalian evolution. Studying the evolution of mammalian cheek tooth crown complexity has benefited our understanding of mammalian evolution in developmental, morphological, and ecological contexts. Most work is focused on individual cheek tooth loci as opposed to considering the premolars and molars as serial homologues. This focus on individual tooth loci has left the exploration of inter-regional phenomena understudied. One such phenomenon is the molarization of premolars across hoofed mammals; some have simple unicuspid premolars while others have premolar crowns that are equal in complexity to their molars. Many developmental models have been proposed to understand cheek tooth evolution, but minimal work has been done to synthesize these models into a holistic understanding of cheek tooth crown complexity evolution. We investigated if applying a synergized theoretical framework of the inhibitory and patterning cascade model to artiodactyl and perissodactyl taxa could be used to study the evolution of molarization in hoofed mammals. We applied an existing 2D landmarking scheme for the upper and lower premolar molar boundaries of hoofed mammals to capture the morphology across this important identity boundary. Shape data were analyzed through phylogenetically informed modularity analyses to capture the covariation structure at the upper and lower premolar-molar boundaries. A-priori modularity hypotheses were proposed based on developmental models including the patterning cascade model and the inhibitory cascade model. Both artiodactyl and perissodactyl results showed support for modularity across the upper and lower premolar molar boundary but showed more variation in the upper premolar molar boundary. Artiodactyls show consistency in support for modularity hypotheses between upper and lower premolar molar boundaries where perissodactyls show significant differences in support for modularity hypotheses between upper and lower premolar molar boundaries. Our results illustrate that the covariation structure at the premolar molar boundary has convergent and divergent elements that both have consequences for our understanding of the evolution of molarization within and between artiodactyls and perissodactyls.

Specimen sampling and imaging

Phylogenetic comparative methods were used to approach our research question (see modularity analyses), so the specimen sampling in our study relied on previously published phylogenies of artiodactyls and perissodactyls (Fraser et al. 2015). These phylogenies are composite trees including extant and extinct taxa for both artiodactyls and perissodactyls (Janis et al. 1998; Price et al. 2005; Maguire and Stigall 2009; Fraser et al. 2015). We modified the tip-dates and time-scaling of these trees to incorporate the most recent updates to first and last occurrence data from the Paleobiology Database (https://paleobiodb.org/#/; (Alroy et al. 2001; Uhen et al. 2013)). Specimens sampled must have been present on these trees for us to include them in our study. Nexus files for these updated trees can be found in our supplemental materials (Supplementary Materials 1).

Several exclusion and inclusion criteria were applied to sample artiodactyls and perissodactyls. Only one representative specimen per species was selected to ensure that the generated morphospace would not be weighted towards species that are overrepresented in museum collections. Only mature specimens with a fully erupted permanent fourth premolar were included. Specimens with minimal occlusal wear on the premolars and molars were prioritized for inclusion to be able to identify cusps and cuspids for landmarking. In addition, specimens had to have four identifiable cusps on the P4/p4 and M1/m1 to ensure landmark standardization. Lastly, the upper or lower premolar molar boundaries must have been in-situ within the maxilla or mandible respectively to minimize error influencing the position of the premolar molar boundary. Our final sample included 45 upper artiodactyl, 18 lower artiodactyl, 30 upper perissodactyl, and 30 lower perissodactyl premolar molar boundaries. Specimens were mostly sampled and imaged from the American Museum of Natural History Zoology and Paleontology collections with additional specimens sourced from our previous work on this topic (Ashbaugh et al. 2025). All specimens used in our study are documented in the supplement (Supplementary Materials 2).

We took occlusal images of each specimen using an IPhone 14 Pro (48MP Main: 24 mm, ƒ/1.78 aperture). Right mandibles with the premolar molar boundaries in focus 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 premolar molar boundary but had an intact left premolar molar boundary that passed all exclusion and inclusion criteria, images were taken of the left premolar molar boundary and then reflected prior to landmarking.

2D GMM, Landmarking, and PCA

The landmark scheme includes 58 landmarks to capture cusp or cuspid covariation across the premolar molar boundary by landmarking the placement of the cusps or cuspids and their respective contours. This landmark scheme has been validated in our previous work (Ashbaugh et al. 2025), but modifications were required to apply the scheme to the upper premolar molar boundary. The landmark scheme was developed for the lower premolar molar boundary with an emphasis on the importance of the trigonid and talonid in the development of the lower cheek teeth (Fig. 4; Ashbaugh et al. 2025). The upper premolars and molars of ungulates have a functional equivalent of the trigonid (trigon) but lack the talonid seen in most ungulate lower molariform cheek teeth. To address this issue, we mirrored the lower landmark scheme to the upper cheek teeth and re-defined landmark placement based on anatomical cusp or cuspid placement to generate the modified upper premolar molar boundary landmark scheme (Fig. 4). This modified landmark scheme retains the same focus as the original – the landmarks are designed to gather information about the functional relationship of the cusps and cuspids across the premolar molar boundary and how they relate to the ICM and PCM. The upper and lower premolar molar boundary landmark scheme with detailed descriptions of landmark placement and visual aids for landmarking can be found in the supplement (Supplementary Materials 3). All specimens were landmarked using the appropriate landmark scheme and then the four shape data sets (upper artiodactyls, lower artiodactyls, upper perissodactyls, lower perissodactyls) were translated, rotated, and scaled to the same size by applying a Procrustes transformation to minimize the impact of size differences on shape differences (Zelditch et al. 2012).

Our sample of species included selenodont, lophodont, and equid or equid-like cusp/cuspid morphology. We ran principal component analyses (PCA) as an exploratory approach to study the principal axes of cusp shape variation across the upper and lower premolar molar boundary in our samples. We found clear separation among the various cheek tooth morphotypes, demonstrating that our landmark scheme could accurately capture the premolar molar boundary for all boundary morphotypes studied (Supplementary Materials 4).

Modularity analyses

 The modularity hypotheses tested across the premolar molar boundaries were similar to those proposed in our previous work to ensure consistency across our studies (Ashbaugh et al. 2025). We will briefly review the tested modularity hypotheses and highlight any modifications required to consider both the upper and lower premolar molar boundary.

Our first hypothesis (Validation hypothesis) was proposed to ensure the placement of the cusp or cuspid apex is associated with the maximum buccal or lingual contour of each cusp or cuspid – strong association between these landmarks shows evidence of the patterning cascade model operating from cusp apex to base of the cusp (Ashbaugh et al. 2025; Fig. 5A).

The Mesial and Distal hypotheses capture the covariance of the anterior and posterior cusps and cuspids with the first molar morphology (Ashbaugh et al. 2025; Fig. 5B and 5C). Previously, these hypotheses were framed as the Trigonid and Talonid modularity hypotheses respectively to capture the two regions of the tribosphenic molar in occlusion (Ashbaugh et al. 2025). To collect biologically meaningful shape data between the upper and lower premolar molar boundary, these hypotheses required revisions. The occlusal interactions orchestrated by the mesial and distal components of the cheek teeth largely determine whether an ungulate is capable of grazing or browsing, raising the question of how the mesial and distal components covary across the upper and lower premolar molar boundary (Schultz et al. 2021; Avedik et al. 2023). The biological context of these hypotheses remains the same in the present study – the mesial and distal components of the upper and lower cheek teeth have different roles in occlusion which could be the focus of selective pressure in the context of the PCM (DeMers and Hunter 2024; Hunter et al. 2024; Stenberg et al. 2024). The Buccal and Lingual modularity hypotheses are unchanged from our previous work (Fig. 5C and 5D; Ashbaugh et al. 2025). The buccal and lingual cusps and cuspids collectively impact the chewing cycle in mammalian mastication, but are of particular importance in the transverse nature of ungulate mastication (Butler 1939, 1956, 1972; Couzens et al. 2021; Grossnickle et al. 2022). The importance of the buccal and lingual cusps and cuspids as a unit in transverse mastication warrants the study of these cusps and cuspids in the premolars and how they may covary with molar morphology within the context of the PCM.. The Mesial, Distal, Buccal, and Lingual modularity hypothesis are proposed to be nested within the PCM as it relates to covariation across different dental regions (Ashbaugh et al. 2025). Collectively, these hypotheses serve to identify if the PCM has potential as a guiding developmental framework for understanding how inter-regional interactions may influence how molarization evolved in hoofed mammals. We therefore predict that all four of these modularity hypotheses should be significantly supported relative to the null of no modularity across the premolar molar boundary.

The Independence modularity hypothesis is also unchanged from our previous work (Fig. 5E; Ashbaugh et al. 2025). If independent patterning cascades are a contributor to the evolution of crown complexity evolution seen between the premolars and molars in mammals, we predict that both artiodactyls and perissodactyls should show significant support in the performance of the Independence hypothesis. However, we also predict that there should be a relative difference between artiodactyls and perissodactyls in the performance of the Independence modularity hypothesis. Perissodactyl molar ratios fail to be explained by the ICM while most artiodactyls have molar ratios that align with the ICM (Polly 2007; Halliday and Goswami 2013). We predict that the conformity to the ICM should manifest in our modularity results – artiodactyl modularity results should show more support for the Independence hypothesis in comparison to perissodactyls. The Independence hypothesis provides a pathway to evaluating the ICM within the framework of the PCM in two ways not attempted before - (1) between the cheek teeth and (2) between the upper and lower cheek teeth.

All modularity hypotheses were tested for the four shape data sets independently. If a significant phylogenetic signal was found in 2D geometric morphometric data using the composite phylogenies mentioned previously (evaluated using K_mult; Adams 2014), modularity analyses were conducted with phylogenetic comparative methods (Klingenberg and Marugán-Lobón 2013). Modularity was quantified using a z-transformed covariance ratio (Z-CR) test statistic using the function geomorph::phylo.modularity()(Adams and Collyer 2019). Further detail on z-score calculation and comparison for modularity hypotheses can be found in Adams and Collyer (2019).

Modularity test results were interpreted in two ways: (1) modularity strength and (2) modularity range. If a modularity hypothesis is stronger or weaker compared to other modularity hypotheses proposed it should show significantly different covariation structure relative to the other modularity hypotheses tested, not just the null hypothesis – we refer to this metric as modularity strength in our results and throughout the remainder of the manuscript. To evaluate this we calculated permuted effect sizes of the Z-CR modularity test statistic for the modularity hypotheses, which were subsequently contrasted in pair-wise comparisons using the function geomorph::compare.pls().

Modularity range is a metric that captures variation between the upper and lower premolar molar boundary. We collected and analyzed the range of effect sizes across all modularity hypotheses studied for both the upper and lower premolar molar boundary. We define modularity range as the breadth of modularity strength effect size across all modularity hypotheses when compared to a null hypothesis of no modularity. Comparing the modularity range of the upper and lower premolar molar boundary allowed study of the covariation structure between the upper and lower premolar molar boundary.

All code was developed using R version 4.4.2with the package Geomorph version 4.0.9(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.