Anatomical Network Analysis Shows Decoupling of Modular Lability and Complexity in the Evolution of the Primate Skull

Modularity and complexity go hand in hand in the evolution of the skull of primates. Because analyses of these two parameters often use different approaches, we do not know yet how modularity evolves within, or as a consequence of, an also-evolving complex organization. Here we use a novel network theory-based approach (Anatomical Network Analysis) to assess how the organization of skull bones constrains the co-evolution of modularity and complexity among primates. We used the pattern of bone contacts modeled as networks to identify connectivity modules and quantify morphological complexity. We analyzed whether modularity and complexity evolved coordinately in the skull of primates. Specifically, we tested Herbert Simon’s general theory of near-decomposability, which states that modularity promotes the evolution of complexity. We found that the skulls of extant primates divide into one conserved cranial module and up to three labile facial modules, whose composition varies among primates. Despite changes in modularity, statistical analyses reject a positive feedback between modularity and complexity. Our results suggest a decoupling of complexity and modularity that translates to varying levels of constraint on the morphological evolvability of the primate skull. This study has methodological and conceptual implications for grasping the constraints that underlie the developmental and functional integration of the skull of humans and other primates.


Introduction
Modularity and complexity are two biological phenomena that occur together in the development and evolution of the skull in humans as well as in non-human primates [1]; however, rather than simplifying it as proposed by previous studies of other authors (see above). This result is a consequence of defining complexity by richness of interaction among bones (e.g. the number of bone contacts and the formation of triangular-loop motifs), instead of merely by the number of bones. Moreover, Esteve-Altava and co-workers [30] proposed a specific mechanism underlying this increase in complexity: the random loss of bones with fewer contacts to other bones, and the preferred fusion of bones with more contacts to each other. Interestingly, this reduction in bone number was also concomitant with an increase in specialization and disparity of skull bones [30], which Gregory [31] called anisomerism. AnNA quantifies anisomerism as heterogeneity in the number of bone's connections [29]. Because anisomerism provides the skull with distinct parts, different from each other, it can also be considered as a measure of complexity.
Other previous studies have used AnNA to identify the modular organization of the human skull [27], and to characterize the modularity-complexity interplay that takes place during ontogeny to produce the adult human skull [28]. These studies identified two connectivity modules in the skull of modern humans that behave as units of allometric growth [27]: a facial module with a hierarchical arrangement of bones in smaller sub-modules (or blocks) around the ethmoid bone, and a cranial module with a pseudo-regular arrangement of bones around the sphenoid bone. Moreover, together the ethmoid and the sphenoid bones participate in almost 40% of the total bone-bone contacts of the human skull, and are thus arguably of developmental and functional significance. It is not yet clear whether the same patterns hold true in the skulls of non-human primates; here, we will analyze these patterns to get a better understanding of the evolutionary processes that ultimately led to the derived form of the human skull.
Our present analysis is therefore, to our knowledge, the first unified study of the phenotypic modularity and complexity of the skull for representatives of all major primate taxa, as well as of the two living groups more closely related to Primates: Dermoptera (or colugos represented by Cynocephalus) and Scandentia (or tree shrews, represented by Tupaia) (Fig 2). We used AnNA to characterize the modular organization of the skull in these taxa, and to test specifically whether modularity correlates positively with complexity (near-decomposability hypothesis) in their skull, in an evolutionary context.

