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

Mirsha Quinto-Sánchez, Francesc Muñoz-Muñoz, Jorge Gomez-Valdes, Celia Cintas, Pablo Navarro, Caio Cesar Silva De Cerqueira, Carolina Paschetta, Soledad de Azevedo, Virginia Ramallo, Victor Acuña-Alonzo, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Mirsha Quinto-Sánchez, Francesc Muñoz-Muñoz, Jorge Gomez-Valdes, Celia Cintas, Pablo Navarro, et al.. Developmental pathways inferred from modularity, morphological integration and fluctuating asymmetry patterns in the human face. Scientific Reports, 2018, 8 (1), 10.1038/s41598-018-19324-y. hal-03268550

HAL Id: hal-03268550

https://hal.science/hal-03268550v1

Submitted on 5 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

OPEN

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

Received: 27 September 2017

Accepted: 15 December 2017

Published online: 17 January 2018

Mirsha Quinto-Sánchez<sup>1,2</sup>, Francesc Muñoz-Muñoz<sup>3</sup>, Jorge Gomez-Valdes<sup>1,4</sup>, Celia Cintas<sup>2</sup>, Pablo Navarro<sup>2</sup>, Caio Cesar Silva de Cerqueira<sup>5</sup>, Carolina Paschetta<sup>2</sup>, Soledad de Azevedo<sup>2</sup>, Virginia Ramallo<sup>2</sup>, Victor Acuña-Alonzo<sup>6,7</sup>, Kaustubh Adhikari<sup>6</sup>, Macarena Fuentes-Guajardo<sup>6,8</sup>, Tábita Hünemeier<sup>9</sup>, Paola Everardo<sup>7,10</sup>, Francisco de Avila<sup>7</sup>, Claudia Jaramillo<sup>11</sup>, Williams Arias<sup>11</sup>, Carla Gallo<sup>12</sup>, Giovani Poletti<sup>12</sup>, Gabriel Bedoya<sup>11</sup>, Maria Cátira Bortolini<sup>13</sup>, Samuel Canizales-Quinteros<sup>14</sup>, Francisco Rothhammer<sup>15</sup>, Javier Rosique<sup>16</sup>, Andres Ruiz-Linares<sup>6,17,18</sup> & Rolando Gonzalez-Jose<sup>2</sup>

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

<footnote n="1"><sup>1</sup>Ciencia Forense, Facultad de Medicina, Universidad Nacional Autónoma de México, Ciudad de México, Mexico.</footnote>

<footnote n="2"><sup>2</sup>Instituto Patagónico de Ciencias Sociales y Humanas. Centro Nacional Patagónico, CONICET, Puerto Madryn, Argentina.
<sup>3</sup>Departament de Biologia Animal, de Biologia Vegetal i d'Ecologia, Facultat de Biociències, Universitat Autònoma de Barcelona, Avinguda de l'Eix Central, Edifici C, E-08193, Bellaterra (Cerdanyola del Vallès), Spain.</footnote>

<footnote n="4"><sup>4</sup>Posgrado en Antropología Física, Escuela Nacional de Antropología e Historia, Ciudad de México, Mexico.</footnote>

<footnote n="5"><sup>5</sup>Superintendência da Polícia Técnico-Científica do Estado de São Paulo. Equipe de Perícias Criminalísticas de Ourinhos, São Paulo, Brazil.
<sup>6</sup>Department of Genetics, Evolution and Environment, and UCL Genetics Institute, University College London, London, UK.
<sup>7</sup>Licenciatura en Antropología Física, Escuela Nacional de Antropología e Historia, Ciudad de México, Mexico.
<sup>8</sup>Departamento de Tecnología Médica, Facultad de Ciencias de la Salud, Universidad de Tarapacá, Arica, Chile.
<sup>9</sup>Departamento de Genética e Biología Evolutiva, Instituto de Biociências, Universidade de São Paulo, São Paulo, Brazil.
<sup>10</sup>Posgrado en Antropología, Instituto de Investigaciones Antropológicas, Universidad Nacional Autónoma de México, Ciudad de México, Mexico.
<sup>11</sup>Universidad de Antioquia, Medellín, Colombia.
<sup>12</sup>Laboratorios de Investigación y Desarrollo, Facultad de Ciencias y Filosofía, Universidad Peruana Cayetano Heredia, Lima, Peru.
<sup>13</sup>Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Porto Alegre, Brazil.
<sup>14</sup>Facultad de Química, UNAM, Mexico City, Mexico.
<sup>15</sup>Facultad Instituto de Alta Investigación Universidad de Tarapacá, Programa de Genética Humana ICBM Facultad de Medicina Universidad de Chile y Centro de Investigaciones del Hombre en el Desierto, Arica, Chile.
<sup>16</sup>Departamento de Antropología, Facultad de Ciencias Sociales y Humanas. Universidad de Antioquia, Medellín, Colombia.
<sup>17</sup>MOE Key Laboratory of Contemporary Anthropology, Fudan University, Shanghai, China.
<sup>18</sup>Aix Marseille Univ, CNRS, EFS, ADES, Marseille, France. Correspondence and requests for materials should be addressed to R.G.-J. (email: rolando@cenpat-conicet.gob.ar)</footnote>
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.

Modularity is an important principle of organization in biological systems, which is also manifest at the morphological level¹. Modularity means that, on a complex structure, some set of traits exhibit greater internal integration in relation to the integration with traits belonging to other modules²,³. Modules observed in complex structures, such as the human face, consequently exhibit some degree of independence, according to their development, evolutionary, population-specific or experimentally induced conditions⁴.

The human face is a bilateral complex phenotype composed of a combination of structures that vary simul-
taneously in an integrated and modular way. For example, it has been reported that as part of the skull, the face
is constrained in some way by variations in the skull base⁵⁻⁹. Conversely, other authors have suggested that it is
closely linked to changes occurring in the neurocranium during brain growth¹⁰. This link has been evidenced in
several experiments made on mouse models. Martínez-Abadías *et al.*¹¹, for instance, analyzed integration of brain
and skull in newborn mouse models of Apert syndromes and demonstrated that skull dysmorphology affects not
only the coronal suture closure, but also the facial skeleton. The authors proposed that shared phenogenetic pro-
cesses, such as regional differentiation by dynamic inductive signaling and repetitive patterning by quantitative
interactions affecting several tissues during head morphogenesis, underlie not only premature suture closure, but
additional cranial anomalies in Apert syndrome, including facial traits. Hallgrímsson *et al.*¹²⁻¹⁴ analyzed pertur-
bations to craniofacial development in mouse models to disentangle the processes that determine covariation
structure. These studies have shown that mutations that influence the growth of the chondrocranium produce
a common pattern of integrated shape change throughout the mouse skull. Furthermore, Martínez-Abadías
*et al.*¹⁵,¹⁶ have shown that this same integrated pattern is also a main component of covariation structure in the
human skull. In other words, shape changes observed in the integrated effects of variation in chondrocranial
growth in mouse mutants also correspond to a significant axis of covariation in humans.

