The Impact of Hybridization on Upper First Molar Shape in Robust Capuchins (Sapajus nigritus x S. libidinosus)

Hybridization, once viewed as rare and universally detrimental (Dobzhansky, 1940; Mayr, 1963), is increasingly viewed as a frequent and innovative evolutionary phenomenon (Ackermann et al., 2019; Arnold, 1997; Taylor & Larson, 2019). Adaptive hybridization and fertile hybrid populations have been observed in a wide array of animals and plants, such as Galápagos finches (Grant & Grant, 2020), toads (C. Chen & Pfennig, 2020), butterflies (Jiggins et al., 2008), and poplar trees (Chhatre et al., 2018). There is an especially rich body of literature on ancient and contemporary primate hybridization. Genetic analyses have demonstrated that hybridization events occurred in many primate lineages during their course of evolution (de Manuel et al., 2016; Fan et al., 2018; Kuhlwilm et al., 2019; Svardal et al., 2017; Tung & Barreiro, 2017; Zichello, 2018), including hominins (Browning et al., 2018; L. Chen et al., 2020; Durvasula & Sankararaman, 2020; Green et al., 2010; Huerta-Sánchez et al., 2014; Reich et al., 2011; Slon et al., 2018). A growing number of primate taxa are proposed to have hybrid origins (Burrell et al., 2009; Detwiler, 2019; Rogers et al., 2019; Roos et al., 2019; Thinh et al., 2010; Tosi et al., 2000; Wang et al., 2015). Hybridization continues to shape genetic and phenotypic variation in present-day primate populations in both natural and anthropogenic contexts (Alberts & Altmann, 2001; Bergman et al., 2008; E. L. Bynum et al., 1997; Cortés-Ortiz et al., 2007; Gligor et al., 2009; Jolly et al., 2011; Malukiewicz et al., 2015; Mather, 1992). While genetic analysis has been crucial for exploring primate hybridization, there is a growing interest in understanding the impact of hybridization on morphology. Studies of hybrid primate morphology offer unique insight into the effect of ABSTRACT To better understand the impact of hybridization on development and morphology, I analyze an understudied phenotype in hybrid morphology research: tooth shape. I apply a 2D geometric morphometric approach to compare variation in first upper molar cusp tip positions and crown outline shape among 31 crested capuchins (Sapajus nigritus), 37 bearded capuchins (S. libidinosus), and 44 hybrids (S. nigritus x S. libidinosus). A principal components analysis shows that group membership accounts for a significantly greater proportion of variance along the first major axis of M1 shape variation than does allometry. While most hybrids have S. nigritus-like M1s, several possess a transgressive M1 shape not observed in either parental species. Procrustes distances are greater in hybrids compared to the parental capuchins, and two-block partial least squares analyses show that hybrids exhibit weaker integration between cusp tip positions and crown outline shape. These results demonstrate that hybridization generates novel M1 shapes and support the hypothesis that destabilized development results in elevated phenotypic variance in hybrids. Further studies of dental shape in hybrid primates will generate important data for on-going efforts to detect potential hybrids in the hominin fossil record and to understand the evolutionary outcomes of anthropogenic hybridization.