Phenotypic modules in skull networks
Phenotypic modularity varies among primate skulls in the number of modules (from 2 to 5), in the bones composing these modules, and, as a consequence, in the strength of modularity. Figs 3-6 show the modular organization of the skull of primates in a phylogenetic context.
We found four types of connectivity modules: (1) neurocranial, (2) midfacial, (3) palatal, and (4) premaxillary. The name of each module refers only to the anatomical region it occupies, and excludes any developmental, functional, and evolutionary interpretation. The neurocranial module is present in all skull networks and always groups together the occipital, the temporals, and the parietals; often, the neurocranial module also includes the sphenoid (in 18 Fig 1. Schema of the skull of primates formalized as a network. An anatomical network represents the bones and physical joints (i.e. sutures and synchondroses) of the skulls as nodes and links in a network model. Because these physical joints are primary sites of bone growth and diffusion of stress forces, topological relations also capture developmental and functional co-dependences among bones [27]. The analysis of these anatomical networks helps uncover the morphological organization of the skull regardless of variation in its shape and size. Labels: eth, ethmoid; fro, frontal; lac, lacrimal; max, maxilla; nas, nasal; nch, nasal concha; occ, occipital; pal, palatine; par, parietal; pmx, premaxilla; sph, sphenoid; tem, temporal; vom, vomer; zyg, zygomatic; l, left; r, right. out of 20 skulls) and/or the zygomatics (15), as well as the frontal/s in Tarsius and Gorilla, and the non-primate genera Tupaia and Cynocephalus. The midfacial module is present in all skull networks and always groups the ethmoid, the lacrimals, and (except in Aotus) the nasal conchae. Sometimes the midfacial module includes the frontal/s (14), the nasals (12), the palatines (6), the vomer (6), and the zygomatics (5). In Strepsirrhini (Lemur, Propithecus, Loris, Nycticebus) and in Gorilla the midfacial module splits into left and right midfacial modules. Here, unpaired bones are always included in the left module, but this is an arbitrary assignment because unpaired bones are equally connected to both sides of the skull. Exceptionally, the left and right palatines are grouped together in the left midfacial module of Loris and Nycticebus. Consequently, in these two taxa, left and right midfacial modules are asymmetrical because the palatines form a strong cluster that prevents their split in left and right modules. In fact, an independent palatal module groups the palatines in 14 out of 20 skull networks, and together with the vomer in 12 of them. The palatal module also includes the premaxillae in Callithrix, Saimiri, Macaca, Papio, and Hylobates; the maxillae in Aotus, Callithrix, Macaca, and Hylobates; the nasals in Callithrix and Macaca; and the sphenoid in Gorilla and Pan. When the palatines do not form an independent module, they form a block within the midfacial module together with the vomer. Finally, an independent premaxillary module is present in 9 out of 20 skull networks grouping the premaxillae; often this module includes the nasals (6) and, less frequently, the frontal in Pithecia and Colobus, the vomer and the maxillae in Lemur and Propithecus, and the nasal conchae in Aotus. In Aotus this premaxillary module splits into left and right premaxillary modules.   In contrast to Strepsirrhini, the frontal bone is unpaired in the skull of Haplorrhini (but see [33]). In Tarsius the frontal bone is included in the cranial module (in red), while in Platyrrhini the frontal belongs to the midfacial module (in blue) or to the premaxillary module (in yellow) in Colobus. Platyrrhini (right) show a conserved bone composition of the cranial module: occipital, sphenoid, parietals, temporals, and zygomatics. The three types of facial modules are also present, but they vary in their bone composition from one skull to another. The hierarchical grouping method that we used to identify connectivity modules also allows us to distinguish two main regions (or macromodules) in primate skulls according to the first split of the dendrogram (see SI), which we called facial and cranial because of their anatomical position. The cranial macromodule always comprises the neurocranial connectivity module, as well as the palatal module in Gorilla, Pan, Pithecia, and Pongo. The facial macromodule always comprises the midfacial module and, except for the aforementioned genera, the palatal module. The premaxillary connectivity module is also included within the facial module in all genera except for Pan and Tarsius, whose premaxillary module splits in the first place from both macromodules. Thus, we observe distinct palatal and premaxillary modules in the skull of primates in addition to the previously reported facial and cranial modules. The near-decomposability hypothesis The skulls of the genera analyzed show a narrow range of variation in complexity and modularity ( Table 1). The Abouheif's test indicates a significant phylogenetic signal in network parameters measuring complexity as number of parts (N: p = 0.003; K: p = 0.019). In contrast, network parameters measuring complexity as richness of interaction do not show a phylogenetic signal Within this group we observe the same variability achieved by the other skulls analyzed: one to three facial modules (in blue, green, and yellow), and variation in bones that integrate the cranial module (in red). The inclusion of the frontal, sphenoid, and zygomatics in the cranial module varies from one skull to another. The human skull shows a two-module division into facial and cranial modules: a unique organization within Primates, but similar to that of Tupaia. These results do not support the existence of a positive feedback loop between modularity and complexity during primate skull evolution, when we measured complexity as richness of interactions. Conversely, if we limit our definition of complexity to the number of bones, then our results support a positive feedback loop between modularity and complexity in this study; however, we view this result with some caution because of the narrow range of variation in the number of skull bones among primates.