Viewed from a different perspective, an analysis of simulated amniote developmental facial morphospaces
suggest that epigenetic factors such as organismal geometry and shape impact facial morphogenesis and alter
the locus of adaptive selection to variation in later developmental events¹⁷. More generally, the face can be repre-
sented after the Palimpsest Model¹³,¹⁴, that propose that patterns of covariation observed in the adult phenotype
emerge from different processes of variance generation, that gradually overlap, add or integrate sequentially dur-
ing the individual's ontogeny¹⁸. Such particularities, then, define the face as a complex phenotype. However, the
modular nature of structures into the face has been scarcely analyzed, to the exception of some recent approaches
to the craniofacial modular patterns using network theory¹⁹.

Modularity can be studied at four different, non-exclusive levels: developmental, genetic, functional and evo-
lutionary²⁰⁻²². Developmental modules are defined by interactions between precursors that ultimately participate
during the development of the adult structure¹². This includes a wide range of processes, such as developmental
switches leading track bifurcations of development, or signaling between tissues through different molecular
mechanisms²³. Development processes can influence each other and therefore achieve a coordinated development
of tissues, organs and the whole body, promoting the so called developmental stability. Also, the relationship
between development stability-instability can mediate the expression of genetic and environmental variation by
transmitting their effects through different routes in development. Therefore, developmental modularity contrib-
utes to all components of phenotypic covariance between traits²⁴.

Functional modularity, by contrast refers to interactions among parts of the organism performing a given function¹. The typical example in the human face includes the array of muscles and bones that perform mastication, which generates intense strains able to modify the osseous phenotypes by, for instance, bone remodeling¹,²⁵.

Genetic modularity refers to patterns of joint effects between genes on traits that can be represented as a
network of relationships between pleiotropic traits²⁶. This type of modularity needs to be placed in the context of
the "genotype-phenotype map"¹,²⁷. Because developmental processes mediate the expression of genetic variation
in phenotypic traits, developmental and genetic modularity are related. As stated by Klingenberg¹, this relation-
ship need not be a perfect congruence, however, because the expression of genetic variation is not exclusively
controlled by developmental interactions¹. Moreover, genetic changes can influence developmental modularity
by causing alterations in the interactions among the developmental pathways that affect the traits of interest¹,²⁸.

Finally, evolutionary modularity is intended as the coordinated evolutionary divergence in different traits. As genetic variation is a critical determinant for evolutionary change by selection and drift²⁹, genetic parcellation contributes substantially to evolutionary modularity. Functional modularity is also an important determinant of evolutionary diversification because it provides a link between the modular structure of morphological traits and selection on performance in organismal functions¹.

As observed, the four levels of modularity (and its counterpart, integration) are interrelated in complex ways.
In many occasions, a common way to infer modularity at different levels departs from the analysis of phenotypic
covariation patterns. As explained by Klingenberg²⁴, identifying the modular components of a configuration of
landmarks is an important task of morphometric analyses in evolutionary developmental biology. Since traits
within modules are tightly correlated with each other, but relatively independent of traits in other modules, then
hypotheses concerning the boundaries of modules in a landmark configuration can be tested by comparing the
strength of covariation among alternative partitions of the configuration into subsets of landmarks. If a subdivi-
sion coincides with the true boundaries between modules, the correlations among subsets should be minimal²⁴.
Furthermore, this variational modularity concept can be applied to both, the symmetric and the asymmetric
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²³ 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³⁰. Therefore, it is a term that incorporates molecular and cellular mechanisms underlying development processes, resulting in complex networks of interactions³⁰. 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¹. As Klingenberg³¹ 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²⁸,³¹,³². 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³³.

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,³⁴ 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³⁵, 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³⁶. Hallgrimsson et al.¹⁴ 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 Hallgrimsson et al. proposal.

## Results

An exploratory Principal Component Analysis of the symmetric and asymmetric 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³⁴,³⁷. 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³⁸. 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³⁸. Results corresponding to the Adams' CR modularity index³⁸ 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.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Anatomical location of the landmarks used in this study and modularity hypotheses studied here. (a) Functional modularity hypothesis (FMH): eyes = orange, mouth = blue, purple = ears, nose = green. (b) Midline modularity hypothesis (MMH): sagittal = black, midsagittal = green-olive, lateral = yellow. (c) Facial thirds modularity hypothesis (FTMH): middle = blue-turquoise, inferior = white. (d) neurocranium-splanchnocranium modularity hypothesis (NSMH): neurocranium = blue-dark, splanchnocranium = red.</figcaption></figure>

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³⁹ 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
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Principal Components Analysis. Scatterplots and morphings of the symmetric (a) and asymmetric (c) component of shape variation. Ellipses represent the 90% of the variance and colored by genetic ancestry (afr = red circule, ame = blue triangle, eur = green triangle, adx = purple square). (b and d) Morphed heads representing the mean shape observed on the four ancestry groups computed through discriminant analysis.</figcaption></figure>