While genetic analysis has been crucial for exploring primate hybridization, there is a growing interest in understanding the impact of hybridization on morphology. Studies of hybrid primate morphology offer unique insight into the effect of rapid genetic recombination on phenotypic development and variation. Many studies have evaluated soft tissue phenotypes in contemporary hybrid populations using traits that are easy to observe in the field and can be measured non-invasively, such as pelage color and distribution, head shape, and tail carriage (Alberts & Altmann, 2001;N. Bynum, 2002;Kelaita & Cortés-Ortiz, 2013;Phillips-Conroy & Jolly, 1986). Researchers have also studied the relationship between hybrid ancestry and hard tissue phenotype, such as skeletodental size, shape, and non-metric trait variation (Ackermann et al., 2006;Ackermann & Bishop, 2010;Boel, 2016;Cheverud et al., 1993;Eichel & Ackermann, 2016;Ito et al., 2015;Kohn et al., 2001;Phillips-Conroy, 1978). The proximate aims of hybrid morphology research are to elucidate how hybrid morphology quantitatively and qualitatively differs from parental morphology and if different kinds of hybrids share diagnosable traits indicative of their hybrid ancestry (Ackermann, 2010;Ackermann et al., 2019).
The data derived from hybrid morphology research has important broader implications for primate conservation and paleoanthropology. Primate conservation biologists observe that the frequency of hybridization will likely increase as primate habitats are disturbed or destroyed by anthropogenic interference (Detwiler et al., 2005; Malukiewicz, 2019; Thompson et al., 2018). Rare, endangered primates may reproduce with more common heterospecifics if conspecific mates are difficult to find. Extensive admixture between divergent taxa may result in loss of genetic and phenotypic diversity and ultimately fuse two lineages (Seehausen et al., 2008), or it may generate novel diversity and prevent inbreeding depression (Arnold & Meyer, 2006). By studying the variation in hybrid phenotypes, conservation biologists may be able to understand if the outcomes of anthropogenic hybridization are harmful, neutral, or adaptive for endangered primate populations.
Paleoanthropologists acknowledge that hybrid hominins are likely present in the fossil record. Fossil evidence demonstrates that multiple hominin taxa cohabited Africa and Eurasia throughout the Pliocene and Pleistocene and could have hybridized where their ranges overlapped (Détroit et al., 2019;Grün et al., 2020;Herries et al., 2020;Spoor et al., 2015). The genetic evidence for hybridization events throughout hominin evolution is substantial (Durvasula & Sankararaman, 2020;Jacobs et al., 2019;Sankararaman et al., 2016;Skov et al., 2020;Villanea & Schraiber, 2019). The hybrid ancestry of several fossilized hominin individuals has been confirmed by ancient DNA analyses (Fu et al., 2015;Slon et al., 2018). However, ancient DNA preservation is rare in most of the hominin fossil record, so analyses of hard tissue phenotypes in extant hybrid primates can be used to assess the feasibility of using morphological indicators to identify hybrid hominin fossils (Ackermann et al., 2019). The identification of hybrid hominin fossils remains an outstanding issue for reconstructing hominin phylogenetic relationships, as most the commonly used phylogenetic frameworks assume evolutionary relationships are hierarchical rather than reticulate (Holliday, 2003).
Quantitative genetic theory states that in firstgeneration (F1) hybrids, phenotypic trait measurements controlled by additive genetic variation will be the midparental value (MPV), or the averaged parental measurements (Falconer & Mackay, 1997). Tests of this theory indicate that while some F1 hybrid primate phenotypes exhibit the expected MPV (Hamada et al., 2012), other traits in the same population may deviate from the expected phenotype (Ackermann et al., 2006;Cheverud et al., 1993;Eichel & Ackermann, 2016). Positive deviations from the MPV in F1 populations is referred to as heterosis, or hybrid vigor, while negative deviations are evidence of dysgenesis, or hybrid breakdown. Later-generation hybrids with higher genetic input from one parental taxon are expected to be more phenotypically like that parent, but some hybrids resemble one parent more than the other, regardless of parental genetic contribution (Boel et al., 2019;Ito et al., 2015). First-and later-generation hybrid populations sometimes exhibit transgressive phenotypes not observed in either parental population, such as extreme trait size, novel combinations of parental traits, or the presence of nonmetric craniodental anomalies Ackermann & Bishop, 2010;Jolly et al., 1997). Research on hybrid morphology in primates has documented a complex array of phenotypic outcomes that vary within and among hybrid populations (Alberts & Altmann, 2001). Importantly, these outcomes are not universally maladaptive (Charpentier et al., 2012) and may help hybrid populations occupy ecological niches unavailable to either parental population, thereby resulting in novel evolutionary lineages (Arnold, 1997;Zinner et al., 2011).
The high morphological variability observed within and among hybrid populations is thought to be the result of destabilized development (Clarke, 1993). The uniquely adapted developmen-tal regimes of two distinct parental taxa are unlikely to merge seamlessly in offspring and could result in perturbations during hybrid morphogenesis. This is supported by the observation that deviations from predicted F1 midparental phenotypes tend to be more pronounced with increasing genetic distance between parental populations (Bernardes et al., 2017;Z. J. Chen, 2013;Stelkens & Seehausen, 2009). Researchers have tested the hypothesis that hybrids experience destabilized development using tests of morphological integration and fluctuating asymmetry (Alibert et al., 1994;Jackson, 1973;Klingenberg, 2003;Klingenberg & McIntyre, 1998). Tightly integrated trait complexes and highly symmetric bilateral trait measurements are hypothesized to reflect stable, canalized development. So, if hybridization results in developmental destabilization, hybrids are expected to exhibit weaker trait integration and greater fluctuating asymmetry between bilateral traits than parental taxa. Some hybrids do meet these expectations Leary et al., 1985;Neff & Smith, 1979), but others do not differ from observed levels of parental trait integration or fluctuating asymmetry (Jackson, 1973;Pallares et al., 2016). In some cases, hybrid samples exhibit stronger trait integration and bilateral trait symmetry than parents, indicating that hybrid development is more stable than parental development (Alibert et al., 1994;Boel et al., 2019;Debat et al., 2000).
Despite growing interest in primate hybrid morphology, the relationships among hybrid ancestry, development, and phenotype remain unclear and difficult to predict. However, one of the most potentially informative anatomical regions for this research has also been one of the most understudied: the dentition. Several lines of evidence suggest that in-depth analyses of dental phenotypic variation will produce valuable data for hybrid morphology research. Anomalous dental non-metric traits are observed at high frequencies in some hybrid populations, such as supernumerary teeth, crown rotation and/or malformation, and dental crowding Ackermann & Bishop, 2010;Goodwin, 1998;Heide-Jorgensen & Reeves, 1993). Intergeneric hybrids of Theropithecus gelada and Papio hamadryas ("geboon") exhibit combinations of parental traits in their dentitions, such as T. gelada-like enamel crenulation on P. hamadryas-like low-crowned molars, resulting in novel dental phenotypes (Jolly et al., 1997). Most of the geboon hybrids also exhibited maxillary cheektooth dimensions that exceeded the parental means. However, hybrids of more closely related baboon species P. hamadryas and P. anubis were not easily differentiable from parental species based on both metric and non-metric dental traits (Phillips-Conroy, 1978). Similarly, dental non-metric trait expression did not discriminate between closely related Macaca fuscata, M. cyclopis, and their hybrids (Boel et al., 2019). Further analyses of dental size, shape, and non-metric trait expression in extant primate hybrids would elucidate if hybrid primates exhibit shared patterns of dental trait variation.
Dental phenotypic analyses of hybrids also could help to understand if deviations from typical parental development generate the high variability observed in hybrid populations. Mammalian dental development is well-studied, and models of dental development have been tested in both extinct and extant primates (Evans et al., 2016;Hlusko et al., 2016;Jernvall & Jung, 2000;Ortiz et al., 2018;Paul et al., 2017). The iterative nature of dental development results in predictable patterns of dental trait integration both within the same tooth crown and among metameres. The patterning cascade model claims that the duration of tooth germ growth and the spatiotemporal distribution and strength of embryonic signaling centers within the germ constrain possible cusp configurations and crown size in the fully formed tooth (Jernvall, 2000). So, differences between parental and hybrid cusp configurations and accessory cusp expression likely reflect deviations in underlying patterning cascade pathways. Similarly, the inhibitory cascade model states that mammalian mandibular molar number and relative size are dictated by embryonic signaling strength and duration of odontogenesis, so differences between hybrid and parental molar size relationships and molar number likely reflect differences in this developmental pathway as well (Kavanagh et al., 2007). Indeed, in a hybrid baboon population, supernumerary mandibular molars are positively correlated with increased molar row length, which suggests that dental development is prolonged in the hybrids compared to parents .
Data derived from studies of hybrid dentitions is especially useful for conservation biologists and paleoanthropologists. Results derived from studies of hybrid skulls and postcrania are not easily applied in living primate populations, but the teeth of primates in hybrid zones can be evaluated, photographed, or molded and cast during trapping expeditions (Kelaita & Cortés-Ortiz, 2013;Phillips-Conroy, 1978). While skeletal data is certainly use-ful for paleoanthropologists interested in determining the feasibility of identifying hybrid ancestry using fossil morphology, teeth tend to be better preserved and comprise most of the hominin fossil record (Bailey, 2002;Gómez-Robles et al., 2007;Martinón-Torres et al., 2012;Wood & Abbott, 1983).
The genus Sapajus is an excellent study taxon for hybridization research. The robust capuchin clade underwent rapid radiation and expansion during the Pleistocene, and species often interbreed where their ranges meet (Lima et al., 2018;Lynch Alfaro, Boubli, et al., 2012), making them an appropriate analog for understanding hominin hybridization. A sample of hybrids of Sapajus nigritus and S. libidinosus are housed at the Smithsonian National Museum of Natural History (NMNH). Sapajus nigritus and S. libidinosus shared a common ancestor approximately 2.6 Ma and belong to different clades within the genus, the former belonging to a more ancient clade endemic to the Atlantic Forest of Brazil, and the latter belonging to a recently evolved clade adapted to Brazilian dry shrublands (Lima et al., 2018;Wright et al., 2015). Both species are listed as 'near threatened' by the IUCN Red List of Threatened Species and both are known to occupy habitats disturbed by agricultural practices Melo, Fialho, et al., 2015). It is possible that anthropogenic hybridization could result in the loss of genetic and phenotypic diversity among robust capuchin species (Lynch Alfaro et al., 2014;Martins et al., 2017). A morphological analysis of hybrid robust capuchins would establish if phenotypic diversity is impacted by hybridization.
Here, I apply 2D geometric morphometric (2DGM) techniques to study variation in first upper molar (M 1 ) crown outline shape and cusp tip configuration among Sapajus nigritus, S. libidinosus, and their hybrids. Dental shape has been used to study population affinity and to characterize extinct and extant primate taxa (Bailey et al., 2016;Gamarra et al., 2016;Gómez-Robles et al., 2007Rizk et al., 2013), including robust capuchins (Delgado et al., 2015), but has not yet been used to study patterns of morphological variation among hybrids and their parental taxa. The primary aims of this study are to explore variation and the factors driving variation in M 1 morphology in hybrids compared to S. nigritus and S. libidinosus; to determine if M 1 morphology can discriminate between hybrids and parental taxa; and to evaluate if hybrids exhibit evidence of destabilized dental development compared to parental taxa. Based on previ-ous hybrid morphology research, I tested the following predictions: 1) M 1 shape is statistically distinct among S. nigritus, S. libidinosus, and their hybrids. 2) The mean shape of hybrid M 1 s is the midparental value (the mean shape of the combined parental sample). 3) There is more variability in M 1 shape within the hybrid sample than within either parental sample. 4) Hybrids exhibit weaker covariation between cusp tip configuration and crown outline shape than parental taxa.