Discussion
Our phylogenetic comparison of modularity in primates suggests a division of the skull into two phenotypic modules: one facial and one cranial. Previous studies have also identified these two regions as two separate connectivity modules in the adult human skull [27] and in newborn human skulls with various craniosynostosis conditions [28]. The division observed resembles the classic division of the mammalian skull into viscerocranium and neurocranium [34], although connectivity modules delimit different boundaries for these two regions. In primate skulls, we observe greater variability in the connectivity modules that form the facial module than in those that form the cranial module. For example, the facial region comprises up to three connectivity modules (midfacial, palatal, and premaxillary), which can be subsumed one into another, split in left and right specular modules, and/or group different bones. In contrast, the cranial region comprises invariably the neurocranial connectivity module, plus the palatal connectivity module in Gorilla, Pan, Pithecia, and Pongo. Phylogenetically, it is more parsimonious to postulate that this feature was acquired in great apes and then lost in humans (two steps; in addition to the independent acquisition in Pithecia) than that it was acquired in Gorilla, Pan and Pongo independently (3 steps). Therefore, this constitutes a synapomorphy of the Hominidae (Pongo, Gorilla, Pan and Homo), and might thus indicate that anatomical networks can also be useful for taxonomic and phylogenetic studies. The evolutionary and ecomorphological implications of this feature need to be studied in future works. In general, the cranial region shows only minor variations between major phylogenetic groups, such as Strepsirrhini, Platyrrhini, and Cercopithecidae (Figs 3 to 5). Exceptionally within Hominidae (Pongo, Gorilla, Pan and Homo), the cranial region varies in how it groups the frontal, sphenoid, and zygomatic bones (Fig 6). Previous studies showed that these bones act specifically connecting (i.e. integrating) the facial and the cranial regions in humans [27].  Connectivity modules derive from the analysis of the overall pattern of bone contacts in a skull, but they capture a deep organization of the growth patterns of the skull because bone contacts (i.e. craniofacial sutures) are primary sites of bone growth. As a consequence, the facial and cranial phenotypic modules identified have further developmental, variational, and functional implications for understanding the evolution of the human head [23]. Notably, the phenotypic modules defined in the present study discriminate the sutures of the skull according to the sequence of closure proposed for hominoids [35]. The cranial module includes the sutures that close earlier (vault, base, and circummeatal regions), whereas the facial region includes the sutures that close later (palatal, facial, and craniofacial regions). Interestingly, the patterns of sutures in the cranial region are conserved even in species as distantly related to humans as zebrafish [36]. This invariance of the suture patterns in the cranial region is surprising in the context of the tremendous increase in brain size, and associated changes in facial shape, during the evolution of primates [37], and suggests a deep conservation of skull morphogenesis across tetrapods.
Characterizing growth processes is essential to understand craniofacial shape variation during the evolution and development of the skull [4]. Thus, in addition to behaving as different units of allometric growth in humans [27], we would expect connectivity modules to be consistent with reports based on phenotypic variation in the size and shape of the skull. Our finding of a unitary cranial module agrees with the presence of an integrated cranial module (vault and base) in humans that shows patterns of covariation during postnatal ontogeny [38]. Along similar lines, a recent study by Villmoare et al. [39] has also reported the presence of premaxillary and palatal modules in the face of extant primates and extinct human relatives, similarly to what we found in this study analyzing connectivity patterns. Other studies have found that indeed the facial and cranial regions vary also in the degree of intra-module integration, which directly affects shape co-variation, functional coordination, and responses to selective pressures [6,40,41]. Moreover, Goswami and Polly [42] have demonstrated that the facial modules differ from cranial modules in their evolutionary plastic and morphological disparity in Primates. As a consequence, any change in the connectivity patterns of the skull, for example a fusion of frontal bones, might provoke specific shape changes in the entire skull. The premature fusion of bones (i.e. craniosynostosis) is a striking example of how a change in a single bone contact affects the shape of the skull, as a consequence of compensatory growth in its vicinity as well as in more distant parts of the skull [43][44][45].
Although deeply conserved developmental processes are likely to be responsible for much of the observed modularity in organisms, Wagner [20] argued that changes in modularity, particularly an increase in modularity, is likely to be the result of selection. Modularity allows more rapid change, and phenotypes that can change more rapidly in response to changing selection pressures have an adaptive advantage. Differences in modularity among primates are likely to be adaptive responses to differing selective pressures. In Platyrrhini, Marroig and Cheverud [46] found that patterns of morphological integration, while clearly influenced by phylogeny, also closely matched shared dietary patterns, a result supported by Makedonska and co-workers [47]. One result of our analysis (Fig 3) shows that, while conserved patterns are apparent in Strepsirrhini, lemurs have acquired a derived pattern. This is consistent with the relatively derived behavioral patterns seen in this clade. Similarly, in our results, within Catarrhinii (Fig 5), Colobus shows a relatively derived pattern of cranial modularity, which suggests that its unique dietary adaptation may influence patterns of modularity. The unique loss of the premaxilla in Homo, and its incorporation within the maxilla, is likely linked to adaptive changes relating to the anterior dentition. Overall, the conserved modularity of the neurocranium compared to the variable patterns of facial modularity seems to support adaptive causes for facial evolution in primates.
As Klingenberg [48] notes, however, changes in modularity and integration will alter the evolutionary constraints of a lineage. Certain phenotypes will no longer be accessible if modularity biases the directions of response to selection. Questions of evolvability and constraint are central to understanding how the varying influences of integration, modularity, and ontogeny interact under the pressure of natural selection to produce the resultant phenotype [49]. Biasing the directions in which responses to selection may operate-by changing the patterns of connectivity in the cranial skeleton-can have strongly advantageous consequences, if adaptive peaks in the fitness landscape fall along the available trajectories [50]. But if integration or modularity constrains the available variation, and peaks cannot be reached, organisms may not be able to respond before becoming extinct [51,52]. Primates are noteworthy for their adaptive diversity, and in future research we intend to investigate the relationship between adaptive diversity, morphological disparity, and modular variation, to determine if modular constraints are at play in other, less diverse clades, which are not present in Primates.
In conclusion, our results showed that the number of modules does not correlate positively with our complexity measures. In Simon's view, a system (here, the skull) is organized hierarchically when it is composed of interrelated subsystem (the modules) that are nested until reaching some lowest level of elementary subsystems (the bones). Taking number of modules as a proxy for the hierarchical arrangement of the skull (more modules = more hierarchy), Simon's near-decomposability hypothesis does not hold for the primate data. Instead, what we observed in Primates is an inverse relationship between these two properties: more complex skulls show lower modularity. This suggests that in primates the formation of new bone contacts, which foster complexity after fusion and loss of bones [30], also integrates modules and lessens the parcellation of the skull. However, we cannot neglect that the narrow range of values of the parameters used to quantify complexity in the skull of primates, if we compare it with the values of tetrapods [29], might conceal this correlation. A promising way of testing and further expanding the findings of this study is to analyze taxonomic groups of higher rank and with greater skull disparity such as the whole Tetrapoda clade. Our study paves the way to undertake such studies, as well as test and discuss other related hypotheses seeking to understand the evolution of phenotypic modularity and complexity in vertebrates.