(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.
<table><caption>**Table 1.** Adams's modularity test for the total sample and subsamples, modularity hypotheses and shape components. Columns presents covariance ratio (CR) and *p*-value. Bold and italic indicate the lowest CR value for asymmetric and symmetric component of the shape variation, respectively.</caption><thead><tr><th>Source</th><th>Module</th><th>Component</th><th>CR</th><th>P-value</th></tr></thead><tbody><tr><td rowspan="8">Total</td><td rowspan="2">Functional</td><td>Asymmetric</td><td><strong>0.4841</strong></td><td>0.001</td></tr><tr><td>Symmetric</td><td><strong>0.6148</strong></td><td>0.001</td></tr><tr><td rowspan="2">Facial third</td><td>Asymmetric</td><td>0.5030</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.8151</td><td>0.005</td></tr><tr><td rowspan="2">Middle line</td><td>Asymmetric</td><td>0.5320</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.6559</td><td>0.001</td></tr><tr><td rowspan="2">Neuro-face</td><td>Asymmetric</td><td>0.7145</td><td>0.003</td></tr><tr><td>Symmetric</td><td>0.7531</td><td>0.001</td></tr><tr><td rowspan="8">American</td><td rowspan="2">Functional</td><td>Asymmetric</td><td><strong>0.5534</strong></td><td>0.001</td></tr><tr><td>Symmetric</td><td><strong>0.6498</strong></td><td>0.001</td></tr><tr><td rowspan="2">Facial third</td><td>Asymmetric</td><td>0.5778</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.8542</td><td>0.006</td></tr><tr><td rowspan="2">Middle line</td><td>Asymmetric</td><td>0.6052</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.6901</td><td>0.001</td></tr><tr><td rowspan="2">Neuro-face</td><td>Asymmetric</td><td>0.8006</td><td>0.004</td></tr><tr><td>Symmetric</td><td>0.7651</td><td>0.001</td></tr><tr><td rowspan="8">European</td><td rowspan="2">Functional</td><td>Asymmetric</td><td><strong>0.4997</strong></td><td>0.001</td></tr><tr><td>Symmetric</td><td><strong>0.6110</strong></td><td>0.001</td></tr><tr><td rowspan="2">Facial third</td><td>Asymmetric</td><td>0.5264</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.8155</td><td>0.001</td></tr><tr><td rowspan="2">Middle line</td><td>Asymmetric</td><td>0.5265</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.6572</td><td>0.001</td></tr><tr><td rowspan="2">Neuro-face</td><td>Asymmetric</td><td>0.6890</td><td>0.002</td></tr><tr><td>Symmetric</td><td>0.7557</td><td>0.001</td></tr><tr><td rowspan="8">Admixed (Heterozygous)</td><td rowspan="2">Functional</td><td>Asymmetric</td><td><strong>0.4903</strong></td><td>0.001</td></tr><tr><td>Symmetric</td><td><strong>0.6238</strong></td><td>0.001</td></tr><tr><td rowspan="2">Facial third</td><td>Asymmetric</td><td>0.5039</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.8388</td><td>0.005</td></tr><tr><td rowspan="2">Middle line</td><td>Asymmetric</td><td>0.5411</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.6670</td><td>0.001</td></tr><tr><td rowspan="2">Neuro-face</td><td>Asymmetric</td><td>0.7292</td><td>0.003</td></tr><tr><td>Symmetric</td><td>0.7532</td><td>0.001</td></tr><tr><td rowspan="8">Homozygous</td><td rowspan="2">Functional</td><td>Asymmetric</td><td><strong>0.5084</strong></td><td>0.001</td></tr><tr><td>Symmetric</td><td><strong>0.6152</strong></td><td>0.001</td></tr><tr><td rowspan="2">Facial third</td><td>Asymmetric</td><td>0.5316</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.8432</td><td>0.006</td></tr><tr><td rowspan="2">Middle line</td><td>Asymmetric</td><td>0.5440</td><td>0.001</td></tr><tr><td>Symmetric</td><td>0.6561</td><td>0.001</td></tr><tr><td rowspan="2">Neuro-face</td><td>Asymmetric</td><td>0.7286</td><td>0.003</td></tr><tr><td>Symmetric</td><td>0.7692</td><td>0.001</td></tr></tbody></table>

<table><caption>**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). Significance values obtained after 10,000 permutations of the original data. The bold show the highest correlation.</caption><thead><tr><th>Source of variation</th><th>Iteration/module</th><th>Correlation</th><th>p value</th></tr></thead><tbody><tr><td rowspan="3">Total sample</td><td>Ind x FA</td><td>0.4952</td><td>&lt;0.0001</td></tr><tr><td>Ind x error</td><td><strong>0.5288</strong></td><td>&lt;0.0001</td></tr><tr><td>FA x error</td><td>0.3829</td><td>&lt;0.0001</td></tr></tbody></table>

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><caption>**Table 3.** Values of total variance, variance of the eigenvectors, eigenvectors variance scaled by the total variance (Evstv), and previously scaled also for the total variance and number of variables (Evstvnv). Data for the total sample are presented, see Table S2 for details on subsamples. Italic indicate morphological integration in the symmetric component, with greater values in bold and underline. The same criteria for the asymmetric component, in bold tones.</caption><thead><tr><th rowspan="2">Source of variation</th><th rowspan="2">Hypothesis</th><th rowspan="2">Module</th><th colspan="2">Total Variance</th><th colspan="2">Eigenvalue Variance</th><th colspan="2">Evstv*</th><th colspan="2">Evstvnv**</th></tr><tr><th>Asymmetric</th><th>Symmetric</th><th>Asymmetric</th><th>Symmetric</th><th>Asymmetric</th><th>Symmetric</th><th>Asymmetric</th><th>Symmetric</th></tr></thead><tbody><tr><td rowspan="12">Total sample***</td><td rowspan="5">Functional</td><td>Total face</td><td>0.0004</td><td>0.0031</td><td>2.20E-10</td><td>1.37E-08</td><td>0.0015</td><td>0.0015</td><td>0.0658</td><td>0.0760</td></tr><tr><td>Mouth</td><td>0.0011</td><td>0.0086</td><td>2.47E-08</td><td>1.97E-06</td><td><strong>0.0191</strong></td><td><strong>0.0264</strong></td><td>0.1372</td><td><strong>0.2416</strong></td></tr><tr><td>Ears</td><td>0.0003</td><td>0.0013</td><td>7.14E-10</td><td>1.43E-08</td><td>0.0077</td><td>0.0083</td><td>0.1003</td><td>0.0998</td></tr><tr><td>Eyes</td><td>0.0004</td><td>0.0013</td><td>1.58E-09</td><td>1.38E-08</td><td>0.0118</td><td>0.0087</td><td><strong>0.1195</strong></td><td>0.0791</td></tr><tr><td>Nose</td><td>0.0005</td><td>0.0051</td><td>2.98E-09</td><td>3.33E-07</td><td>0.0141</td><td>0.0127</td><td>0.0633</td><td>0.0792</td></tr><tr><td rowspan="4">Middle line</td><td>Sagital</td><td>0.0007</td><td>0.0062</td><td>2.05E-09</td><td>1.43E-07</td><td>0.0041</td><td>0.0037</td><td>0.0539</td><td><strong>0.0667</strong></td></tr><tr><td>Midsagital</td><td>0.0004</td><td>0.0013</td><td>1.58E-09</td><td>1.38E-08</td><td>0.0118</td><td>0.0087</td><td>0.1195</td><td>0.0791</td></tr><tr><td>Lateral</td><td>0.0003</td><td>0.0013</td><td>7.14E-10</td><td>1.43E-08</td><td>0.0077</td><td>0.0083</td><td>0.1003</td><td>0.0998</td></tr><tr><td>Middle</td><td>0.0005</td><td>0.0039</td><td>1.84E-09</td><td>1.05E-07</td><td>0.0063</td><td>0.0071</td><td>0.1018</td><td>0.1280</td></tr><tr><td rowspan="2">Facial third</td><td>Inferior</td><td>0.0009</td><td>0.0090</td><td>8.92E-09</td><td>1.33E-06</td><td>0.0123</td><td>0.0164</td><td>0.1006</td><td>0.1819</td></tr><tr><td rowspan="2">Neuro-Splach</td><td>Neurocranium</td><td>0.0003</td><td>0.0020</td><td>4.91E-10</td><td>1.98E-08</td><td>0.0042</td><td>0.0052</td><td>0.0709</td><td>0.0879</td></tr><tr><td>Splanchnocranium</td><td>0.0005</td><td>0.0037</td><td>4.67E-10</td><td>3.14E-08</td><td>0.0023</td><td><strong>0.0023</strong></td><td>0.0598</td><td>0.0736</td></tr></tbody></table>