Materials and Methods
My sample includes Sapajus nigritus (n = 31), S. libidinosus (n = 37), and a hybrid sample of S. nigritus x S. libidinosus (n = 44). The dental sample comprises 112 right M 1 s (Table 1). Only specimens with unworn or minimally worn M 1 s were included.
I used a Nikon D500 DSLR digital camera fitted with a macro lens and attached to a copy stand to photograph M 1 occlusal surfaces. I positioned the M 1 cementoenamel junction (CEJ) parallel to the lens and included a scale placed at the same level as the occlusal plane (Bailey, 2004;Gómez-Robles et al., 2007). While directly referencing the specimen, I marked each of the four main M 1 cusp tips (the paracone, protocone, metacone, and hypocone) on the digital images using the GNU Image Manipulation Program version 2.10.12 (The GIMP Development Team, 2019). If a specimen exhibited slight wear, I marked the cusp tip in the center of the wear facet.
I uploaded the photographs to TpsDig2 version 2.31 to digitize a series of 2D landmarks (points of biological homology among specimens) and semilandmarks (non-homologous points of morphological interest; Bookstein, 1997). Landmarks 1 through 4 were placed on the tips of the paracone, protocone, metacone, and hypocone ( Figure 1). These landmarks capture variation in the position of the main cusps relative to each other and relative to the crown outline. In order to examine variation in the shape of M 1 crown outlines, I placed 30 semilandmarks around the perimeter of the occlusal surface, starting at the point of maximum curvature where the buccal and mesial margins intersect. I drew a closed curve around the crown outline, and then appended 29 additional equidistant semilandmarks to the curve (see Figure 1). Finally, I exported all landmark and semilandmark coordinates as a .tps file to RStudio version 1.2.5033 for analysis.
All geometric morphometric analyses were performed using the R package geomorph . First, I defined semilandmarks 5 through 34 as sliding semilandmarks. Sliding semilandmarks can move along the crown outline between neighboring semilandmarks to optimize their position with respect to the average shape of the entire sample. This process removes random variation from the coordinate data introduced by the initial arbitrary placement of semilandmarks around the crown margin and converts semilandmarks 5 through 34 to homologous points statistically comparable to landmarks 1 through 4 (Gunz & Mitteroecker, 2013). Next, I performed a generalized Procrustes analysis (GPA) on the landmark and sliding semilandmark coordinates to remove the effects of specimen size, orientation, and position, leaving only variation related to shape. The GPA superimposes specimens by translating, scaling, and rotating the coordinates to generate an average shape, or consensus configuration, for the entire sample (Bookstein, 1997;Zelditch et al., 2012). I used the results of the GPA for all subsequent analyses.
To explore and compare major axes of M 1 shape variation among S. nigritus, S. libidinosus, and their hybrids, I conducted a principal components analysis (PCA). To visualize variation in shape space, I plotted PC 1 against PC 2 and PC 2 against PC 3. I included 95% confidence ellipses for each taxon to illustrate within-taxon variability. Then, to test the effect of allometry on M 1 shape, I constructed twelve linear models using results from the PCA (Table 2a). Each model tested the association between scores on PCs 1, 2, and 3 with taxonomic designation, logarithm-transformed centroid size, or a combination of both variables. The best-fitting model for PCs 1, 2, and 3 were selected using the function for Akaike's information criterion (AIC) in the R package bbmle (Bolker et al., 2020).
I used warp grids representing the mean shape for each taxon to visually evaluate if and how M 1 shape varies among groups. To statistically evaluate the extent to which M 1 shape is morphologically distinct among these groups by maximizing intergroup differences, I extracted the first ten PCs derived from the PCA (encompassing the majority of shape variation within the sample) for a discriminate function analysis (DFA). Then I used the results from the DFA for a cross-validated assignment test.
I measured within-and between-group variance using pairwise Procrustes distances (the Euclidean distance between two sets of shape coordinates; Spoor et al., 2015). A Procrustes distance equal to zero represents a pair of individuals with identical M 1 shape, while increasing distance reflects increasing dissimilarity in shape. I evaluated statistical differences in Procrustes distances among taxa using pairwise t-tests using Bonferroni correction for multiple comparisons.
I performed a two-block partial least squares analysis (2B PLS) to evaluate the level of covariation between the position of cusp tips (block 1: landmarks 1 through 4) and the shape of the crown outline (block 2: sliding semilandmarks 5 through 34), and implemented a permutation procedure (n = 1,000 permutations) to test the r-PLS correlation coefficients generated by the 2B PLS for statistical significance. Because calculation of the r-PLS statistic is dependent on sample size, I employed a standardized z-score converted to pairwise effect sizes to compare the strength of integration among groups (Adams & Collyer, 2016). Large pairwise effect sizes indicate that the level of morphological integration differs between the two samples.

