Developmental pathways inferred from modularity, morphological integration and fluctuating asymmetry patterns in the human face

Facial asymmetries are usually measured and interpreted as proxies to developmental noise. However, analyses focused on its developmental and genetic architecture are scarce. To advance on this topic, studies based on a comprehensive and simultaneous analysis of modularity, morphological integration and facial asymmetries including both phenotypic and genomic information are needed. Here we explore several modularity hypotheses on a sample of Latin American mestizos, in order to test if modularity and integration patterns differ across several genomic ancestry backgrounds. To do so, 4104 individuals were analyzed using 3D photogrammetry reconstructions and a set of 34 facial landmarks placed on each individual. We found a pattern of modularity and integration that is conserved across sub-samples differing in their genomic ancestry background. Specifically, a signal of modularity based on functional demands and organization of the face is regularly observed across the whole sample. Our results shed more light on previous evidence obtained from Genome Wide Association Studies performed on the same samples, indicating the action of different genomic regions contributing to the expression of the nose and mouth facial phenotypes. Our results also indicate that large samples including phenotypic and genomic metadata enable a better understanding of the developmental and genetic architecture of craniofacial phenotypes.

SCIENTIFIC REPORTS | (2018) 8:963 | DOI: 10.1038/s41598-018-19324-y components of total shape variation, and the signals of modularity in both morphospaces can be potentially informative of modularity acting at specific levels of organization. For instance, since fluctuating asymmetry has been postulated as an indicator of developmental instability (DI), the conservation or disruption of modularity and integration patterns across several genomic backgrounds can be informative of aspects of the developmental and genetic architecture of the human face.
The metaphor of paths of development 23 is useful for understanding the basis for the development of the covariance on morphological features. A development path denotes the set of processes that generate a trait 30 . Therefore, it is a term that incorporates molecular and cellular mechanisms underlying development processes, resulting in complex networks of interactions 30 . To characterize developmental integration/modularity it is important to distinguish the different origins of morphological covariation. A range of different mechanisms of development can produce interactions between developing traits and, therefore, generate covariation between them 1 . As Klingenberg 31 proposed, a way to identify developmental covariation or disruption simultaneously controlling for the effects of genetic and environmental covariation is to analyze patterns of fluctuating asymmetry (FA). Because the left and right sides of an individual share the same underlying genome, the study of FA is an effective way to detect developmental noise and minimize the effects of among-individuals genetic and environmental variation. In other words, covariation in the FA of different traits is due solely to the covariation of direct development pathways 28,31,32 . Therefore, it is possible to evaluate the role of direct developmental interactions that generate covariance by comparing between versus among-individuals FA covariation patterns. Thus, by quantifying FA one can test whether patterns of developmental interactions (observed in the asymmetric component) correspond to the main components of among-individual variation, as would be expected if the integration of development is a significant constraint on evolutionary variation 33 .
Here we aim to explore several modularity hypotheses in the human face, comparing their stability on the asymmetric and symmetric components of shape variation, and testing for covariation disruptions across several subsamples of Latin American mestizos differing in their Native American, European or African genomic ancestry. Specifically, we have used a composite sample of admixed Latin Americans including landmark data aimed to recover facial shape, along with genomic estimations of ancestry. Then, we have used a covariational statistical approach to test for modularity patterns in both, the symmetric and asymmetric component of facial shape variation, as a way to identify associations among developmental and among-individual variation. By organizing the original landmark configuration into several arrays, here we evaluate four general facial modularity hypotheses: a) functional modularity hypothesis (FMH), b) midline modularity hypothesis (MLMH), c) facial thirds modularity hypothesis (FTMH), and d) neurocranium-splachnocranium modularity hypothesis (NSMH) (see Fig. 1). Since we have recently reported that genetic ancestry is related to the patterns of directional and fluctuating asymmetry in the same sample of Latin-American mestizos, 34 we replicated all the analyses on ancestry-based subsamples, to see if the recent microevolutionary history of the studied groups affect the covariation patterns, or else remain stable across the different genetic backgrounds.
To further explore the modular structure of the face, we investigated the magnitude of integration of different facial structures. Theoretically, integration is stronger if all the variation is concentrated in a single dimension of the morphospace, indicating a perfect correlation of all traits, and is absent if the variation is uniformly distributed across all available dimensions. Accordingly, morphological integration in geometric morphometric data can be measured as the scaled variance of the eigenvalues (SVE) of a principal component analysis 35 , computed on the phenotype under study. When only a few eigenvalues are large relative to the rest, then the variance of all eigenvalues will be higher than if all eigenvalues explain similar amounts of variation. In the case that the variance is large, it is considered that the trait analyzed is strongly integrated, as variation is confined to a limited morphospace in relation to the theoretical total morphospace 36 . Hallgrimsson et al. 14 argued that the increase in the variance does not necessarily imply an increase in integration, and propose that to verify whether integration (as measured using SVE) is accompanied by an increase in the phenotypic variance (measured as the trace of the variance/covariance TVC) then the regression of SVE on TVC should be large and significant when the structure strongly covariates and phenotypic variance is high, as well. Thus, we complete our analyses by exploring putative disruptions in the internal integration of facial traits according to the Hallrimsson et al. proposal.