<table><caption>**Table 4.** Results of the regression of morphological integration on the total variance (Evstv) for symmetric and asymmetric shape components. R², RMSE, df, SS, MS, F, and p values are shown. Bolded values indicates significance at α = 0.05.</caption><thead><tr><th>Component</th><th>Effect</th><th>r<sup>2</sup></th><th>RMSE</th><th>df</th><th>SS</th><th>MS</th><th>F</th><th>p</th></tr></thead><tbody><tr><td rowspan="5">Asymmetric</td><td>Total sample</td><td>0.0700</td><td>0.0280</td><td>1</td><td>0.0005</td><td>0.0005</td><td>0.6775</td><td>0.4317</td></tr><tr><td>American</td><td>0.3016</td><td>0.0236</td><td>1</td><td>0.0022</td><td>0.0022</td><td>3.8865</td><td>0.0802</td></tr><tr><td>European</td><td>0.0383</td><td>0.0320</td><td>1</td><td>0.0004</td><td>0.0004</td><td>0.3586</td><td>0.5640</td></tr><tr><td>Admixed/Heterozygous</td><td>0.0894</td><td>0.0331</td><td>1</td><td>0.0010</td><td>0.0010</td><td>0.8839</td><td>0.3717</td></tr><tr><td>Homozygous</td><td>0.0037</td><td>0.0277</td><td>1</td><td>0.0000</td><td>0.0000</td><td>0.0337</td><td>0.8585</td></tr><tr><td rowspan="5">Symmetric</td><td>Total sample</td><td>0.4892</td><td>0.0410</td><td>1</td><td>0.0145</td><td>0.0145</td><td>8.6191</td><td><strong>0.0166</strong></td></tr><tr><td>American</td><td>0.2203</td><td>0.0316</td><td>1</td><td>0.0025</td><td>0.0025</td><td>2.5427</td><td>0.1453</td></tr><tr><td>European</td><td>0.4213</td><td>0.0432</td><td>1</td><td>0.0122</td><td>0.0122</td><td>6.5528</td><td><strong>0.0307</strong></td></tr><tr><td>Admixed/Heterozygous</td><td>0.3579</td><td>0.0351</td><td>1</td><td>0.0062</td><td>0.0062</td><td>5.0163</td><td>0.0519</td></tr><tr><td>Homozygous</td><td>0.3328</td><td>0.0405</td><td>1</td><td>0.0074</td><td>0.0074</td><td>4.4886</td><td>0.0632</td></tr></tbody></table>

## 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⁴⁰. 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³³. 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²⁸,³¹,³².

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⁴¹,⁴². 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⁴³, 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⁴⁴⁻⁴⁷. 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
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Linear regression of the variance of the eigenvectors scaled by the total variance phenotypic variance (trace), and scaled by both, the number of variables and the trace for the total sample and ancestry sub-samples. (a) Asymmetric component, (b) symmetric component. Solid black line represents the fit of the model and the dashed line curve $\alpha$ confidence interval 0.01. The reference for the identification of each module appears in the bottom.</figcaption></figure>

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¹⁹,⁴⁸. Our results are in accordance with previous ones¹⁵ in the sense that a traditional modular division into neurocranium, base and splanchnocranium⁹,⁴⁹ 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¹⁴. In addition, Esteve-Altava et al.⁵⁰, 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¹⁹. However, note that the Esteve-Alta anatomical networks paper¹² 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³⁴ the modular pattern seems to remain unaltered across all of them, reinforcing the idea that phenotypic evolution can be canalized throughout stable patterns of covariation⁵¹,⁵². In this regard, González-José et al.⁵¹ 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⁵¹. 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⁴⁴,⁴⁶,⁵³⁻⁵⁵.

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.⁵¹), even when individuals with a disease-altered facial development are considered. For instance, Martínez-Abadías et al.⁵⁶ 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.⁵⁷. 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)⁵⁸.

Regarding the stasis/disruption of covariation patterns, Klingenberg²¹,²⁸,³¹ 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⁵⁶, 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.⁵⁹ 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³⁹,⁶⁰ for flies, and ref.⁶¹ 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.⁵⁶ 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⁶²,⁶³, 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⁶⁴. 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⁶⁵. 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⁶⁵. The
genetic basis of regulation of the eye is well known and is associated with PAX6⁶⁴, with intervention of the sonic hedgehog (SHH) to the up regulation of PAX2 and down regulation of PAX6. Moreover, fibroblast growth factors (FGF), and transcription factors MITF, CHX10, SOX2, BMP-4, LMAF, SIX3 and PROX1 determine the course of ocular substructures.

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⁶⁶. 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⁶⁴. 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⁶⁷. Thus, the hits in different and non-overlapping genetic systems also correlates with different phenotypes: the ear structures nucleated in the pinna⁶⁶, versus localized nasal traits⁶⁷. In some way, these combined 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⁴⁹.