Mean M 1 shape in parental and hybrid taxa
The mean shapes for S. nigritus, S. libidinosus, and the hybrids compared to the pooled-sample consensus configuration are shown in Figures 2a, 2b, and 2c, respectively. Differences in mean shape are magnified by a factor of three to assist in visual interpretation. The average crown outline shape in S. nigritus is rhomboid, while that of S. libidinosus is more ovoid. The average crown outline in the hybrid sample is more mesiobuccally skewed than either parental taxon and has a waisted lingual margin. The two parental taxa exhibit similar intercusp distances relative to the crown outline, but the mean S. nigritus paracone, metacone and hypocone (landmarks 1, 3, and 4, see Figure 1) are buccally displaced compared to the consensus cusp tips. The average hybrid protocone, paracone, and metacone (landmarks 1 through 3, see Figure 1) are slightly mesially displaced compared to the consensus configuration. The average hybrid M 1 shape differs from the expected midparental shape. The midparental M 1 crown outline does not have the waisted lingual margin that is present in the hybrid mean outline, and the hybrid protocone, paracone, and metacone are mesially displaced compared to the expected midparental M 1 cusp configuration.

Principal components analysis
The first three PCs account for approximately half (48.1%) of the variation in M 1 shape among S. nigritus, S. libidinosus, and S. nigritus x S. libidinosus. Principal component 1 explains 20.8% of shape variation, while PC 2 explains 16.6% (Figure 3a). The warp grids representing M 1 shape at extreme ends of variation along each PC illustrate that M 1 s with low PC 1 scores have a mesiobuccally skewed rhomboid crown outline which tapers distally and a waisted lingual margin. The cusp tips are displaced towards the buccal margin. First molars with high PC 1 scores have squared, symmetrical outlines and roughly equidistant cusp tips. Along PC 2, M 1 s with low scores have mesiobuccally skewed rhomboid outlines with cusp tips displaced towards the buccal margin, while M 1 s with high scores have more symmetrical crown outlines, increased buccolingual distance between the two mesial cusps and between the two distal cusps, and lingual displacement of the lingual cusps. There is substantial overlap among the three taxa, but hybrids tend to have low PC 1 scores and high PC 2 scores, while the parental taxa tend to have high PC 1 scores and low PC 2 scores. The 95% confidence ellipse for hybrids is much broader than those of the parental taxa, reflecting greater variation in shape space.
Principal component 3 accounts for 11.2% of M 1 shape variation (Figure 3b). First molars with low PC 3 scores have symmetrical and ovoid crown outlines and a rhomboid cusp tip configuration. High scores on PC 3 correspond to M 1 s with mesiobuccally skewed, rhomboid crown outlines with a waisted lingual margin, wide intercusp spacing, and all cusp tips displaced towards the periphery. There is very little separation among taxa in PC 2 vs. PC 3 shape space. The range of variation among S. nigritus individuals is almost entirely subsumed within the range of the hybrids. S. libidinosus tends to cluster on the low end of PC 2 away from S. nigritus and the hybrids. Based on shape and size of the 95% confidence ellipses projected onto tangent space for each taxon, the hybrids exhibit the highest variation in M 1 shape along PCs 2 and 3.