Skull data and phylogeny
We studied the skull of 20 Euarchonta taxa: 18 primate taxa representing all the major extant groups of primates, and representatives of the two closest living relatives of primates: Dermoptera (Cynocephalus) and Scandentia (Tupaia) [53]. We coded for articulations among bones to model skull networks by inspecting natural skulls and casts housed in the zoology collection of the Dept. of Archaeology, Univ. of Saskatchewan. Data were collected from the following specimens: Aotus  Taxonomy ( . We calibrated the phylogeny of primates according to molecular estimations [32]; branch lengths indicate time of divergence in millions of years.

Modeling anatomical networks
Bone-bone articulations were coded in three-dimensional unweighted network models, in which network nodes represent skull bones and connections represent craniofacial sutures and synchondroses. For each skull, we coded 1/0 for the presence/absence of contact between two bones in a symmetric adjacency matrix of size equal to bone number. We performed the modeling and analysis of skull networks in R (R Core Team 2014) using the package igraph [54]. Adjacency matrices of skulls are available at http://dx.doi.org/10.6084/m9.figshare. 1221540. The protocol used to analyze each skull network is detailed in the S1 Protocol.

Identifying connectivity modules in anatomical networks
A connectivity module in a skull network is a group of bones with more contacts among them than to other bones outside the module [27]. We identified the number and bone composition of connectivity modules using a hierarchical cluster analysis of the pair-wise bone similarity, which is quantified as the generalized topological overlap (TO) between each pair of bones. TO is a normalized measure of similarity that quantifies pair-wise common neighbors between the nodes of a network: TO i,j = TO j,i = J(n i , n j )/min k (i, j), where J(n i ,n j ) is the total amount of neighbors in common between two nodes and min k (i,j) is the lowest connectivity of both nodes. TO ranges from 0 to 1: two nodes that connect to all the same other nodes have TO = 1, whereas two nodes without any neighboring node in common have TO = 0. This similarity index has been extensively used to analyze modularity in different types of networks [55,56]. Hierarchical clustering gathers together bones with higher TO into single branches until all bones form one group. For each potential partition we quantified Q: a quality index that quantifies how well a potential partition groups the nodes of the network compared to other possible partitions [57]: where K is the number of connections, A ij is the adjacency matrix, k i is the connections of i, k j that of j, m i is the module of i, and m j that of j. If m i = m j then δ(m i ,m j ) = 1, else δ (m i ,m j ) = 0. If the number of connections among nodes in the same module is not higher than expected at random then Q = 0, otherwise Q > 0: the higher the Q, the better the partition. For each skull network, the partition with the highest Q identifies its connectivity modules. We quantified the modularity of the skull network (M) as the product of the number of modules and Q. We identified connectivity modules using the functions hclust (package stats), GTOMmdist1 (labs.genetics.ucla.edu/horvath/GTOM/old/gtom.R), and modularity, package igraph.

Estimating complexity in anatomical networks
We estimated the complexity of skull using five network parameters: number of nodes (N), number of connections (K), density of connections (D), average clustering coefficient (C), and heterogeneity of connections (H). We chose these parameters as proxies to estimate diverse features that have been associated in the past with the morphological complexity of the skull [29]. N and K are the basic descriptors of the network: the number of nodes and the number of links. D is the number of existing node-node connections with respect to the total maximum possible according to the total number of nodes: D = 2k/N(N − 1). C is the average of the clustering coefficients of all nodes in the network: where τ i is the number of triangular motifs that include node i and k i is the number of connections of node i. H is the ratio between the variance and the mean of connectivity: H = σ K /μ K , where σ K is the variance of the number of connections of all nodes in the network and μ K is the mean of the number of connections.
Estimating the phylogenetic signal of modularity and complexity We tested the phylogenetic signal of the parameters estimating modularity (M) and complexity (N, K, D, C, H) using the Abouheif's test with 1000 permutations [58]. A significant phylogenetic signal would indicate a phylogenetically constrained evolution of modularity and/or complexity in Primates. We estimated the phylogenetic signal with the function abouheif.moran, package phytools [59]. We also performed a Blomberg's test (see SI).

Testing the near-decomposability hypothesis
According to Simon's near-decomposability hypothesis modularity and complexity evolve in a positive feedback. We tested the null hypothesis that there is no correlation between modularity (M) and complexity (N, K, D, C, H). Our alternative hypothesis was that modularity and complexity are positively correlated. We performed Pearson's product-moment correlations of phylogenetic independent contrasts using the method FIC4 [60]. We performed these tests with the function pic, package ape [61], and the function cor.test (package stats).
Supporting Information S1 Protocol. R Markdown file used to analyze the data. (HTML)