Moreover, natural selection acting quite independently upon the modules⁵⁶, the complexity of structures⁶⁸, 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⁶⁹, maybe due to the fact that the masticatory apparatus grows even up to 25 years⁶⁵,⁷⁰. 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⁷¹⁻⁷³. 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¹⁵,⁵¹,⁵²,⁵⁶,⁷⁴⁻⁷⁶ 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³³, and similar results are corroborated in mandibles of shrews (for details see refs⁷⁷,⁷⁸). 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⁷⁷,⁷⁸. 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³⁴. 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⁷⁹. 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⁸⁰,⁸¹, thus promoting modularity at the developmental level. While previous analyses have hypothesized that selection makes functional and genetic modularity converge²²,⁴⁴,⁴⁵,⁸²,⁸³, other studies indicate that different morphological structures can perform equivalent functions and, therefore, there can be considerable flexibility in the neutral divergence⁸⁴,⁸⁵ or among-trait connectivity, independent of their function¹⁹,⁴⁸. Whatever the case, our results indicate that the modular landscape of the
<table><caption>**Table 5.** Regression results of fluctuating asymmetry (FA) on morphological integration (Evstv) in both components of the shape, for the total sample and subsamples (Native American, European, Admixed/ Heterozygous, Homozygous). R², root mean square error, correlation, degree of freedom, sum squares, mean squares, F value and p value are presented. The bold represent the unique test with statistical significance with α = 0.05. Italic indicate the higher correlations.</caption><thead><tr><th>Component</th><th>Effect</th><th>r<sup>2</sup></th><th>RMSE</th><th>Correlation</th><th>df</th><th>SS</th><th>MS</th><th>F</th><th>p</th></tr></thead><tbody><tr><td rowspan="6">Asymmetric</td><td>Total sample</td><td>0.0709</td><td>0.0280</td><td>0.2662</td><td>1</td><td>0.0005</td><td>0.0005</td><td>0.6864</td><td>0.4288</td></tr><tr><td>American</td><td>0.2011</td><td>0.0252</td><td>0.4484</td><td>1</td><td>0.0014</td><td>0.0014</td><td>2.2649</td><td>0.1666</td></tr><tr><td>European</td><td>0.0051</td><td>0.0326</td><td>-0.0715</td><td>1</td><td>0.0000</td><td>0.0000</td><td>0.0463</td><td>0.8344</td></tr><tr><td>Admixed/Heterozygous</td><td>0.2757</td><td>0.0295</td><td><strong>0.5251</strong></td><td>1</td><td>0.0030</td><td>0.0030</td><td>3.4257</td><td>0.0972</td></tr><tr><td>Homozygous</td><td>0.0021</td><td>0.0277</td><td>0.0457</td><td>1</td><td>0.0000</td><td>0.0000</td><td>0.0189</td><td>0.8938</td></tr><tr><td>Total sample</td><td>0.5104</td><td>0.0401</td><td>0.7145</td><td>1</td><td>0.0151</td><td>0.0151</td><td>9.3839</td><td><strong>0.0135</strong></td></tr><tr><td rowspan="4">Symmetric</td><td>American</td><td>0.5407</td><td>0.0242</td><td>0.7353</td><td>1</td><td>0.0062</td><td>0.0062</td><td>10.5963</td><td><strong>0.0099</strong></td></tr><tr><td>European</td><td>0.4043</td><td>0.0438</td><td>0.6359</td><td>1</td><td>0.0117</td><td>0.0117</td><td>6.1087</td><td><strong>0.0355</strong></td></tr><tr><td>Admixed/Heterozygous</td><td>0.7377</td><td>0.0224</td><td><strong>0.8589</strong></td><td>1</td><td>0.0127</td><td>0.0127</td><td>25.3056</td><td><strong>0.0007</strong></td></tr><tr><td>Homozygous</td><td>0.3437</td><td>0.0401</td><td>0.5862</td><td>1</td><td>0.0076</td><td>0.0076</td><td>4.7125</td><td><strong>0.0580</strong></td></tr></tbody></table>

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.³⁷). 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³⁴,³⁷,⁸⁶ and conforms the main database of the CANDELA consortium³⁴,³⁷,⁶⁶,⁶⁷,⁸⁶,⁸⁷. 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³⁴ and⁸⁶. 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³⁴,⁸⁶, 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.³⁷). The SNPs were pruned to remove Linkage Disequilibrium (LD) 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⁸⁸,⁸⁹, 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⁸⁶ 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⁵⁹. 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⁵⁹. 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⁹⁰,⁹¹. Moss and Young⁹¹ 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¹⁹.

MLMH (Fig. 1b) is aimed to explore if the farther a module is of the sagittal midline, the more asymmetry will display⁹². 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⁶⁸. 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⁹³⁻⁹⁵. 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⁴¹ 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³⁸. 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³⁸. 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³⁸.

**Covariance matrix similarity: exploring developmental pathways.** Similarity between FA and individuals covariance matrices was tested following Klingenberg and McIntyre³⁹. 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⁹⁶,⁹⁷. Statistical significances were determined through matrix permutation tests, with 10,000 iterations, against the null hypothesis of complete dissimilarity between the covariance matrices concerned³⁹.

**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³⁵, 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³⁶. Hallgrimsson et al.¹⁴ 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³⁹.

**References**

1. Klingenberg, C. P. Morphological Integration and Developmental Modularity. *Ann Rev Ecol Syst* **39**, 115–132 (2008).

2. Simon, H. The architecture of complexity: hierarchic systems. *Proc Amer Phil Soc* **106**, 467–482 (1962).

3. Bolker, J. A. Modularity in Development and Why It Matters to Evo-Devo. *Amer Zool* **40**, 770–776 (2000).

4. Mitteroecker, P. The Developmental Basis of Variational Modularity: Insights from QuantitativeGenetics, Morphometrics, and Developmental Biology. *Evol. Biol.* **36**, 377–385 (2009).

5. Enlow, D. H. *Facial growth*. (Saunders, 1990).

6. Lieberman, D. E., Washington, T. G. & Kenneth, M. Basicranial influence on overall cranial Osbjorn M. Pearson. *J Hum Evol*. https://doi.org/10.1006/jhev.1999.0335 (2000).

7. Rosas, A. & Bastir, M. Thin plate spline analysis of allometry and sexual dimorphism in the human craniofacial complex. *Am. J. Phys. Anthropol.* **117**, 236–45 (2002).