Regression analysis and allometry
The regression analysis demonstrated that taxonomic designation explains more variation in PC 1 scores than does M 1 size (Tables 2a and 2b). Approximately 20% (p < 0.001) of variation in PC 1 scores is explained by taxonomic designation. A post-hoc pairwise t-test using Bonferroni adjustment for multiple comparisons showed that hybrids had significantly lower scores on PC 1 than  . Scatter plots of (a) PC 1 against PC 2 scores, and (b) PC 2 scores against PC 3 scores derived from the PCA of M 1 shape. The warp grids illustrate the transformation of the consensus configuration (gray) into the shape of M 1 s with the lowest (blue) and highest (red) scores along PCs 1 and 2. Ellipses represent 95% confidence intervals for each group. Note that, while there is considerable overlap among the groups, the hybrids tend to exhibit lower PC 1 and higher PC 2 scores than the parental taxa, corresponding to M 1 s with skewed crown outlines, a waisted lingual margin, and wider intercusp distances. M: mesial; D: distal; B: buccal; L: lingual.

a)
The best-fitting model for each PC is in bold.

b)
The best-fitting model for each PC is in bold.  Figure 4a). More complex models testing the effect of taxonomic designation on the relationship between PC 1 scores and M 1 size were nonsignificant. Change in M 1 shape along the main axis of variation is not driven by size alone. However, the next-best-fitting model according to AIC suggested that the average PC 1 score estimated from M 1 size varies by taxon. A more complex model is required to explain variation in PC 2 scores. First molar size and taxonomic designation only explain 6% (p = 0.008) and 11% (p = 0.002) of variation in PC 2 scores, respectively, and a comparison of the two models indicates that PC 2 ~ taxon is a better fit than PC 2 ~ log(centroid size) (F = 5.93, p = 0.02). A post-hoc comparison of differences in PC 2 scores by taxon indicates that S. libidinosus has significantly lower PC 2 scores than the hybrids (p = 0.001; Figure 4b); all other pairwise comparisons are non-significant. A multivariate model combining the effect of taxonomic designation and M 1 centroid size explains 17% (p < 0.001) of variation in PC 2 scores and is a significantly better fit than PC 2 ~ taxon. There is no significant increase in explanatory power with the addition of an interaction term describing change in the slope of the relationship between M 1 shape and size among taxa (F = 1.12, p = 0.33). So, variation in shape along PC 2 is partly driven by size, but the average PC 2 score estimated from M 1 size differs by taxon.

Model
As with PC 1, taxonomic designation explains the most variation in PC 3 scores rather than M 1 size. However, the amount of variation in PC 3 scores explained by taxonomic designation is small (R 2 = 0.07, p = 0.02), and a post-hoc comparison average PC 3 scores by taxon indicates that the only significant difference among taxa is between S. libidinosus and the hybrids (p =0.016, Figure 4c). Comparisons with more complex models accounting for different size/shape relationships by taxon do not add significant explanatory power.

Discriminant function analysis
Results for the discriminant function analysis are illustrated in Figure 5. The DFA maximized differences in between-group variation, but there is little separation among parental species and their hybrids along linear discriminant functions (LDs) 1 and 2. Along LD 2, there is some separation between S. libidinosus, which clusters at the positive end, and S. nigritus and hybrids, which both cluster toward the negative end. Most of the hybrids have negative loadings on LD 1 and LD 2.
The results of the cross-validated assignment test are presented in Table 3. The percentage of individuals correctly classified to their a priori assigned taxon ranges from only slightly better than chance in S. nigritus (61.3%) to moderate in S. libidinosus (73.0%). In the hybrid group 70.5% were correctly assigned as such. Both S. nigritus and S. libidinosus misclassified individuals were more frequently assigned to the hybrid group than to the wrong parental taxon, reflecting higher variation in M 1 shape among hybrids.