Results
An exploratory Principal Component Analysis of the symmetric and assymetric morphospaces is presented in Fig. 2. Morphings presented in the figure depict shape variation among the centroids of each ancestry subsample. Further details concerning shape variation in this sample can be consulted in refs 34,37 . In the symmetric morphospace, the first PC show facial shape changes associated with eyes placement in relation to other structures, being more inferiorly placed in the positive scores. Also, the positive values exhibit more anteriorly and superiorly placed mouths and chins. The perifrontal region is displaced forwardly in the positive scores, whereas the ears are placed slightly higher and posteriorly in the positive values. The second PC describes shape changes related to the eversion of the ears, along with a reduction of the distance between them. The asymmetric morphospace can be decomposed on a first PC describing changes associated to directional asymmetry towards the right side in the positive scores. The eyes seem to be the more symmetric structures.
Evaluations of the hypotheses of modular organization and the magnitude of integration between the subsets of landmarks was quantified as the covariance ratio 38 . For random sets of variables, the covariation ratio (CR) has an expected value of one. While CR values lower than one will indicate some degree of modularity within the structure, CR values higher than one will indicate greater covariation between regions than within them 38 . Results corresponding to the Adams' CR modularity index 38 indicated that all the modularity hypothesis are significant with p < 0.006 (Table 1). This indicates both, that modular organization is perceived as a patent process in the human face, and that the modular hypotheses analyzed here are, to some extent, coincident in their conception. However, specific particularities can be observed (Table 1), that inform about some modular patterns that seem to be supported by the CR coefficients as seen collectively. For instance, the "Functional" hypothesis, representing a modular organization inspired in the craneo-functional matrix theory and reflecting a key role of soft tissues in the modular formation of the face, exhibit the lowest CR values in both, the symmetric and the asymmetric components of shape, thus supporting a stronger modularity signal. Furthermore, this modularity signal is conserved across all the ancestry subsamples, suggesting that the genetic-phenotypic map intervening in the determination of adult facial phenotypes were not affected by the recent admixture dynamics. Klingenberg and McIntyre 39 proposed that if at given developmental process is operating on a particular morphology, the covariance matrices of the FA and individuals should be similar and proportional. The correlations between covariance matrices of FA and individual variation show, in general, a significant signal of proportionality among these matrices ( Table 2). A more detailed inspection of matrix correlations performed on the different subsamples indicate some degree of variation among the proportionality of the developmental noise, depicted by the FA covariance matrix, and the among-individuals genetic/environmental covariance matrices (Table S1). For instance, modules such as the mouth and nose tend to show low, non-significant correlations among the FA and the individual matrices, suggesting a disruption among developmental pathways and the genetic/environmental basis for such traits. Conversely, regions such as the neurocranium and the eyes tend to exhibit higher than the average and significant correlations, indicating an alignment across the genetic-development-phenotype map.
Results concerning the eigenvector variance (in both versions, scaled by total variance, and scaled by total variance and number of variables, Table 3) indicate that in both, the asymmetric and symmetric components, the more integrated landmark configurations are the mouth, the eyes, the sagittal plane, and the splachnocranium (see Supplementary Fig. S1). When such proxies to integration are corrected for total variance, as recommended by Hallgrimmson et al., the signal of integration is limited to specific cases. In fact, the regression of morphological integration, intended here as the eigenvectors variance scaled by the total variance (Evstv) on the variance of the structure, intended here as the trace of the covariation matrix, indicates a significant relationship only for the symmetric component in the total sample and the European subsamples (see Table 4 and Fig. 3). Besides significance, there are recurrent behaviors across the different subsamples in both asymmetric and symmetric variation. Note, for instance, that the landmark configurations depicting the sagittal plane, the ears, and the splachnocranium tend to be below of the confidence limit for the regression of morphological integration on variation. This indicates that these traits exhibit greater variation than expected for their degree of internal integration, and thus suggest a lower degree of canalization of such phenotypes (see Fig. 3).  Table 2. Values of correlation matrices of individual variation (Ind), fluctuating asymmetry (FA), and measurement error (error) for the total sample (see Table S1 for details of subsamples

Discussion
In this paper we evaluated the most classical modularity hypothesis in order to detect if one of them fits better the observed variation in two different morphospaces: the asymmetric and the symmetric components of shape variation. Also, we have explored if the detected patterns remain stable across different subsamples differing in their genomic ancestry. Our analyses are aimed to infer the relationship among diverging environmental or genetic perturbations assumed to be promoters of phenotypic variation, and to identify the perturbations due to putative changes in the development pathways underlying the facial phenotype 40 . In this context, integration/modularity processes acting at the developmental level can suggest an important role of development as a promoter of stability or disruption of covariation patterns 33 . Finally, we intended to test for differences in the integration/ variation ratio across different ancestry contexts and different facial anatomical structures, in order to corroborate if levels of developmental, genetic, or environmental signals are congruent or, conversely, show some degree of variation 28,31,32 . So far, studies exploring modularity patterns on the external, soft facial structures are absent or scarce. Classical studies of phenotypic modularity focused primarily on the analysis of functional or developmental skeletal units 41,42 . For example, a recent study has indicated that the semi-independency of modules is not only characterized by its spatial aspects, but also for its own temporary ontogenetic structure 43 , which reinforces the recognition of the integrative and complex role of development. Modularity patterns observed in our samples confirmed the intrinsic complexity of the human face suggested in previous reports.
All the modularity hypotheses obtained statistical support, in part because they overlap in terms of landmark allotting into similar structures, but also because of the complexity of the phenotype under study. As already outlined, the differences between the symmetric and asymmetric modular patterns may be related to differential responses to gene expression during the development of the left and right sides. Some authors have hypothesized that patterns of genetic modularity and development must evolve to match the patterns of functional requirements [44][45][46][47] . If so, patterns of modularity should be equivalent in both, the asymmetric and symmetric components. In the samples analyzed here, it seems clear that modularity at the asymmetric level, reflecting patterns of  covariation due exclusively to developmental factors, do not differ from the modularity pattern arising from the symmetric morphospace. Thus, a likely explanation to this observation is that developmental and genetic/environmental parcellation are somewhat aligned. Regarding the performance of the different modularity hypotheses, in both morphospaces the functional model performs better than the remaining modularity schemes, which indicates that modularity driven by organic functioning of the craniofacial complex is observed at the developmental level and maintained at the across-individuals adult phenotypes. This suggest that the modular pattern that operates within the asymmetric and symmetric morphospace are related to function, primarily (e.g. reflecting a relative independence of the eyes, nose, mouth and ears). Such functional patterning operating at the developmental scale have been reported previously 19,48 . Our results are in accordance with previous ones 15 in the sense that a traditional modular division into neurocranium, base and splanchnocranium 9,49 can be more complex than previously thought. The difficulty in exploring the relationship between craniofacial modules can be discussed from the perspective of the Palimpsest model, which holds that several covariation-generating factors influence each other over time, making the inference and detection of such factors from the mere observation of adult's phenotypes a difficult task 14 . In addition, Esteve-Altava et al. 50 , reported that modular patterns differ in relation to laterality, especially in the jawbone. In the same line, they define until four facial subunits with different modular variation in the sides of the face when compared to bony and muscular tissues 19 . However, note that the Esteve-Alta anatomical networks paper 12 deals with a notion of modularity that is not related to variational modularity and thus differs from our study of phenotypic covariance patterns.
Another remarkable point is that the abovementioned results remain stable across different subsamples diverging in terms of ancestry/heterozigosity. It is interesting to note that, even when we previously reported differences in the asymmetric facial phenotypes of the different ancestry sub-samples used here 34 the modular pattern seems to remain unaltered across all of them, reinforcing the idea that phenotypic evolution can be canalized throgouth stable patterns of covariation 51,52 . In this regard, González-José et al. 51 showed that inter-population differences in the matrices of variance/covariance are not associated with matrices of molecular or morphological distances among modern human populations, suggesting that the stability of integration is independent of the history and structure of populations. They also argued that a possible explanation for this result is that patterns of integration could be limited to intraspecific level and, in the skull, some speciation events may involve large rearrangements of integration patterns that would facilitate developmental scenarios on different regions of the skull 51 . The above assumes that morphological integration of functionally or developmentally related traits will be consistent and they will respond to selective pressures on a coordinated way 44,46,[53][54][55] .
To sum up, our results suggest that a modular organization directed by the functionally-based apportionment of covariation patterns can be seen as a regular pattern maintained across both, the genotype-phenotype map, and across several genomic backgrounds.
That the modularity/integration patterns are stable across sub-samples with different genetic ancestry is not that surprising. Previous research has already reported stability when considering different human populations (see González-José et al. 51 ), even when individuals with a disease-altered facial development are considered. For instance, Martínez-Abadías et al. 56 suggested that FGF/FGFR signaling is a covariance-generating mechanism established early in vertebrate evolution that acts as a global factor modulating coordinated development of various skull components. Thus, this as well as other unknown global genetic and non-genetic factors can be seen as covariance-generating processes that constrain both, normal or disease-altered craniofacial phenotypes Martínez-Abadías et al. 57 . Recent research has also clarified that the signaling activity of mutations at systems such as the FGF/FGFR are not restricted to skeletal tissues, but also occurs at cellular levels, where signaling may trigger specific patterns of gene expression and/or cell/tissue differentiation Martínez-Abadías et al. (2013) 58 .
Regarding the stasis/disruption of covariation patterns, Klingenberg 21,28,31 argued that if there is a significant correlation between the covariance matrices of individual variation and fluctuating asymmetry, then it is likely that both were modeled similarly during development. Thus, any difference between the individual and FA patterns would be the result of a chronological divergence between the times in which each morphospace is modeled during development. In this context, our results suggest that the same processes are involved in the generation of the observed variation in the two levels. In other words, genetic and environmental factors that cause variation between individuals are suspected to produce patterns of change that are similar to the processes acting at the within individuals level, creating differences between the sides of the body.
An inspection of the literature on this topic appears as inconclusive regarding the congruence among the individual and the FA covariation matrices. For instance, Debat et al. reported a null congruence in the case of the mouse skull 56 , and suggest that it could operate some buffering against the impact of different sources of variation operating on distant regions such as the nasal capsule, the zygomatic arch, the orbit, and the cranial vault, etc. Similarly, Klingenberg et al. 59 found no congruence between individual covariance matrices and FA in pharyngeal jaws of cichlid fishes, which is associated with a polymorphism identified as a source of discrepancy. Similar tests performed on insects, however, report congruence among both covariation matrices (see refs 39,60 for flies, and ref. 61 for bumblebees). Several factors can be listed to explain these differences. One is the obvious discrepancy between genetic and developmental pathways in insects, mice, fish and humans. Phylogenetic analyses at the macro-evolutionary level are then necessary to answer these questions on an evolutionary context. Moreover, Debat et al. 56 suggest that natural selection acting on the shape of the wings of insects could dilute or hide the differences between within and among-individuals developmental processes. In the case of the human face, selection processes have also been postulated 62,63 , thus part of its expression pattern could be mediated by past evolutionary events.
A more detailed inspection of the discrepancies among FA and individual covariance matrices (Table S1) enabled us to investigate more fine differences among covariance patterns. Note, for instance, that the eyes, ears and neurocranium modules would be affected by parallel or equivalent noise-stability development processes. The complexity of the mosaic expression of developmental pathways is suggested indirectly on such results: the eye region develops around day 22, with the appearance of a pair of shallow grooves on both sides of the forebrain 64 . In addition, the medial migration of the eyes from their initial side locations results on the massive growth of the cerebral hemispheres and enlargement of the head, and enable the normal movement of eyes toward the sagittal line 65 . The largest migratory eye movements occur between the fifth and ninth week of the development. Thereafter, the angulation of the optical axes (between 71°-68°) is stabilized during postnatal phases 65  Regarding ears' development, they begin its formation around day 22 and its development extends to week six. Specifically, their tissues derive from ectoderm precursors. The pavilion and the ear meatus move from the base of the neck (cervical region) to its normal adult location on the side of the head. This process is largely due to mandibular growth. On a recent GWAS analysis, we demonstrated that there are seven genomic regions significantly associated with the size and the attachment of the lobe, folding of antihelix, rotation of the helix, the protrusion of the ear and the antitragus size 66 . These features are associated with variants of EDAR, TBX15 and CART1 genes. In summary, the parallelism between individual and FA covariation patterns seems to be supported by the synchronic timing of developmental events in the growing embryo, and some candidate genes have been detected whose role must be further investigated.
In contrast, the mouth and the nose arise from a rather complex combination of developmental processes, i.e.: they condense the combined action of several mechanisms, tissue origins and differentiation types. Notably, they come from the same layer, the endoderm, forming during the middle of week three and keeping their development until birth. During the half of the sixth week only the oronasal membrane divides the oral and nasal cavities. The mouth results from the interrelation of the maxillary and mandibular prominence, and besides depends on the closure of the palate, the vomeronasal system, the configuration of the nose and in general of all the pharyngeal cavity 64 . The complex processes described above may explain the discrepancy (or lack of proportionality) of the covariance matrices. In another recent GWAS study, four genomic regions (4q31, 6p21, 7p13, and 20p11) were significantly associated to features related to the nose morphology: the inclination of the columella, the amplitude of the nasal bridge, and the alar amplitude of the nose (with p-values = 3 × 10.9 to 9 × 10-9). The reported SNPs are associated to DCHS2, RUNX2 and GLI3 genes, respectively, while the region 20p11 overlaps with PAX1 67 . Thus, the hits in different and non-overlapping genetic systems also correlates with different phenotypes: the ear structures nucleated in the pinna 66 , versus localized nasal traits 67 . In some way, these combinated results can be seen as a preliminary signal of genetic modularity, detected by different methods as the ones used here.
Note that mouth and nose, as already mentioned, are structures where several functional, developmental, and growth processes simultaneously converge. Maybe, this complexity trigger the observed differences in the covariance patterns observed on eyes, ears and neurocranium. A point to note is that all the comparisons made in the middle third, which include three landmarks associated with the neurocranium, nose, and mouth (Figs 1 and 3) were statistically significant. This may endorse the reported influence of neurocranium in shaping splachnocranium 49 .
Moreover, natural selection acting quite independently upon the modules 56 , the complexity of structures 68 , and the phenotypic plasticity of the nose and mouth, as a whole, can be a plausible explanation to the observed differences in terms of matrix proportionality. For example, the maxillary region exhibits large craniofacial plasticity 69 , maybe due to the fact that the masticatory apparatus grows even up to 25 years 65,70 . Similarly, in the case of the nose, there is evidence of plasticity associated with the climate and it is possibly a structure that evolved under a regime of natural selection to cold climates [71][72][73] . Considering all the above, it is likely that the mouth and nose do not achieve the integration levels exhibited by other facial modules.
Most, if not all evidence provided so far by previous studies 15,51,52,56,[74][75][76] focuses on the symmetric component form. Instead, our approach was twofold, focusing also in the asymmetric dimension of variation. In this regard, a noticeable result is that the eyes, mouth, and ears tend to present larger-than-average morphological integration (Table 3). Previous asymmetry studies suggested that integration at the development level can play an important role in canalization of variation through periods of environmental change. Such statements seems to be supported by a study of skulls of Late Pleistocene carnivores and the relationship between fluctuating asymmetry, phenotypic integration, variance and environmental change or stress 33 , and similar results are corroborated in mandibles of shrews (for details see refs 77,78 ). In the abovementioned examples, the effects of environmental stress artificially induced the integration of development and variation in the jaws of shrew, and also showed that fluctuating asymmetry and the variation of the stressed population increased and that this increasing was canalized along the same direction of between-species variation 77,78 . In this context, we have obtained a recurrent correlation between the pattern of fluctuating asymmetry and morphological integration across all the subsamples in the symmetric component of shape, which is not the case when the asymmetric morphospace is explored for such correlation (Table 5). In this sense the regression trends corresponding to each genomic subsample are roughly equivalent (i.e. are positive), excepting for the European descent sample, that shows a pattern of negative correlation between asymmetry and integration in the asymmetric morphospace (Table 5). This suggest that the relationship between asymmetry and integration differs at the symmetric and asymmetric morphospaces. Theoretically, we could infer that the developmental pathways operating on individuals of (mostly) European-descent are divergent, to some extent, a statement that supported also by their differential expression on the pattern of facial asymmetries that we have reported recently 34 . Since functional structures are expressed early during prenatal development, it is expected that developmental modularity can potentially influence their functional counterpart, although more research is needed to disentangle the robustness of the relationship between these two modularity types 79 . Of course, the inverse explanation cannot be discarded, and it could be the case that functional modularity influences developmental processes such as bone remodeling and other forms of plasticity in which the loading mechanisms influence the rates and direction of tissue growth 80,81 , thus promoting modularity at the developmental level. While previous analyses have hypothesized that selection makes functional and genetic modularity converge 22,44,45,82,83 , other studies indicate that different morphological structures can perform equivalent functions and, therefore, there can be considerable flexibility in the neutral divergence 84,85 or among-trait connectivity, independent of their function 19,48 . Whatever the case, our results indicate that the modular landscape of the SCIENTIFIC REPORTS | (2018) 8:963 | DOI:10.1038/s41598-018-19324-y human face can be affected by both, genetic ancestry and the morphospace under study (asymmetric versus symmetric). Furthermore, these results are of importance to future Genome Wide Association Studies based on CANDELA (and other) datasets, since they provide a first identification of integration and modularity at the phenotypic scale, that can be used to submit more refined traits to the association analyses, as a preferred approach instead of heuristic searches based on analyzing non-supervised, raw characters.

Materials and Methods
The sample. As part of the CANDELA initiative, we recruited 4,104 volunteers from six Latin-American cities: Mexico City (Mexico), Medellin (Colombia), Lima (Perú), Arica (Chile), Porto Alegre, and Jequié (Brazil). The CANDELA consortium aims to evaluate the genetic basis of nonpathological phenotypes differentiated between European, American, and African populations through the analysis of admixed populations (see ref. 37 ). Volunteers with antecedents of craniofacial dysmorphologies, orthodontics treatments or severe facial trauma were not considered in this study. Further sample details are provided in Supplementary Table S3. Approvals provided by the Ethics Committees of the Universidad Nacional Autónoma de México and Escuela Nacional de Antropología e Historia (México), Universidad de Antioquia (Colombia), Universidad Peruana Cayetano Heredia (Perú), Universidad de Tarapacá (Chile), Universidade Federal do Rio Grande do Sul/Universidade Estadual do Sudoeste da Bahia (Brazil), and University College London (UK) were obtained prior the data collection, and an informed consent were signed by each participant before genetic, socioeconomic, and facial phenotypes data was collected. The same sample was used on a recent contribution 34,37,86 and conforms the main database of the CANDELA consortium 34,37,66,67,86,87 . All methods and procedures used here were performed in accordance with relevant guidelines and regulations (see below). Validation dataset are available to download from https://laofunam.com/data.
Facial shape was captured following scientific photographic protocols described in detail in references 34 and 86 . Upon these images, two observers (MQS and LC) placed a set of 34 standard facial landmarks using Photomodeler (http://www.photomodeler.com/ Eos Systems Inc, Vancouver, Canada). As described elsewhere 34,86 , this platform corrects for any lens distortion automatically, and we have followed the standard recommendations for quality and accuracy provided by the software. A scale factor was assessed using the nasion-gnathion distance measured directly on the individuals using a standard anthropometric caliper.
On each individual, blood samples were collected and DNA extraction was performed following standard laboratory procedures. Genomic data involving 730,525 marker SNPs was obtained from these samples (see further details in ref. 37 ). The SNPs were pruned to remove Linkage Disequilibrium, and after removing correlated SNPs, 90,000 SNPs were left for analysis. Ancestry estimation was performed with this SNP data. Genome-wide average heterozygosity was estimated from this data using PLiNK 88,89 , which provides a measure of excess heterozygosity compared to the overall sample. It is calculated as 1-excess homozygosity, while excess homozygosity is estimated using the inbreeding coefficient as the average excess of homozygous alleles across all SNPs for an individual as compared to the overall sample. To enable among-subsample comparisons, individuals were allotted to American, European, or Admixed groups if their population-specific ancestry markers were above 80% in the first two cases, and below that percentages in the third group (admixed).
Preliminary statistical analyses: error measurement and acquisition of the symmetric and asymmetric components of variation. Since on a previous paper 86 we have demonstrated that measurement error in size and shape was significantly lower than variation in FA, and thus negligible, the rest of the analyses were based on a single digitization of landmarks.
For the multilevel analysis of morphological integration and modularity, we first decomposed the total shape variation into its symmetric and asymmetric components by conducting a Procrustes ANOVA on the total dataset. The symmetric component is the variation between individuals in terms of the averages of the left and right configurations and corresponds to the phenotypic variation 59 . Conversely, the asymmetric component is the variation within individuals in terms of differences between configurations from the left and the right sides of each individual and depicts variation arising from direct developmental interactions 59 . Main shape changes in the symmetric and asymmetric morphospaces were preliminary explored using Principal Component Analysis and morphing of a scanned face to depict shape variation among the centroid of the of the genetic ancestry subsamples, via Discriminant Function Analysis.
Modularity hypotheses. The FMH (Fig. 1a) is based on the craneo-functional matrix theory, and is aimed to represent a key role of soft tissues in the modular formation of the face which would facilitate its stability, performance and evolution 90,91 . Moss and Young 91 argued that the head is divided into functional components (e.g. modules) determined by the soft tissues and cranial cavities, and limited by their bony surroundings. In this regard, the soft tissue components would guide the development of skeletal units. It should be noted that the functional components are based on assumptions around form and function links rather than the results of quantitative analyses or empirical evidences 19 . MLMH (Fig. 1b) is aimed to explore if the farther a module is of the sagittal midline, the more asymmetry will display 92 . Recently, it has been suggested that the expression of FA not only depends on developmental stability, but also on the cost of growth of the trait, defined as the amount of structural components necessary to form a unit of length of a given character 68 . In accordance with this argument, a trait with more structural components per unit of length should show lower asymmetry than a simpler one. In other words, this hypothesis seeks to verify whether those landmarks that are located closer to the sagittal plane versus those placed on the sides follow a modular pattern.
The FTMH (Fig. 1c) is aimed to test a traditional classification of facial anatomy, mainly used in medicine studies [93][94][95] . It is important to check how much these regions respond to a modularity pattern with empirical support.
Finally, the NSMH (Fig. 1d) is an attempt to verify if the differences in development timing, tissular precursors, and epigenetic stimuli that the splachnocranium (face) and the neurocranium 41 experience are enough to generate a modularity pattern observed in the soft structures studied here.
All the above-mentioned hypotheses will be tested in the whole sample as well as in sub-samples organized according to their genomic ancestry. Also, the hypotheses are evaluated in both the symmetric and the asymmetric components of shape variation.
The CR coefficient is a ratio of the covariation between modules over the covariation within modules, and consequently it ranges between zero and positive values 38 . For random sets of variables, the CR has an expected value of one. Whereas CR values lower than indicate some degree of modularity within the structure, CR values higher than one depict greater covariation between regions than within them 38 . Thus, we calculated the CR values on the several modularity hypotheses tested here, and across the different ancestry subsamples. The significance values were obtained by comparing the observed CR values with permutational distributions of 999 CR values obtained by assigning the landmarks randomly to modules. The proportion of permuted values lower than the observed CR value was used as the significance of the test 38 .
Covariance matrix similarity: exploring developmental pathways. Similarity between FA and individuals covariance matrices was tested following Klingenberg and McIntyre 39 . Similarity of the covariance matrices resulting from the Procrustes ANOVA was tested by computing pairwise matrix correlations, as a way to compare FA versus individual covariance patterns. Matrix correlation is a measure of the overall similarity of covariance matrices and its use is ubiquitous in geometric morphometrics 96,97 . Statistical significances were determined through matrix permutation tests, with 10,000 iterations, against the null hypothesis of complete dissimilarity between the covariance matrices concerned 39 .

Facial morphological integration.
Theoretically, integration is stronger if all the variation is concentrated in a single dimension, indicating a perfect correlation of all traits, and is absent if the variation is uniformly distributed across all available dimensions. Accordingly, morphological integration in geometric morphometric data can be measured as the scaled variance of the eigenvalues (SVE) of a principal component analysis 35 , made on the structure under study. When only a few eigenvalues are large relative to the rest, then the variance of all eigenvalues will be higher than if all eigenvalues explain similar amounts of variation. In the case that the variance is large, it is considered that the trait analyzed is strongly integrated, as variation is confined to a limited morphospace in relation to the theoretical total morphospace 36 . Hallgrimsson et al. 14 argued that the increase in the variance does not necessarily imply an increase in integration, and propose that to verify whether integration (as measured using SVE) is accompanied by an increase in the phenotypic variance (measured as the trace of the variance/covariance TVC) then the regression of SVE on TVC should be large and significant when the structure strongly covariates and phenotypic variance is high, as well. Analyses were performed on both, the symmetric and asymmetric component 39 .