8. Goodrich, J. T. Skull base growth in craniosynostosis. *Childs. Nerv. Syst.* **21**, 871–9 (2005).
9. Bastir, M. & Rosas, A. Mosaic evolution of the basicranium in Homo and its relation to modular development. *Evol Biol* **36**, 57–70 (2008).

10. Marcucio, R. S., Young, N. M., Hu, D. & Hallgrimsson, B. Mechanisms that underlie co-variation of the brain and face. *Genesis* **49**, 177–89 (2011).

11. Martínez-Abadias, N. *et al*. Beyond the closed suture in apert syndrome mouse models: Evidence of primary effects of FGFR2 signaling on facial shape at birth. *Dev. Dyn.* **239**, 3058–3071 (2010).

12. Hallgrimsson, B. *et al*. The brachymorph mouse and the developmental-genetic basis for canalization and morphological integration. *Evol. html_ent glyph = '@amp;' ascii = '@&amp;';/Dev.* **8**, 61–73 (2006).

13. Hallgrímsson, B., Lieberman, D. E., Liu, W., Ford-Hutchinson, A. F. & Jirik, F. R. Epigenetic interactions and the structure of phenotypic variation in the cranium. *Evol. Dev.* **9**, 76–91 (2007).

14. Hallgrímsson, B. *et al*. Deciphering the Palimpsest: Studying the Relationship Between Morphological Integration and Phenotypic Covariation. *Evol Biol* **36**, 355–376 (2009).

15. Martínez-Abadias, N. *et al*. Pervasive genetic integration directs the evolution of human skull shape. *Evolution* (N. Y). **66**, 1010–23 (2012).

16. Martínez-Abadias, N. *et al*. The Developmental Basis of Quantitative Craniofacial Variation in Humans and Mice. *Evol. Biol.* **39**, 554–567 (2012).

17. Young, N. M. *et al*. Embryonic bauplans and the developmental origins of facial diversity and constraint. *Development* **141**, 1059–63 (2014).

18. Hallgrimsson, B., Lieberman, D. E., Young, N. M., Parsons, T. & Wat, S. Evolution of covariance in the mammalian skull. *Tinkering microevolution Dev.* **164**–190 (2007).

19. Esteve-Altava, B., Diogo, R., Smith, C., Boughner, J. C. & Rasskin-Gutman, D. Anatomical networks reveal the musculoskeletal modularity of the human head. *Sci. Rep.* **5**, 8298 (2015).

20. Klingenberg, C. P. & Leamy, L. J. Quantitative genetics of geometric shape in the mouse mandible. *Evolution* **55**, 2342–2352 (2001).

21. Klingenberg, C. P., Mebus, K. & Auffray, J. C. Developmental integration in a complex morphological structure: how distinct are the modules in the mouse mandible? *Evol Dev* **5**, 522–531 (2003).

22. Klingenberg, C. P. *In Phenotypic integration: studying the ecology and evolution of complex phenotypes* (eds Pigliucci, M. and Preston, K.) 213–230 (Oxford University Press, 2004).

23. Wilkins, A. *The Evolution of Developmental Pathways*. (Sinauer Associates, 2002).

24. Klingenberg, C. P. Morphometric integration and modularity in configurations of landmarks: tools for evaluating a priori hypotheses. *Evol Dev* **11**, 405–421 (2009).

25. Enlow, D. H. & Hans, M. G. *Essentials of facial growth*. (Saunders, 1996).

26. Nadeau, J. H. *et al*. Pleiotropy, homeostasis, and functional networks based on assays of cardiovascular traits in genetically randomized populations. *Genome Res.* **13**, 2082–91 (2003).

27. Wagner, G. P. & Altenberg, L. *Perspectives: Complex Adaptations and the Evolution of Evolvability*. *Evolution* (N. Y). **50**, 967–976 (1996).

28. Klingenberg, C. In *Variation: A Central Concept in Biology* 219–247 (Elsevier, 2005).

29. Lande, R. Quantitative Genetic Analysis of Multivariate Evolution, Applied to Brain: Body Size Allometry. *Evolution* (N. Y). **37**, 402–416 (1979).

30. Alonso, C. In *Evolving Pathways: Key Themes in Evolutionary Developmental Biology* (eds Minelli, A. & Fusco, G.) 80–99 (Cambridge Univ Press, 2008).

31. Klingenberg, C. & Polak, M. In *Developmental instability: causes and consequences* 427–442 (Oxford University Press, 2003).

32. Klingenberg, C. In *The Biology of Genetic Dominance* (ed. Veitia, R.) 37–51 (LandesBioscience, 2004).

33. Goswami, A., Binder, W. J., Meachen, J. & O'Keefe, F. R. The fossil record of phenotypic integration and modularity: A deep-time perspective on developmental and evolutionary dynamics. *Proc. Natl. Acad. Sci.* **112**, 4891–4896 (2015).

34. Quinto-Sánchez, M. *et al*. Facial asymmetry and genetic ancestry in Latin American admixed populations. *Am. J. Phys. Anthropol.* **157**, 58–70 (2015).

35. Pavlicev, M., Cheverud, J. M. & Wagner, G. P. Measuring Morphological Integration Using Eigenvalue Variance. *Evol. Biol.* **36**, 157–170 (2009).

36. Wagner, G. P. A comparative study of morphological integration in Apis mellifera (Insecta, Hymenoptera). *J. Zool. Syst. Evol. Res.* **28**, 48–61 (1990).

37. Ruiz-Linares, A. *et al*. Admixture in Latin America: geographic structure, phenotypic diversity and self-perception of ancestry based on 7,342 individuals. *PLoS Genet.* **10**, e1004572 (2014).

38. Adams, D. C. Evaluating modularity in morphometric data: Challenges with the RV coefficient and a new test measure. *Methods Ecol. Evol.* **7**, 565–572 (2016).

39. Klingenberg, C. P. & McIntyre, G. S. Geometric morphometrics of developmental instability: analyzing patterns of fluctuating asymmetry with Procrustes methods. *Evolution* (N. Y). **52**, 1363–1375 (1998).

40. Murren, C. J. The integrated phenotype. *Integr. Comp. Biol.* **52**, 64–76 (2012).

41. Bastir, M. & Rosas, A. Hierarchical nature of morphological integration and modularity in the human posterior face. *Am J Phys Anthr.* **128**, 26–34 (2005).