Procrustes distances within and among taxa
The mean pairwise Procrustes distances in M 1 shape are listed in Table 4, and frequency distributions of pairwise Procrustes distances within and between taxa are visualized in Figure 6. The average distance for the entire sample is 0.06. Withintaxon shape variability is highest for the hybrids (distance = 0.059) compared to the parental taxa (S. nigritus = 0.055, S. libidinosus = 0.053). All three taxa have significantly different mean Procrustes distances (p < 0.001). The between-taxa comparisons show a greater degree of similarity between S. libidinosus and S. nigritus M 1 shape (distance = 0.058) than between each parental taxon and the hybrids, and there is approximately equal distance between the parental taxa and the hybrids (S.  . Based on Akaike's information criterion, the variation in PC 1 and 3 scores is best explained by taxonomic designation (Table  2b), while average PC 2 scores estimated from log(Centroid size) significantly differ among taxa (indicated by separate lines of best fit for each taxon; PC 2 ~ taxon + log(Centroid size)).   . Frequency distributions of (a) within-group, and (b) between-group pairwise Procrustes distances, reflecting degree of similarity in M 1 shape between specimen pairs. Vertical dashed lines represent the mean pairwise Procrustes distance for each group. Note that the hybrids exhibit elevated within-group Procrustes distances, reflecting the higher morphological variability in this group compared to the parental taxa. Also, between-group comparisons between one parental taxon and the hybrids exhibit higher mean Procrustes distances compared to the pairwise distances between parental taxa.

Covariation between cusp tip configuration and crown outline shape
The results of the 2B PLS analysis are summarized in Table 5. There is a strong and statistically significant correlation between cusp tip configuration (block 1) and crown outline shape (block 2) for the combined sample (r-PLS = 0.66, p = 0. 001). There are differences in covariation between blocks 1 and 2 among the three taxonomic groups. The two parental taxa exhibit high and significant correlations between blocks 1 and 2 (S. nigritus r-PLS = 0.76, p = 0. 002; S. libidinosus r-PLS = 0.75, p = 0.001), while the hybrids exhibit weaker correlation between cusp configuration and crown outline shape (r-PLS = 0.63, p = 0.015). Effect sizes for each sample indicate that the hybrids, compared to the parental taxa, exhibit weaker integration than expected based on its permutated sampling distribution. However, pairwise statistical comparisons indicate that there is no statistically significant difference in the strength of integration among groups (although at p = 0.07, the difference in integration between the hybrids and S. libidinosus does approach significance; Table 5b).

Discussion
The impact of hybridization on primate hard tissue morphology is difficult to predict. While traits under additive genetic control are expected to exhibit the midparental state in F1 hybrid populations, studies reveal that F1 morphology often deviates from expectations (Ackermann et al., 2006;Ito et al., 2015). Frequently, F1 hybrid trait morphology is polytypic and individuals exhibit novel phenotypes not observed in either parental population (Bergman et al., 2008;Fuzessy et al., 2014;Jolly et al., 1997). Many commonly measured phenotypic traits are not under additive genetic control. Hybridization may affect non-additive trait expression, resulting in heterosis or dysgenesis (Z. J. Chen, 2013). The recombination of two divergently adapted parental genomes in hybrids may disrupt the interaction and expression of non-additive genes that control complex physiological and metabolic networks, including growth and development. This ultimately relaxes the constraints observed in parental developmental pathways and is associated with increased morphological variability in hybrids. Deviations from expected midpa- Table 5. Results of the two-block partial least squares analysis.

a)
The r-PLS value reflects the degree of covariation between configuration of the cusp tips (block 1, landmarks 1 through 4) and the shape of the crown outline (block 2, sliding semilandmarks 5 through 34). Larger effect sizes are associated with stronger observed covariation between cusp tip and crown outline shape than expected based on the permutated sampling distribution.  (Allen et al., 2020;Bernardes et al., 2017;Stelkens & Seehausen, 2009). For example, there are fewer instances of cranial and postcranial trait heterosis in hybrids of recently diverged tamarin subspecies than between crosses of more anciently diverged tamarin subspecies (Cheverud et al., 1993;Kohn et al., 2001). Similarly, non-metric indicators of disrupted skeletodental development tend to be more frequently observed in primates with increasing parental divergence Boel et al., 2019). Beyond the first generation, hybrid morphology is expected to more closely resemble that of the parental population into which the hybrids have backcrossed (Falconer & Mackay, 1997). Continuous trait values in the backcrossed offspring of an F1 hybrid and an individual from the parental population are predicted to be the average of the parental value and the MPV. Some novel phenotypes observed in F1 hybrids persist in latergeneration hybrids regardless of parental genetic contribution. Macaca fuscata x M. cyclopis macaques have enlarged, M. fuscata-like sinus size even in backcrossed individuals who derive most of their ancestry from M. cyclopis (Ito et al., 2015), and transgressive non-metric dental traits are observed in backcrossed P. cynocephalus x P. anubis individuals . The morphology of individuals in multigenerational hybrid zones depend on a combination of physiological, reproductive, and ecological selective pressures (Charpentier et al., 2008;Fourie et al., 2015;Jolly et al., 2011;Mourthe et al., 2019). These selection pressures structure the distribution of hybrid phenotypes across contact zones. For example, hybrids from the contact zone between P. anubis and P. cynocephalus in Amboseli, Kenya exhibit a continuous distribution of phenotypes ranging from more P. anubis-like to intermediate to more P. cynocephalus-like, while the phenotypic distribution of hybrids in the P. anubis x P. hamadryas contact zone in Awash, Ethiopia is bimodal with very few intermediate phenotypes (Alberts & Altmann, 2001;Wango et al., 2019). Phenotypically intermediate hybrids in Awash also exhibit reproductive behaviors intermediate to those observed in parental taxa, and are therefore thought to be at a reproductive disadvantage when backcrossing with P. anubis or P. hamadryas compared to hybrids with predominantly parental phenotypes and behaviors (Bergman et al., 2008). Hybrids in recently formed anthropogenic contact zones show more continu-ous phenotypic distributions and symmetrical contribution of parental genes into the contact zone (Malukiewicz, 2019). So, a biologically relevant understanding of phenotypic outcomes in hybrid populations requires information regarding a variety of endogenous and exogenous variables.

