Cranial Anatomical Integration and Disparity Among Bones Discriminate Between Primates and Non-primate Mammals

The primate skull hosts a unique combination of anatomical features among mammals, such as a short face, wide orbits, and big braincase. Together with a trend to fuse bones in late development, these features define the anatomical organization of the skull of primates—which bones articulate to each other and the pattern this creates. Here, I quantified the anatomical organization of the skull of 17 primates and 15 non-primate mammals using anatomical network analysis to assess how the skulls of primates have diverged from those of other mammals, and whether their anatomical differences coevolved with brain size. Results show that primates have a greater anatomical integration of their skulls and a greater disparity among bones than other non-primate mammals. Brain size seems to contribute in part to this difference, but its true effect could not be conclusively proven. This supports the hypothesis that primates have a distinct anatomical organization of the skull, but whether this is related to their larger brains remains an open question.


Introduction
The skull of primates shows a peculiar combination of anatomical features. Compared to other mammals, primates in general have shortened snouts and wide forward-facing eye orbits (Fleagle et al., 2010). Primates also tend toward a more vertical posture regardless of locomotion behavior, which affected their cranial base morphology (Fleagle, 1999). One of the most remarkable feature of primates' evolution is the enlargement of the brain relative to body size, especially in some groups like anthropoids (Boddy et al., 2012). Evolving bigger brains imposes functional and anatomical trade-offs among the different parts of the skull. For example, the fact that the back of the face is physically integrated with the cranial base limits changes in the orientation of the cranial base and in encephalization (Lieberman, 2011;Lieberman et al., 2000;McCarthy & Lieberman, 2001).
The loss of ossification centers and the fusion of bones during development has also influenced the anatomy of the primate skull (Esteve-Altava et al., 2013;Goodrich, 1931;Kardong, 2005;Koyabu et al., 2012Koyabu et al., , 2014Rager et al., 2014;Richtsmeier et al., 2006). Modern mammals have lost bones that were once present in their synapsid ancestors, such as the prefrontal, postfrontal, or postorbital Sidor, 2001). In turn, primates fuse bones that tend to remain unfused in other mammals, resulting in bones such as the occipital, temporal, and sphenoid ; haplorhines also fuse the left and right frontals into a single frontal bone (Rosenberger & Pagano, 2008), and in humans the premaxilla fuses to the maxilla (Lang, 1995). Comparative studies in tetrapods (Esteve-Altava et al., 2013), synapsids (Sidor, 2001), and archosaurs (Lee et al., 2020;Plateau & Foth, 2020;Sookias et al., 2020;Werneburg et al., 2019) have shown that the loss and fusion of bones can modify the anatomical organization of the skull-which bones articulate to each other and the pattern of interactions this creates-in a way that increases the anatomical integration of the skull and the disparity among bones through evolution. In this context, integration and disparity refer to the aforementioned anatomical organization: integration refers to how strongly skull bones are interconnected (as in number of connections among elements and of 3-node loops), and disparity refers to how different or heterogeneous 1 3 bones are in terms of number of connections. These general properties can be proxy-measured in different ways (see Box 1). In turn, we can understand modularity as a pattern of structured integration, with groups or clusters of bones more integrated (connected) among them than to bones in other clusters, which allows different regions of the skull to still evolve semi-independently (Esteve-Altava, 2017a). Indeed, studies using network models found that anatomical integration is strong in the skull of primates (Esteve-Altava et al., 2015;Powell et al., 2018). However, these studies did not compare primates with non-primate mammals. Curiously, changes in anatomical organization that increase integration and disparity among bones seem to be more common in lineages that evolved bigger brains and globular neurocrania (e.g., birds compared to archosaurs (Marugán-Lobón et al., 2016;Plateau & Foth, 2020), which is probably a consequence of the strong ties between brain and skull development (Richtsmeier & Flaherty, 2013). This relationship is even stronger for the neurocranial bones which intramembranous growth is linked to brain growth (Opperman, 2000). Building on that, we can propose the working hypothesis that primates have an anatomical organization distinct from that of other mammals, which is related to bigger brains.
Here, I quantified the anatomical organization of the skull of primates and non-primate mammals to test: (i) whether primates show a distinctive organization of the skull that sets them apart from other mammals; if so, (ii) what specific anatomical features create this difference; and (iii) how such anatomical features relate to the increase in brain size relative to body size characteristic of primates. To this end, the skulls of 17 primates and 15 non-primate mammals were modelled as networks and their anatomical organization measured using six topological variables. See Fig. 1 for a visual schema of the topological variables and the list of species compared. Topological variables were compared in a phylogenetic context to identify differences between primates and other nonprimate mammals and to assess the relationship between cranial anatomical organization and brain size.

Box 1. Morphological Interpretation of Topological Variables
Density of connections (D) An aspect of anatomical complexity and integration related to biological burden (sensu Riedl (Schoch, 2010;Rasskin-Gutman & Esteve-Altava, 2018)), so that the more interactions among bones a skull has the more complex and integrated its functioning and the more constrained to change.
Mean clustering coefficient (C) An aspect of anatomical integration related to the interdependency among bones. The more interdependency (A links to B, B links to C, and C links to A), the more integrated are the bones of the skull.
Mean shortest path length (L) An aspect of anatomical integration related to functional effectiveness. The closer two parts are in a topological sense (fewer links separating them), the more probably they will interact and function together.
Parcellation (P) An aspect of anatomical organization related to the level of modularity, how evenly distributed are the bones into as many modules as possible. The more parcellated are the bones into many modules, the more modular is the skull.
Heterogeneity of connections (H) An aspect of disparity or anisomerism among bones related to the irregularity of the skull. The more diverse the number of interactions among bones, the more irregular is the anatomical organization (e.g., a honeycomb has a nearly regular organization, almost all cells have six neighboring cells).
Assortativity (A): an aspect of disparity or anisomerism among bones related to their tendency to articulate to other bones with similar number of articulations. The more similar the number of connections among connected bones, the more regular is the anatomical organization.

Anatomical Network Models
Anatomical network models of the typical young adult skull of 17 primates and nine non-primate mammals were gathered from existing literature (Esteve-Altava et al., 2013. The original studies defined age by dentition and 'typical' by the frequency of presenting or not an articulation, which tends to be a relatively stable feature. Six additional non-primates were coded to broaden the comparison: Antilope cervicapra, Dasypus novemcinctus, Equus caballus, Felis silvestris, Oryctolagus cuniculus, and Sus scrofa (see details in Supplementary Materials). Figure 1B shows the complete list of species. Following the modeling approach of previous studies (Rasskin-Gutman & Esteve-Altava, 2014): every bone was coded as a network node and every suture or synchondrosis between two bones was coded as a link connecting two nodes. Networks available at https:// doi. org/ 10. 6084/ m9. figsh are. 14662 257. v1.

Topological Variables
Six topological variables were measured: density of connections (D), mean clustering coefficient (C), mean shortest path length (L), heterogeneity (H), assortativity (A), and parcellation (P). D is the actual number of connections divided by the maximum theoretically possible given the number of nodes; D ranges from 0 to 1. C is the actual number of 3-node loops divided by the maximum theoretically possible given the number of nodes; C ranges from 0 to 1. L is the averaged minimum number of links separating all pairs of nodes divided by the number of pairs; L equals 1 when all nodes are directly connected and gets greater with longer distances between pairs of nodes. H is the coefficient of variation (σ/μ) of the number of connections per node; H ranges from 0 (all nodes have the same number of connections) to √ N − 1 (maximum disparity among nodes' connections). A is a preference for the nodes of a network to attach to other nodes that have a similar number of neighbors.; A ranges from -1 (dissimilar connectivity) to 1 (similar connectivity). Finally, P measures the extent to which nodes are evenly distributed into modules; P ranges from 0 (no modules) to 1 (as the number of modules increases). Networks have a higher P when they have more modules and also when nodes are similarly spread among them (e.g., a network with four modules of five nodes each has a higher parcellation than a network with one module of 17 nodes and three modules of two nodes). For a mathematical description of topological variables see . Topological variables were measured using functions from the R package igraph (Csardi & Nepusz, 2006) and each variable rescaled for the whole sample using min-max normalization.

Calibrated Phylogeny and Brain Size Residuals
A consensus calibrated phylogeny for the 32 species was built from 100 trees pruned from the VertLife project (Upham et al., 2019) (Supplementary Data), using the mean edges method as implemented in the function consensus. edges of the package phytools (Revell, 2012). Brain size residuals were collected from a recent large-scale compilation by Burger et al. (2019), which includes all species in the study except Cynocephalus volans. This study reports an allometric relation of Brain = −1.26 × Body 0.75 , which I used to estimate the residual of C. volans using brain and body mass measurements from Lewitus et al. (2014). Raw data is available at https:// doi. org/ 10. 6084/ m9. figsh are. Their actual meaning depends on the nature of the relations that are being modelled as links (Butts, 2009). In skull networks links represent sutures and synchondroses, which embody growth and functional relations among bones (Herring & Teng, 2000;Opperman, 2000;Rafferty et al., 2003;Di Ieva et al., 2013), therefore, topological variables translate into morphological features that depend on growth and function, such as integration, modularity, and disparity (see Box 1). The six topological variables used here are illustrated with colors. D is the number of links (in gray) given the maximum number of links possible, which depends on the total number of nodes. C is the number of 3-node loops given the maximum number of loops possible (in blue, the 3-node loop between temporal, parietal, and occipital). L is the mean of all shortest paths (in green, shortest path of length two between the left maxilla and temporal). P is the even partition of nodes into many modules (dashed red lines indicate the three modules identified using an optimization algorithm, see Methods). H is the heterogeneity of connections among nodes and A is the correlation of connectivities: note how bones vary in their connections (the black circle size within nodes is proportional to their number of connections relative to the most connected node, the ethmoid with 13); in yellow, the connections of the four links of the lacrimal and of the 12 links of the sphenoid illustrate the two extremes of the variation. B Evolutionary relations of sampled mammals. The skulls of primates (including representatives of the Strepsirrhini and Haplorhini suborders) was compared with a diverse set of mammals spanning 12 orders, including marsupials and monotremes (Color figure online) 1 3 14662 257. v1 (phylogeny: output2.nex; brain data: Burger-etal2019brainallometry.tsv).

Statistical Analysis
All statistical analyses and visualizations were performed in R (R Core Team, 2019) (available within Supplementary Results). A principal components analysis (PCA) was performed (without rotation) on the six topological variables with the function principal from the package psych (Revelle, 2019). The relationship between number of components, % of explained variance, and the load of variables on components were examined to decide the final number of components, as a compromise between capturing as much % of explained variance as possible but with each variable loading mostly onto only one component (Online Resource ,  Fig. S3). As a result, two obliquely rotated components were extracted (RC1 and RC2). Results were visualized using a phylomorphospace projection using the homonym function of the package phytools (Revell, 2012). Differences between primates and non-primates were tested with phylogenetic permutational multivariate analyses of variance (PER-MANOVA) on raw topological variables and RCs, using the function aov.phylo of the package geiger (Harmon et al., 2008) with 10,000 permutations.
The association between RCs and brain size residuals was tested for the whole sample and for each group separately using phylogenetic methods under the Pagel's λ model. Phylogenetic Pearson's product moment correlations were calculated with a custom script, using the function phyl.vcv in phytools to measure the variance-covariance matrix of variables. Phylogenetic generalized least squares (PGLS) regressions were performed using the function gls from the package nlme (Pinheiro et al., 2019). To test whether a single regression for the whole sample fits the data better than separate regression models for primates and non-primates, a phylogenetic analysis of covariance (PANCOVA) was performed using the function gls.ancova of the package evomap (Smaers & Mongle, 2020).
The assumption of normality of the residual error for PGLS models was tested with a Lilliefors (Kolmogorov-Smirnov) test. Outliers were first visualized with boxplots (interquartile rule) and tested using PANCOVA, which has been shown to be the only correct way to directly test whether a single species is an outlier (Smaers & Rohlf, 2016).

Anatomical Network Analysis
The 32 skull network models analyzed range in number of bones between 21 and 32 and in number of articulations between 47 and 99 (Online Resource, Table S1), producing a diverse set of networks, as captured by topological variables (Online Resource, Table S2). The phylogenetic PERMANOVA test supported a statistical difference in topological variables between primates and non-primates (Wilks' lambda = 0.1596, F 6,25 = 21.9322, p-value = 0.0097). A PCA of the topological variables was performed to visualize the overall pattern of variation across skulls and to select an adequate number of components. Two PCs were enough to explain 78.71% of the variance without increasing variable loads on a single PC above 1.5 (Online Resource, Fig.  S3), which makes possible to interpret components anatomically because one variable is captured largely by only one component.
A PCA with oblique rotation was performed to extract the two main rotated components (RC). RC1 and RC2 explained 49.96% and 28.74% of the variance, respectively (Online Resource, Table S4). RC1 captured topological variables D, C, L, and P (Online Resource, Table S5), which are related to anatomical integration (Box 1). Greater values of RC1 align with greater values of D and C, and lower values of L and P; thus, RC1 ranks skull networks from less to more anatomically integrated. In turn, RC2 captured the variables H and A (Online Resource, Table S5), which are related to anatomical disparity (Box 1). Greater values of RC2 align with greater values of H and lower values of A; thus, RC2 ranks skull networks from less to more disparate. When analyzed in a phylogenetic context (Fig. 2), RCs also discriminate significatively between primates and non-primates (Wilks' lambda = 0.3865, F 2,29 = 23.012, p-value = 0.0298). These results support the working hypothesis that primates have an anatomical organization distinct from that of other mammals.

Network Anatomy and Brain Size
There was a significant positive phylogenetic correlation between RC1 and brain size residuals for the species in the sample (r = 0.5467, t = 3.5766, p-value = 0.0012). The PGLS of brain size residuals on RC1 confirmed bigger brains predict greater values of RC1 (slope = 1.3769 ± 0.6327, t = 2.1763, p-value = 0.0375; Online Resource, Fig. S6) in mammals. This baseline model fitting a single slope and intercept for the whole sample was compared with alternative models where primates and non-primates fit to different intercepts and/or slopes (Fig. 3). PANCOVA tests showed that none of the alternative models fitted the data significantly better than the baseline model (Online Resource, Tables S10, S11, S12). A closer examination of the relation between variables for primates and non-primates (Online Resource, Tables S13, S14, S15, S16) showed that only primates had a significant correlation between RC1 and brain size residuals (r = 0.6356, t = 3.1888, p-value = 0.0061), with a linear relation similar to that of the baseline model (slope = 2.2078 ± 0.8786, t = 2.5129, p-value = 0.0239). This may indicate that the correlation seen in the whole sample is mainly driven by primates alone. These results support an alternative working hypothesis, in which brain size affects anatomical integration overall, but brain size cannot explain differences of integration between primates and non-primates.
PGLS residuals of brain size on RC2 failed the Lilliefors test of normality (D = 0.1788, p-value = 0.0107) and showed two outliers (Pan and Cercopithecus), visually identified using the interquartile rule and confirmed with a PAN-COVA for differences in intercept between the outliers and the whole sample (F 2,3 = 29.9505, p-value < 0.0001). After removing the confirmed outliers the PGLS residuals passed the normality test and showed no additional outliers (Online Resource, Table S19). However, no significant relationship was found between brain size residuals and RC2 for the whole sample (Online Resource, Tables S19 and S20) nor for primates and non-primates separately (Online Resource, Tables S21, S22, S23, S24). Thus, rejecting the working hypothesis that anatomical disparity of the skull is driven by brain size.

Discussion
These results support three points: (i) topological variables measured on skull network models can discriminate between primates and other non-primates mammals based on their anatomical organization; (ii) two features distinguish primates: a greater anatomical integration of the skull and a greater disparity among bones; and (iii) brain size seems to influence anatomical organization differently in primates and non-primates, but its true effect could not be conclusively demonstrated. Results support the working hypothesis that primates have a distinct anatomical organization of the skull compared to non-primate mammals, but whether this is related to primates evolving bigger brains remains an open question.
Primates have skulls with greater anatomical integration and disparity than the other mammals studied. This finding builds on the morphological interpretation of topological variables (Box 1) and their respective loading on each RC (Online Resource, Table S5). As introduced before, topological variables capture anatomical integration and disparity because connections in skull network models embody growth and functional relations between bones (Herring & Teng, 2000;Opperman, 2000;Rafferty et al., 2003;Di Ieva et al., 2013). Having a greater cranial anatomical integration means that, collectively, the bones of the skull interact more with one another, which affects form and function. Because a connection is a shared growth site between two bones (Opperman, 2000), changes in the shape of one bone, caused by genetic or epigenetic factors modifying growth, are more likely to produce a collateral change in an adjacent Fig. 2 Phylomorphospace of skull network models defined by the two rotated components. Primates and non-primates occupy different regions of the space, with primates having larger values of RC1 and RC2 meaning that they have more anatomically integrated skulls and a greater disparity among bones than the other non-primate mammals sampled. Interestingly, the bat Pteropus directly occupies a position within the space occupied by primates (dashed line circle). The fact that bats share with primates a pervasive pattern of bone fusions during development may explain having a similar degree of anatomical disparity and disparity among bones Fig. 3 PGLS fits of brain size on anatomical integration. The baseline model (dashed black line) fitted the data better than the alternative models where primates (in green) and non-primates (in blue) have different slope and/or intercept. Note that the fit represented by the blue like is not significant and it has been included only for comparison and to show the alternative models compared (Color figure online) bone (Pearson & Woo, 1935;Woo, 1931), and indirectly to bones adjacent to that one (and so on). Therefore, the more structurally integrated (i.e., interconnected) is the whole skull, the stronger its shape integration and the weaker its modularity. This type of relationship between anatomy and shape would be a baseline expectation (actually, one that remains untested), with other factors also weighing in to determine the final phenotypic variation of the skull.
As related as they may seem, anatomical or structural integration and morphological integration are not interchangeable concepts, because the actual relation between organizational units (i.e., network anatomy) and variational units (i.e., shape co-variation) is not yet fully understood (Eble, 2005;Esteve-Altava, 2017b;Rasskin-Gutman & Buscalioni, 2001). Patterns of connectivity emerge from, and act on, growth and function: two of the main drivers of shape variation and morphological adaptation. The same is true for morphological integration (as shape co-variation), which is also characteristically strong in the skull of primates (Neaux et al., 2018(Neaux et al., , 2019. From a more functional perspective, bone connections act also as sites of energy absorption and deformation (Herring, 2008), working together to diffuse strains across the skull (Moazen et al., 2009). Thus, we could speculate that the more interconnected is the skull, the more efficient its biomechanical functioning in terms of redistributing stress loads. Additionally, the manner in which topological variables capture anatomical integration (Box 1) is related to a kind of structural compactness or robustness (see Sidor, 2001 for a discussion of this argument in relation to bone losses in the skull of synapsids). Thus, as primates evolved bigger brains the need for more robust and protective skulls might have co-evolved.
On the other hand, having a greater disparity (in number of connections) among bones means that there are some bones that bear many interactions while others bear few. Results indicate this is what happens in primates, as it happens in modern mammals compared to stem mammals (Esteve-Altava et al., 2013).. In mammals, the number of connections seems to be related to the midline location of bones, which allows to articulate with bones from both sides of the skull; however, there are enough exceptions that prevent us to establish a constant relationship between location and connectivity. For example, in the human skull the ethmoid, frontal, and sphenoid have 13, 12, and 12 links, while most bones have 4-6 links (Fig. 1); the other two unpaired bones, vomer and occipital, have an average number of connections. Network studies have shown that systems with a few highly connected nodes and many poorly connected nodes are more robust to accidental damage (e.g., injury or loss), because accidents are more likely to happen to the most abundant type of node, which are the less connected ones. In contrast, the same systems are sensible to damage targeting the most connected nodes (Albert et al., 2000;Allesina et al., 2009), because they are more essential. This phenomenon has also been demonstrated in anatomical networks (Esteve-Altava et al., 2013Murphy et al., 2018). Because more connected bones bear more growth and functional interactions, they are more constrained to change (Rasskin-Gutman & Esteve-Altava, 2018;Schoch, 2010). In contrast, bones with fewer connections can vary, specialize, or get lost more freely without compromising the integrity of the skull. However, when considering the potential effects of connections in constraining bones variation, it is important to also consider their relations with the surrounding soft tissues and functional matrices, in particular, the central neural and circulatory systems (Esteve-Altava & Rasskin-Gutman, 2014).
In mammals, a greater anatomical integration of the skull is related to bigger brains. The mechanisms by which bigger brains may produce more integrated skulls is unknown. A plausible explanation points to the role of the brain as a functional matrix in development (Moss & Young, 1960), and the way in which the bones and articulations of the neurocranium respond to brain growth (Lesciotto & Richtsmeier, 2019;Richtsmeier & Flaherty, 2013;Richtsmeier et al., 2006). It is also possible that bigger brains impose a selective pressure for more robust and protective skulls, which can be achieved by closing sutures and fusing bones together (Sidor, 2001), incidentally increasing anatomical integration and sometimes anisomerism . Given that many lineages of primates also show trends toward bigger orbits and shorter faces , bigger brains may also produce a tighter packing of midfacial and basicranial bones (Bastir et al., 2010;Lesciotto & Richtsmeier, 2019). At the level of skull sutures and synchondroses, compression creates an environment that favors suture closure, for example, it enhances osteogenesis, narrows suture space, and immobilizes bones (Herring, 2008), which, as said before, is a process known to increase anatomical integration . The correlation between anatomical integration and bigger brains is mostly driven by primates; non-primate mammals show no relationship between brain size and cranial integration. This is mostly because primates have a statistically significant linear fit between brain size and integration, while nonprimates do not. To sustain this preliminary result, future studies will need to look at this relationship with similarly large samples in other orders of mammals.
In summary, network models and topological variables paint a picture of the skull that is not accessible to, and complements that of, other methods in morphology (Eble, 2005; Rasskin-Gutman & Esteve-Altava, 2014). Primates evolved a particular anatomical organization-in liaison with bigger brains?-that distinguish them from non-primate mammals. Future studies should address whether other lineages with relatively big brains, such as dolphins within cetaceans (Marino, 2009) or corvids within birds (Uomini et al., 2020), convergently evolved a similar relationship between cranial anatomy and brain size. On a more general note, the results presented here add up to a growing body of research that seeks to better understand morphological evolution, development, and function by using methods that allow us to studying the body's anatomy as a complex biological system (Esteve-Altava et al., 2013Saucède et al., 2015;Dos Santos et al., 2017;Kerkman et al., 2018;Murphy et al., 2018;Werneburg et al., 2019;Fernández et al., 2020;Plateau & Foth, 2020;Sookias et al., 2020).