42. Bastir, M. A systems-model for the morphological analysis of integration and modularity in human craniofacial evolution. *J. Anthropol. Sci.* **86**, 37–58 (2008).

43. Bastir, M., Rosas, A. & O'higgins, P. Craniofacial levels and the morphological maturation of the human skull. *J. Anat.* **209**, 637–54 (2006).

44. Cheverud, J. Developmental Integration and the Evolution of Pleiotropy. *Integr. Comp. Biol.* **36**, 44–50 (1996).

45. Wagner, G. P. & Altenberg, L. Complex Adaptations and the Evolution of Evolvability. *Evolution* (N. Y). **50**, 967–976 (1996).

46. Cheverud, J. M. Quantitative genetics and developmental constraints on evolution by selection. *J theor Biol* **110**, 155–171 (1984).

47. Wagner, G. Homologues, Natural Kinds and the Evolution of Modularity. *Integr. Comp. Biol.* **36**, 36–43 (1996).

48. Esteve-Altava, B. & Rasskin-Gutman, D. Beyond the functional matrix hypothesis: a network null model of human skull growth for the formation of bone articulations. *J. Anat.* **225**, 306–16 (2014).

49. Lieberman, D. *The evolution of de human head*. (The Belknap Press of Harvard University Press, 2011).

50. Esteve-Altava, B., Marugán-Lobón, J., Botella, H., Bastir, M. & Rasskin-Gutman, D. Grist for Riedl's mill: a network model perspective on the integration and modularity of the human skull. *J. Exp. Zool. B. Mol. Dev. Evol.* **320**, 489–500 (2013).

51. González-José, R., Van Der Molen, S., González-Pérez, E. & Hernández, M. Patterns of phenotypic covariation and correlation in modern humans as viewed from morphological integration. *Am. J. Phys. Anthropol.* **123**, 69–77 (2004).

52. Marroig, G. & Cheverud, J. M. A comparison of phenotypic variation and covariation patterns and the role of phylogeny, ecology, and ontogeny during cranial evolution of new world monkeys. *Evolution* **55**, 2576–600 (2001).

53. Olson, E. C. & Miller, R. L. *Morphological integration*. (University of Chicago Press, 1958).

54. Cheverud, J. M. Relationships among ontogenetic, static, and evolutionary allometry. *Am J Phys Anthr.* **59**, 139–149 (1982).

55. Cheverud, J. Morphological integration in the saddle-back tamarin (Saguinus fuscollis) cranium. *Am Nat* **145**, 63–89 (1995).

56. Debat, V., Alibert, P., David, P., Paradis, E. & Auffray, J. C. Independence between developmental stability and canalization in the skull of the house mouse. *Proc. Biol. Sci.* **267**, 423–30 (2000).
57. Martínez-Abadias, N. et al. FGF/FGFR Signaling Coordinates Skull Development by Modulating Magnitude of Morphological Integration: Evidence from Apert Syndrome Mouse Models. *PLoS One* **6**, e26425 (2011).

58. Martínez-Abadias, N. et al. Tissue-specific responses to aberrant FGF signaling in complex head phenotypes. *Dev. Dyn.* **242**, 80–94 (2013).

59. Klingenberg, C. P., Barluenga, M. & Meyer, A. Shape analysis of symmetric structures: quantifying variation among individuals and asymmetry. *Evolution* (N. Y.) **56**, 1909–1920 (2002).

60. Klingenberg, C. P. & Zaklan, S. D. Morphological integration between development compartments in the Drosophila wing. *Evolution* (N. Y.) **54**, 1273–1285 (2000).

61. Klingenberg, C. P., Badyaev, A. V., Sowry, S. M. & Beckwith, N. J. Inferring Developmental Modularity from Morphological Integration: Analysis of Individual Variation and Asymmetry in Bumblebee Wings. *Am Nat* **157** (2001).

62. Grammer, K. & Thornhill, R. Human (Homo sapiens) facial attractiveness and sexual selection: the role of symmetry and averageness. *J. Comp. Psychol.* **108**, 233–42 (1994).

63. Perrett, D. I. et al. Symmetry and Human Facial Attractiveness. *Evol. Hum. Behav.* **20**, 295–307 (1999).

64. Sadler, T. *Medical embryology*. (Wolters Kluwer, Lippincott Williams & Wilkins, 2012).

65. Sperber, G. H. *Craniofacial development*. (BC Decker Inc, 2001).

66. Adhikari, K. et al. A genome-wide association study identifies multiple loci for variation in human ear morphology. *Nat. Commun.* **6**, 7500 (2015).

67. Adhikari, K. et al. A genome-wide association scan in admixed Latin Americans identifies loci influencing facial and scalp hair features. *Nat. Commun.* **7**, 10815 (2016).

68. Aparicio, J. M. & Bonal, R. Why do some traits show higher fluctuating asymmetry than others? A test of hypotheses with tail feathers of birds. *Heredity* (Edinb). **89**, 139–44 (2002).

69. Paschetta, C. The influence of masticatory loading on craniofacial morphology: A test case across technological transitions in the Ohio valley. *Am J Phys Anthr.* **141**, 297–314 (2010).

70. Sperber, G. H. New Insights in Facial Development. *Semin. Orthod.* **12**, 4–10 (2006).

71. Hubbe, M., Hanihara, T. & Harvati, K. Climate signatures in the morphological differentiation of worldwide modern human populations. *Anat Rec A 292*, 1720–1733 (2009).

72. Noback, M. L., Harvati, K. & Spoor, F. Climate-related variation of the human nasal cavity. *Am J Phys Anthr.* **145**, 599–614 (2011).

73. de Azevedo, S. et al. Nasal airflow simulations suggest convergent adaptation in Neanderthals and modern humans. *Proc. Natl. Acad. Sci. USA* **114**, 12442–12447 (2017).

74. Strait, D. S. Integration, phylogeny, and the hominid cranial base. *Am J Phys Anthr.* **114**, 273–297 (2001).

75. Ackermann, R. R. Ontogenetic integration of the hominoid face. *J Hum Evol* **48**, 175–197 (2005).

76. DeLeon, V. B. & Richtsmeier, J. T. Fluctuating asymmetry and developmental instability in sagittal craniosynostosis. *Cleft Palate. Craniofac. J.* **46**, 187–96 (2009).

77. Badyaev, A. V., Foresman, K. R. & Fernandes, M. V. Stress and development stability: vegetation removal causes increased fluctuating asymmetry in shrews. *Ecology* **81**, 336–345 (2000).