Matrix of pairwise differences in 2B PLS effect size measuring difference in the strength of integration between samples in the lower
This study assumes that the NMNH is correct in its taxonomic designations of the specimens used. However, the hybrids and parental taxa studied here have not been genotyped, as is preferable in analyses examining the relationship between phenotype and degree of hybridity (Ackermann et al., 2006;Boel et al., 2019;Cheverud et al., 1993;Hamada et al., 2012). The assumption that the parental taxa are not themselves admixed may be particularly problematic for robust capuchins, as Sapajus species have a complex history of hybridization in secondary contact zones (Lima et al., 2018). In addition, there may be cryptic hybrids in the sample with a high degree of genetic admixture but no phenotypic indication of hybridity (Ackermann, 2010;Kelaita & Cortés-Ortiz, 2013). Regardless, the results of the analyses presented here, combined with results from previous research, allow for some predictions to be made regarding the genetic makeup of the robust capuchin hybrids.
My analyses indicate that, while hybrid M 1 shape largely falls within the range of variation observed in S. nigritus and S. libidinosus, some aspects of hybrid M 1 shape are unique compared to parental morphology. While PC 1 typically captures the allometric component of shape variation (Zelditch et al., 2012), in this study taxonomic designation explained a greater proportion of variation in PC 1 scores than did molar size (Table 2). Hybrids significantly differed from S. nigritus and S. libidinosus on the main axis of M 1 shape variation in the PCA. Hybrids had significantly lower PC 1 scores than both parental taxa and higher PC 2 scores than S. libidinosus, corresponding to M 1 s with increased buccolingual distance between cusps and a waisted lingual crown margin (see Figures 3 and 4). However, hybrids and S. nigritus did not exhibit significantly different PC 2 or PC 3 scores. So, some hybrids exhibit a unique molar morphotype compared to parents, while others cluster with S. nigritus. This was reflected by the reasonably accurate classifications generated by the DFA assignment test (Table 3). Sapajus nigritus had the lowest correct assignment (61.3%), with more individuals misclassified as hybrids than S. libidinosus (Table 3). Sapajus libidinosus specimens also were more often misclassified as hybrids than S. nigritus but exhibited the highest percentage of correctly classified individuals (73.0%).
The results of the PCA and DFA support recent revisions in capuchin taxonomy. Based on genetic and morphological data, the capuchins are proposed to contain two genera: the gracile Cebus capuchins and the robust Sapajus capuchins (Lynch Alfaro, de Sousa e Silva-Júnior, et al., 2012). The IUCN recognizes eight Sapajus species that can be subdivided into a more ancient clade that evolved in the Brazilian Atlantic forest and a clade that recently left the Atlantic Forest to spread throughout the Amazon (Lima et al., 2018). Sapajus nigritus belongs to the more ancient clade, and retains morphological features indicative of arboreal living, such as longer limbs and tails. Sapajus libidinosus belongs to the Amazonian clade but has recently evolved morphological traits for terrestrial life in the dry shrublands of the Brazilian Cerrado-Caatinga, including thickened molar enamel and shorter, more robust limbs (Wright et al., 2015). Sapajus libidinosus is therefore the most morphologically derived robust capuchin species (Wright et al., 2015). The results of the PCA and DFA presented here indicate that S. nigritus and S. libidinosus exhibit statistically significant differences in M 1 shape. Hybrids cluster more with S. nigritus rather than with the more derived S. libidinosus. Based on the tendency of the hybrids and S. nigritus to cluster along PCs 2 and 3, I would expect the hybrids to exhibit greater genetic affinity with S. nigritus. Tail length has been shown to track degree of hybridity in macaques (Hamada et al., 2012), so it would be interesting to test this in S. nigritus x S. libidinosus hybrids.
Mean hybrid M 1 shape in this study is not the MPV (see Figure 2). Compared to the expected shape, the observed hybrid M 1 mean shape exhibits buccolingual expansion and a waisted lingual margin. However, the MPV is expected only in F1 hybrids and only for traits under additive genetic control (Falconer & Mackay, 1997). It is highly unlikely that wild hybrid populations contain only F1 individuals (Kelaita & Cortés-Ortiz, 2013;Phillips-Conroy & Jolly, 1986). Additionally, it is known that the genetic architecture controlling M 1 size and shape is partly non-additive (Hardin, 2019;Hlusko et al., 2016). Combined, these observations indicate that the deviation from the expected midparental M 1 shape observed in this study are likely caused by the disruption of non-additive gene expression or epigenetic interactions in latergeneration S. nigritus x S. libidinosus hybrids. This suggests that the morphological impact of hybridization persists beyond early hybrid generations, as has been demonstrated in baboons and macaques Ito et al., 2015).
The S. nigritus x S. libidinosus hybrids exhibit evidence of destabilized dental development. Measured by pairwise Procrustes distances, hybrids exhibit statistically significant elevation of within-taxon variation in M 1 shape compared to both parental taxa. This variation may be driven by relaxed constraints during dental development (Fuzessy et al., 2014). Indeed, hybrids exhibit lower mean correlation between cusp tip configuration and crown outline shape (r-PLS = 0.63) compared to S. nigritus and S. libidinosus (r-PLS = 0.76 and r-PLS = 0.75, respectively). Hybrids tend to have wider intercusp distances and cusps positioned closer to the crown periphery than the parental taxa. Cusp tips correspond to the position of embryonic signaling centers in developing tooth germs. The distance between cusp tips is controlled by the relative strengths of activator and inhibitor molecules excreted by each signaling center and the duration of germ growth. Increased inhibitory signaling and/or prolonged germ growth are expected to result in fully formed teeth with widely spaced cusp tips (Guatelli-Steinberg et al., 2013;Jernvall, 2000). So, the wide intercusp distances and weaker correlation of cusp configuration and crown outline shape observed in S. nigritus x S. libidinosus hybrids are likely the result of prolonged dental development and/or deviation in levels of signaling molecules compared to those observed in parental dental development. Similarly, Ackermann et al. (2014) found that the presence of supernumerary distomolars is associated with increased molar row length in F1 hybrid P. cynocephalus x P. anubis individuals, suggesting that dental development is prolonged in hybrids compared to parents. Among other papionin hybrids, Papio hamadryas x P. anubis hybrids exhibit unique molar size relationships compared to parental taxa, suggesting that developmental pathways controlling hybrid baboon molar size may be destabilized compared to unadmixed baboons (Phillips-Conroy, 1978). However, based on frequencies of dental non-metric trait expression and fluctuating asymmetry of bilateral cranial traits, there is no evidence for destabilized dental development in earlygeneration M. fuscata x M. cyclopis macaques (Boel et al., 2019). It is possible that these observations support the prediction that the degree of developmental destabilization observed in hybrids is associated with parental divergence. Sapajus nigritus and S. libidinosus shared a common ancestor around 2.6 Ma (Lima et al., 2018); P. cynocephalus and P. anubis diverged approximately 1.5 Ma while P. hamadryas and P. anubis diverged approximately 800 ka (Rogers et al., 2019); and M. fuscata and M. cyclopis are estimated to diverge as recently as 170 ka (Chu et al., 2007). A comparison of dental phenotypic variation and integration among these different hybrid populations would confirm the relationship between the degree of parental divergence and destabilized development in hybrids.
While non-metric dental anomalies are observed at high frequencies in some mammalian hybrid populations, this pattern is not shared by all extant primates. This calls into question the suggestion that certain dental non-metric traits, especially supernumerary distomolars or dental crowding, are evidence of significant hybrid ancestry in extinct hominins (Ackermann, 2010;Ackermann et al., 2019). However, continuous dental trait variation remains understudied, even though non-metric dental traits are often correlated with continuous trait variation (Ortiz et al., 2018) and a Homo sapiens fossil with substantial H. neanderthalensis ancestry exhibits extremely large upper third molars (Fu et al., 2015). The results presented here suggest that transgressive M 1 morphology that falls outside of the range of variation observed in well-defined hominin taxa may be indicative of hybrid ancestry in hominin fossils. Further analyses comparing molar shape variation in other extant primate hybrids would confirm if this is a valid prediction. In terms of primate conservation, this study did not indicate that hybridization reduced phenotypic variation among hybrids of S. nigritus and S. libidinosus. Rather, hybridization generated novel phenotypes not observed in either parental population. It remains to be determined if expanded intercusp distances in these hybrids facilitate ecological niche separation from other robust capuchin populations.

Conclusions
The dentition has been an anatomical region of interest in hybrid research, but previous work has predominantly studied non-metric dental trait variation rather than tooth shape. The results presented here suggest that a more in-depth analysis of the impact of hybridization on continuous dental phenotypes and development is warranted. The shape of the first upper molar is statistically distinct among S. nigritus, S. libidinosus and their hybrids, and hybrids exhibit morphological evidence of destabilized development, including elevated within-sample variance and weaker correlation between cusp tip configuration and crown outline shape. The same analyses used here applied to the rest of the postcanine teeth would likely uncover other significant differences between the hybrids and parental taxa. A more comprehensive understanding of the impact of hybridization on dental development could be gained by further comparisons of continuous trait integration between metameres and between occluding upper and lower molars; and by comparing levels of fluctuating asymmetry of continuous traits in left and right antimeres. The data derived from such studies would offer crucial information for attempts to diagnose hybrid ancestry from fossil morphology and to understand the evolutionary outcomes of hybridization among endangered primates in degraded habitats.