78. Badyaev, A. V. & Foresman, K. R. Evolution of morphological integration. I. Functional units channel stress-induced variation in shrew mandibles. *Am. Nat.* **163**, 868–79 (2004).

79. Breuker, C. J., Debat, V. & Klingenberg, C. P. Functional evo-devo. *Trends Ecol. Evol.* **21**, 488–92 (2006).

80. Enlow, D. H. In *Essentials of Facial Growth*. (eds Enlow, D. H. & Hans, M.) 57–78 (Saunders, 1996).

81. West-Eberhard, M. Developmental plasticity and evolution. 794 (2003).

82. Leamy, L. Is developmental stability heritable? *J. Evol. Biol.* **10**, 21 (1997).

83. Mezey, J. G., Cheverud, J. M. & Wagner, G. P. Is the Genotype-Phenotype Map Modular?: A Statistical Approach Using Mouse Quantitative Trait Loci Data. *Genetics* **156**, 305–311 (2000).

84. Wainwright, P. C., Alfaro, M. E., Bolnick, D. I. & Hulsey, C. D. Many-to-One Mapping of Form to Function: A General Principle in Organismal Design? *Integr. Comp. Biol.* **45**, 256–62 (2005).

85. Young, R. L., Haselkorn, T. S. & Badyaev, A. V. Functional equivalence of morphologies enables morphological and ecological diversity. *Evolution* **61**, 2480–92 (2007).

86. Quinto-Sánchez, M. et al. Socioeconomic Status Is Not Related with Facial Fluctuating Asymmetry: Evidence from Latin-American Populations. *PLoS One* **12**, e0169287 (2017).

87. Adhikari, K. et al. A genome-wide association scan implicates DCHS2, RUNX2, GLI3, PAX1 and EDAR in human facial variation. *Nat. Commun.* **7**, 11616 (2016).

88. Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. *Am. J. Hum. Genet.* **81**, 559–75 (2007).

89. Yang, J., Lee, S. H., Goddard, M. E. & Visscher, P. M. GCTA: a tool for genome-wide complex trait analysis. *Am. J. Hum. Genet.* **88**, 76–82 (2011).

90. Diogo, R. & Wood, B. The broader evolutionary lessons to be learned from a comparative and phylogenetic analysis of primate muscle morphology. *Biol. Rev. Camb. Philos. Soc.* **88**, 988–1001 (2013).

91. Moss, M. & Young, R. A functional approach to craniology. *Am. J. Phys. Anthropol.* **18**, 281–92 (1960).

92. Özener, B. Does urban poverty increase body fluctuating asymmetry? *Coll. Antropol.* **35**, 1001–5 (2011).

93. Farkas, L. G., Katic, M. J., Hreczko, T. A., Deutsch, C. & Munro, I. R. Anthropometric proportions in the upper lip-lower lip-chin area of the lower face in young white adults. *Am. J. Orthod.* **86**, 52–60 (1984).

94. Larrabee, W., Makielski, K. & Henderson, J. In *Surgical anatomy of the face* (eds Larrabee W. F., Makielski, K. & Henderson, J.) 22–28 (Lippincott Williams & Wilkins, 2004).

95. Anic-Milosevic, S., Mestrovic, S., Prlić, A. & Slaj, M. Proportions in the upper lip-lower lip-chin area of the lower face as determined by photogrammetric method. *J. Craniomaxillofac. Surg.* **38**, 90–5 (2010).

96. Klingenberg, C. P., Debat, V. & Roff, D. A. Quantitative genetics of shape in cricket wings: developmental integration in a functional structure. *Evolution* (N. Y.) **64**, 2935–2951 (2010).

97. Muñoz-Muñoz, F. et al. Drosophila wing modularity revisited through a quantitative genetic approach. *Evolution* (N. Y.) **70**, 1530–1541 (2016).

## Acknowledgements

We are grateful to the volunteers of the CANDELA sampling for their enthusiastic support for this research. We thank Alvaro Alvarado, Mónica Ballesteros Romero, Ricardo Cebrecos, Miguel Ángel Contreras Sieck, Joyce De la Piedra, María Teresa Del Solar, William Flores, Martha Granados Riveros, Ilich Jafet Moreno, Jodie Lampert, Paola León-Mimila, Francisco Quispealaya, Diana Rogel Diaz, Ruth Rojas, Norman Russell, Vanessa Sarabia, Rosilene Paim, Ricardo Gunski, Sergeant João Felisberto Menezes Cavalheiro and Major Eugenio Correa de Souza Junior for assistance with volunteer recruitment, sample processing and data entry. The following institutions kindly provided facilities for the assessment of volunteers: Escuela Nacional de Antropología e Historia and Universidad Nacional Autónoma de México (México); Pontificia Universidad Católica del Perú, Universidad de Lima and
Universidad Nacional Mayor de San Marcos (Perú); Universidade Federal do Rio Grande do Sul (Brazil). This work was funded by grants from: the Leverhulme Trust (F/07 134/DF to ARL), BBSRC (BB/1021213/1 to ARL, and the CONICET Latin American grant program (No. 41488 to M.Q-S). Raw data regarding symmetric and asymmetric coordinates, along with associated metadata is available at: laofunam.com/data

**Author Contributions**

M.Q.S., R.G.J. designed the study, M.Q.S., C.C., P.N., C.C.S.C., C.P., V.R., S.A., V.A.A., K.A., L.C., J.G.V., P.E., F.A., T.H., C.J., W.A., M.F.G., C.G., G.P., L.S.F., M.C.B., S.C.Q., F.R., G.B., J.R., A.R.L., R.G.J., F.M.M. data acquisition-curation, M.Q.S., F.M.M., R.G.J., J.G.V., C.C., P.N., C.C.S.C., C.P., S.A., V.R., V.A.A., K.A., M.F.G., T.H., P.E., F.A., C.J., W.A., C.G., G.P., G.B., M.C.B., S.C.Q., F.R., J.R., A.R.L. formal analysis, A.R.L., R.G.J. funding acquisition, M.Q.S., F.M.M., R.G.J., C.C. methodology, C.C. software, M.Q.S. visualization, photography and design, R.G.J., M.Q.S. writing original draft.

**Additional Information**

**Supplementary information** accompanies this paper at https://doi.org/10.1038/s41598-018-19324-y.

**Competing Interests:** The authors declare that they have no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_6.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018