Morphological consequences of artificial cranial deformation: Modularity and integration.

The cranium is an anatomically complex structure. One source of its complexity is due to its modular organization. Cranial modules are distinct and partially independent units that interact substantially during ontogeny thus generating morphological integration. Artificial Cranial Deformation (ACD) occurs when the human skull is intentionally deformed, through the use of different deforming devices applied to the head while it is developing. Hence, ACD provides an interesting example to assess the degree to which biomechanical perturbations of the developing neurocranium impact on the degree of morphological integration in the skull as a whole. The main objective of this study was to assess how ACD affects the morphological integration of the skull. This was accomplished by comparing a sample of non-deformed crania and two sets of deformed crania (i.e. antero-posterior and oblique). Both developmental and static modularity and integration were assessed through Generalized Procrustes Analysis by considering the symmetric and asymmetric components of variation in adults, using 3D landmark coordinates as raw data. The presence of two developmental modules (i.e. viscero and neurocranium) in the skull was tested. Then, in order to understand how ACD affects morphological integration, the covariation pattern between the neuro and viscerocranium was examined in antero-posterior, oblique and non-deformed cranial categories using Partial Least-Squares. The main objective of this study was to assess how ACD affects the morphological integration of the skull. This was accomplished by comparing a sample of deformed (i.e. antero-posterior and oblique) and non-deformed crania. Hence, differences in integration patterns were compared between groups. The obtained results support the modular organization of the human skull in the two analyzed modules. The integration analyses show that the oblique ACD style differentially affects the static morphological integration of the skull by increasing the covariance between neuro and viscerocranium in a more constrained way than in antero-posterior and non-deformed skulls. In addition, the antero-posterior ACD style seems to affect the developmental integration of the skull by directing the covariation pattern in a more defined manner as compared to the other cranial categories.

Introduction

There is longstanding scientific tradition that has analyzed the complex nature of human skull using developmental [1], functional [2] and evolutionary [3] approaches, among others. The diverse components of the skulls are integrated with respect to each other because they have developed, functioned and evolved jointly [4]. The morphological integration of the skull is an inexorable process because cranial components are packed together, they have a common developmental background, functional demands exerting reiterative pressures and all the cranial portions share an evolutionary history [5,6]. This integration among different components of the cranium can occur through many biological processes such as hormonal influences, pleiotropy, scaling, regional interactions, among others [7]. However it is important to notice that even though integration is a common biological phenomenon manifested at different levels [7-9], the underlying causes generating it are in most situations not directly observable. Although a recently proposed model-bound approach for understanding morphological evolution of the human skull considers not only statistical analysis of form, but also quantitative genetics and explicit evolutionary hypotheses, like neutral theory [10], morphological integration is typically studied by analyzing the covariation among morphological traits [11]. Nonetheless, this integration is not absolute but organized in units that are relatively independent while participating to generate a structure that act as functional whole. This fact has been recognized for decades [12-16]. These independent units are nevertheless integrated internally, and are operationally known as modules [7]. Even though the majority of the studies on modularity and integration have focused on variation among individuals within populations, there are more levels of variation that exhibit modularity and integration [9], deriving from distinct sources such as genetic variation, phenotypic plasticity, fluctuating asymmetry, evolutionary change, among others [10,17].

Morphological integration and modularity has been studied in skulls by defining modules based in its differential embryological origin: the viscerocranium -derived from branchial arches-, and the neurocranium—derived from neural crest and paraxial mesoderm cells [4,18-26]. There is some evidence that the human skull has a modular behavior divided in two main modules based in its differential embryological origin: the viscerocranium or facial skeleton and the neurocranium or braincase [24,25,27-31]. Based on a similar criterion [32], the latter has been also divided in two additional regions: i.e. the basicranium with endochondral ossification from a cartilaginous mesodermic precursor (chondrocranium) [3,33-37], and the calvaria whose bones derive from the desmocranium, which has an intramembranous ossification from paraxial mesoderm and neural crests cells [38]. In fact, different studies have shown that the patterns of covariation and correlation between different parts of the cranium are partially independent, thus suggesting that they behave as partially independent units [15,39-42]. On the other hand, it has been also suggested that the development of the skull is extremely integrated both functionally and ontogenetically, so any mechanic perturbation could affect not only the immediately adjacent tissues but also other structures that covary with them [43]. For example, children with cartilage growth defects in the cranial base often develop abnormal calvarial morphologies, malocclusion and concave faces [43,44]. Likewise, children with craniosynostoses usually develop malformations of the skeletal face and base [45,46].

In the present study artificial cranial deformation (ACD) is used as a proxy for assessing the morphological integration and modularity of the human skull. ACD consists in the modification of the magnitude and direction of the normal vectors of growth and development of the skull, by using compressive forces generated by deforming devices during the early years of post-natal life [47]. It is a cultural practice which modifies the newborn’s natural head shape to achieve a certain desired morphology, being since long considered either as an identity symbol and/or social status marker in different populations around the world [48-50]. Depending on the particular tradition, ACD comprises the application of different methods, such as bandages, boards, stones, pads, or a combination of these to the neonate skull [48]. Ethnographic evidence shows even that intentional head deformation can be even achieved by manual modification [51], and has been carried out not only in pre-historical but in historical [52], and present time societies too [53]. This cranial modification practice was widely distributed both geographically and temporally [1,49,51,54-56]. Despite findings about the presence of ACD in Central Asia [57-59], in Near and Middle East [60-62], in Eastern [63-67], and Western Europe [68-71], ACD exhibits particularly high frequencies in the Americas as exemplified in [72-74], especially in the South-Central Andes [1,51,55,73,75-80]. Regarding the effect that ACD has on the skull the available evidence shows that ACD is not confined to the calvaria but also affects the cranial base. Despite the fact that deforming devices are mostly placed on the neurocranium [47,81-84], as a result of ACD practice cranial base decreases its normal angle and increases the anterior-inferior projection of the viscerocranium. It seems that ACD influences the complete skull by affecting its globularity, which could mean that the cranial vault behaves as a highly constrained trait [85]. Under certain circumstances ACD could even cause the death of the newborn, as is evidenced by the pathognomonic periosteal reaction at the occipital and parietals bones found in the artificially deformed skull of a one year old children from a late Inca burial [86].

Considering that the ACD constrains the skull’s normal ontogenetic growth direction, it would be expected that the expansion of the encephalon will continue generating internal pressures that will be redirected towards those areas of the cranial vault that are not restricted by the application of the deforming device. In fact, it has been shown that ACD typically generates some amount of compensatory growth in the sutures [87,88]. It seems likely that when the deforming device is applied for a long period, the internal pressures within the skull will continue, however they will be more focused to other cranial regions. This can probably increase or decrease the covariation between the different cranial regions, modifying the pattern of shape/size variation of other bones like the zygomatic bone [89], the shape and volume of craniofacial cavities [90,91], and the basicranial, and the relative position of the mandibular fossae [92]. Therefore, ACD provides an interesting example to assess the degree to which biomechanical perturbations of the developing neurocranium impact on the degree of morphological integration in the cranium as a whole. The null hypotheses that lead this research are:

h0a: The human skull does not exhibit a modular structure corresponding to the viscero and neurocranium.

h0b: ACD, independently of the type of deformation, does not change the degree of morphological integration of the skull compared with non-deformed individuals.

Materials and method

The sample comprised 269 archaeological skulls from Northern Chile (Table 1). Individual numbers and complete repository information, including museum name and geographic location are provided in S1 Table. No permits were required for the described study, which complied with all relevant regulations. The skulls were scanned using a NextEngine Desktop 3D Scanner (model 2020i; Santa Monica, CA) in two different positions, each of them composed by ten divisions to allow a complete three-dimensional (3D) reconstruction. The number of divisions controls the degree of rotation between scans and the total number of scans, so each division correspond to 36° of the total rotation. 3D models were generated by trimming the unwanted artifacts, aligning the two different positions and fusing the two surfaces of each skull using both Meshlab v.1.3.3 (http://meshlab.sourceforge.net/) and ScanStudio™ (http://www.nextengine.com/). All the individuals analyzed in this study were adults. The following inclusion criteria were followed: third molar eruption and the fusion of the spheno-occipital synchondrosis. Deformed skulls were classified as antero-posterior or oblique deformations based on a simplified version of the classification proposed by [93] (Fig 1). Antero-posteriorly deformed skulls were characterized by an antero-posterior vault compression, resulting in the flattening of the frontal and occipital bones, alongside a lateral expansion of the head. On the other hand, oblique skulls were defined according to their conical and lengthened cranial vault, with a relatively narrowed medial-lateral dimension.

Table 1

Sample used in the present study.

Origin	Antero-posterior	Non-deformed	Oblique
Caspana	5	5	16
Catarpe 2	4	31	10
Chorrillos	0	0	20
Chunchuri	8	4	5
Coyo Oriente	8	10	2
Coyo 3	3	5	0
Larache	4	11	0
Quitor 5	0	2	1
Quitor 6	4	3	2
Sequitor Alambrado	0	5	1
Sequitor Alambrado Oriente	4	5	0
Solcor	8	3	1
Solcor 3	4	8	5
Solor 3	1	2	3
Tchecar	12	9	1
Tchilimoya	4	11	4
Toconao	0	3	0
Toconao Oriente	2	1	1
Toconce	1	0	2
Topater	2	0	3
	74	118	77

Origin of the analyzed sample (archaeological site), and type of artificial cranial deformation. Skulls are housed at Museo Arqueológico R.P. Gustavo Le Paige, San Pedro de Atacama

Chile, excluding Chorrillos (Corporación Cultural y de Turismo, Calama, Chile), and

Chunchuri (Musée de l’Homme, Paris, France).

Fig 1

Representative 3D models.

3D surfaces of the different cranial categories analyzed in the present paper and their corresponding anatomical views.

This first step was followed by a geometric morphometrics generalized Procrustes superimposition of landmark data. This branch of shape analysis has been usually understood as the quantitative study of shape and its covariates [94]. Landmark acquisition was carried out by one of us (TP) in Landmark Editor v.3.6 software (IDAV) [95] by collecting 31 homologous and well-defined 3D points using the nomenclature of [96] (Table 2; Fig 2). These raw coordinates are provided in S1 Table. Geometric morphometrics and statistical analysis were carried out in MorphoJ v. 1.06d [97] and the R package ‘geomorph’ [98]. One surface file was warped using Landmark Editor v.3.6 [95] in order to visualize shape changes depending on the applied analysis. Measurement error (ME) has a critical importance when performing morphometric analyses, therefore to assess it, a sub-sample of 170 skulls were digitized twice and compared via a Procrustes ANOVA [99].

Table 2

Landmarks used in the present study.

Number	Name	Definition	Module
1	Prostion (pr)	The most anterior point on the alveolar ridge of the maxilla between the central incisors in the median sagittal plane.	Viscerocranium
2	Nasospinale (ns)	The lowest point of the lower border of the piriform aperture, projected into the median sagittal plane.	Viscerocranium
3	Nasion (n)	Crossing point of the frontonasal suture with the median sagittal.	Viscerocranium
4	Glabella (g)	Most anterior point, in the mediansagittal, betwen the superciliary arches	Neurocranium
5	Bregma (b)	The point at which the sagittal suture meets the coronal suture.	Neurocranium
6	Lambda (l)	The point at which the two parts of the lambdoidal suture meet the sagittal suture.	Neurocranium
7	Inion (i)	The point in the median sagittal plane, in which the two superior nuchal lines join.	Neurocranium
8	Opistion (o)	The point at which the posterior border of the foramen magnum is cut by the medial sagittal plane.	Neurocranium
9	Basion (ba)	The point at which the anterior margin of the foramen magnum is cross by the median sagittal plane.	Neurocranium
10	Sphenobasion (sphba)	Intersection of the sphenooccipitalsynchondrosis with the median sagittal plane.	Neurocranium
11,18	Zygomaxillare (rzm)	Lowest point of the right zygomaticomaxillary suture	Viscerocranium
12, 19	Zygoorbitale (rzo)	Upper point of the right zygomaticomaxillary suture.	Viscerocranium
13, 20	Frontomalare orbitale (rfmo)	The only point on the lateral right orbital rim, at which it is cross by the frontozygomatic suture.	Viscerocranium
14, 21	Auriculare (rau)	The point that is perpendicular to the center of the right porus acusticus externus located on the zygomatic root.	Neurocranium
15, 22	Entomion (ren)	The point at which the right squamous suture passes into the right parietomastoidea suture.	Neurocranium
16, 23	Asterion (rast)	The point at which the right lambdoidea, occipitomastoid and parieto-mastoid the sutures meet.	Neurocranium
17, 24	Mastoideale (rms)	The most inferior point of the right mastoid process.	Neurocranium
25, 29	Krotaphion (rk)	The posterior point of the right sphenoparietalis suture	Neurocranium
26, 30	Sphenion (rsphn)	The anterior point of the rigth sphenoparietalis suture	Neurocranium
27, 31	Stephanion (rst)	The point where the coronal suture meets the right temporal line	Neurocranium
28	Obelion (ob)	The point between a transverse line connecting the two parietal foramina and the sagittal suture.	Neurocranium

Number, name, definition, and location of the landmarks used for carrying out the geometric morphometrics workflow.

Fig 2

Landmarks used to carry out the GM analyses.

An oblique artificially deformed skull from the archaeological site of Chorrillos is shown (blue = neurocranium; fuchsia = viscerocranium).

Human skulls exhibit object symmetry (i.e. the axis of symmetry passes through mid-plane, and hence the left and right halves are mirror images of each other), therefore the recommendations of [100,101] were followed. A generalized Procrustes analysis taking into account object symmetry was performed to remove the differences due rotation, scale and translation. Due to the presence of object symmetry, two separate matrices were generated, representing the symmetric and asymmetric components of variation respectively [101]. The symmetric component represents shape variation among individuals in what could be regarded as a left-right average, thus being useful to study static integration and modularity. On the other hand, the asymmetric component represents the differences between the original and mirrored configurations, being used in the present study to analyze developmental integration and modularity. These two components represent the shape variables that were used in the subsequent analysis. Due to the fact that there is a morphological continuum between deformed and non-deformed skulls (with slightly deformed skull in between), in order to test if it was possible to distinguish between deformed (i.e. antero-posterior and oblique) and non-deformed skulls a canonical variate analysis (CVA) of the symmetric component was performed using the Procrustes coordinates as raw data. The posterior probabilities from the cross-validated CVA are available in S3 Table (Linear Discriminant Analysis using the Procrustes coordinates of the first 11 PCs after a broken-stick model for 3 classes /'Antero-Posterior', 'Non.deformed', 'Oblique'/, Accuracy = 0.8215613, Kappa = 0.7222987). A principal component analysis (PCA) was also carried out in order to quantify skull shape variation. In addition, pairwise PERMANOVA tests with a Holm-Bonferroni correction for multiple comparisons were carried out to test for shape differences between the three cranial categories [102] (i.e., non-deformed, oblique and antero-posterior). The PERMANOVA corresponds to a non-parametric test of significant difference between two or more classes based on a distance measure [103]. In the present work Euclidean distances calculated from the PC scores obtained from the PCA were used as similarity index and 10.000 permutations were carried out. These tests were carried out by adapting the adonis() function from the R package ‘vegan’ v. 2.5–6 to perform multiple comparisons.

As previously mentioned, depending on the processes responsible for integration and modularity, several levels of integration can be distinguished (e.g. static, developmental, evolutionary, among others) [9]. In the present study two main integration levels were analyzed: a) static integration, which is merely the level of variation among individuals in a homogeneous sample, where all specimens are from the same species and ontogenetic stage (this level was analyzed using the symmetric component), and b) developmental, which refers to the fact that morphological traits tend to be associated statistically when they share a specific function and/or a synchronic appearance during embryonic development (this level was analyzed using the asymmetric component) [9].

From a morphometric point of view, modularity is manifested as a strong correlation within the components of a module, versus a weak covariation between modules [7]. There are several methods to statistically test modularity hypotheses (see for example [4] for a review), including the vectorial correlation (RV) coefficient [104], which is a multivariate analogue of a correlation, where the landmark configurations are divided into the hypothetical modules. However, this coefficient has been criticized because it can be adversely affected by the attributes of the data under analysis (i.e. sample size and the number of variables) [105]. Therefore, the covariance ratio (CR) test of modularity was used instead. This test corresponds to a ratio of between to within modules covariances and can range between zero and more than one and has been proposed to overcome some of the problems of RV coefficient [105].

In the present study, two hypothetical modules were proposed based on developmental grounds for the whole sample: the a) viscerocranium and the b) neurocranium. The two proposed modules are parts of the cranium, being consequently spatially adjacent and within the same structure. Therefore, the analysis was restricted to the generation of only spatially contiguous configurations [106]. Based on the same argument, the modularity test was performed using a single Procrustes fit, which automatically considers the information about the connection and spatial arrangement of the subsets [106]. This analysis was carried out both for the symmetric and asymmetric components of variation separately to analyze both static and developmental modularity.

The morphological integration of the skull components was analyzed by means of a partial least-squares analysis (PLS) [107]. According to the available information (cf. Introduction section, this work), deforming devices were normally applied on the occipital, frontal, parietals and temporal bones (although this varied depending on the specific deforming tradition). Consequently, the PLS analysis was carried out to study the covariation patterns between neurocranium and viscerocranium shape variables. The PLS method finds the principal components of covariation between the sets of variables by means of decomposing the covariance matrix of the two blocks into two sets of eigenvectors (one for each set of variables) and eigenvalues, using a singular value decomposition in order to generate a matrix of singular values (the square-roots of the eigenvalues) [107], a matrix of eigenvectors for the first set of variables and the transpose of the matrix of eigenvectors for the second set of variables [108]. The PLS was carried out within the complete landmark configuration, which means that the analysis considered all the existing covariation between the two blocks, including the part related to relative positions, orientations and sizes of the blocks [7]. This approach is well-suited for morphological integration studies that consider the covariation between parts in the context of an entire structure, where the relative arrangement of the components could produce a significant contribution to the covariation patterns, such as in the case of the human cranium [7]. This procedure was repeated separately for each one of the analyzed categories (i.e. non-deformed, antero-posterior and oblique deformed crania) and for both the symmetric and asymmetric components of variation in order to analyze both static and developmental integration. To compare the differences in the covariation patterns between deformed (i.e. antero-posterior and oblique) and non-deformed skulls, pairwise permutation tests of the angular differences between the PLS vectors were estimated [109]. In addition, an angular comparison between the PC’s and the PLS axes was performed to test if the observed patterns of morphological variation (i.e. PC’s) followed a similar trajectory as the observed covariation pattern between the neuro and the viscerocranium (i.e. PLS axes) [110]. In a similar fashion, from all the possible comparisons only a few of them were significant (Symmetric component: θ PC1 v/s PLS1 = 13.361°, P-value <0.00001; Asymmetric component: θ PC1 v/s PLS1 = 34.190°, P-value <0.00001; θ PC2 v/s PLS2 = 49.545°, P-value <0.00001).

Finally, in order to quantify the overall similarity of covariance matrices, matrix correlations and the associated permutation tests were computed, using procedures adapted for geometric morphometrics and object symmetry [7,99,101]. These procedures were carried out to analyze whether the observed patterns of variation were similar between the different categories under study (i.e. non-deformed, antero-posterior and oblique crania). This procedure was carried out as well for both the symmetric and the asymmetric component of variation. All the analyses were carried out considering alpha = 0.05.

Results

In geometric morphometrics, hypothesis testing is based on the linear properties of the vectors obtained after applying Procrustes analysis onto the raw data (2 o 3 dimensional landmark coordinates). Later those vectors (i.e. the shape and size components of form) are used as analogs of the traditional interlandmark ("point to point") vectors. As a result, the linearity of the shape vectors allow to apply the toolkit of standard multivariate statistical analyses in order to explore (i.e. PCA), and contrast (i.e. CVA, PERMANOVA, PLS) hypotheses for assessing the pattern of observed shape variation.

Measurement error

The Procustes ANOVA used to measure intra-observer error in the sub-sample showed that the mean square for individual variation exceeded measurement error, so the effect of measurement error was negligible (see S2 Table for further details). Measurement error was also quantified as repeatability using a ratio of the among-individual to the sum of the among-individual and measurement error components as explained in [33]. Shape repeatability was 0.955, which indicates a minimal ~4% error.

Canonical variate analysis

Antero-posteriorly deformed and non-deformed skulls were significantly different (Mahalanobis distance: 3.2198; P-value: <0.0001; 10.000 perm.), as well as antero-posterior and oblique deformed crania (Mahalanobis distance: 3.2177; P-value: <0.0001; 10.000 perm.) and non-deformed and oblique skulls (Mahalanobis distance: 2.8743; P-value: <0.0001; 10.000 perm.). There was only a slight overlap between them (Fig 3) that could perhaps be attributed to either misclassification or due to the presence of skulls showing only slight deformations. The skulls presenting an antero-posterior deformation seem to be the most different ones, as observable in the higher distance between these individuals and the other cranial categories (i.e. antero-posterior, and non-deformed skulls).

Fig 3

CVA of the symmetric component of the complete dataset.

The ellipses represent the 90% confidence intervals.

Principal component analysis

Morphological variation was relatively widespread in the different PCs obtained from the symmetric component. The first two components accounted for 24.39% and 10.17% of the total shape variation, respectively, depicting some group separation along PC1between deformation types and non-deformed skulls (Fig 4). PC1 noticeably separated between the deforming styles. The skulls with an antero-posterior deformation were located at the left of the scatter plot, while the oblique ones were situated at the right. The non-deformed crania were positioned at the center between the two deforming styles, although separated from them by PC2. The shape changes associated with PC1 can be described as an elongation of the cranial vault towards the right of the scatter plot, while the skulls located at the left exhibit more flattened occipital and frontal bones. On the other hand, PC2 can be related to a relative elevation of the occipital landmarks and a retraction of the upper cranial vault coordinates.

Fig 4

PCA of the symmetric component of the complete dataset.

The ellipses represent the 90% confidence intervals.

Pairwise PERMANOVAs

There were significant differences between the shape of the different cranial categories when comparing them using the PERMANOVA (Table 3). Therefore, the analyzed cranial categories are distinguishable in spite of the slight overlap observed in the PCA (Fig 4).

Table 3

Shape differences.

Cranial category	F-test	R2	Adjusted P-value(Holm-Bonferronicorrection)
Antero-posterior vs. oblique	29.77004	0.16652702	0.003
Antero-posterior vs. non-deformed	26.02839	0.12048598	0.003
Oblique vs. non-deformed	19.4387	0.09150263	0.003

Pairwise PERMANOVA tests used to test for shape differences between cranial categories.

Covariance ratio test

The CR tests applied to test the null hypothesis that the human skull does not show a modular behavior between viscero and neurocranium was rejected for both the symmetric (CR: 0.846; P-value: 0.007; 10.000 perm.) and the asymmetric (CR: 0.637; P-value: 0.023; 10.000 perm.) components of variation. This null hypothesis was also rejected when running the analyses separately for each one of the categories under analysis (i.e. antero-posterior, non-deformed and oblique crania) (see S1 Fig for further details).

Partial least-squares analysis

Symmetric component

The analysis identified the characteristics of shape variation that most covary between the two blocks and highlighted their relative contribution to the total amount of covariation between blocks (Tables 4 and 5). The overall integration measured using the r-PLS was 0.851 (P-value: 0.0001; 10.000 perm.), whereas the PLS1 axes explained 73.55% of the total covariance of the sample (P-value: <0.0001; 10.000 perm.) (Fig 5). This percentage of explained total covariance was similar to the values obtained when analyzing each cranial category separately. The anterior-posterior (r-PLS: 0.845; P-value: 0.0001; 10.000 perm.), the oblique deformed skulls (r-PLS: 0.857; P-value: 0.0001; 10.000 perm.) and the non-deformed skulls also showed relatively similar values, although they were slightly higher in non-deformed skulls. (r-PLS: 0.888; P-value: 0.0001; 10.000 perm.). However, the antero-posterior (PLS1: 52.96% of explained covariance; P-value: <0.0001; 10.000 perm.) and the non-deformed skulls (PLS1: 51.95% of explained covariance; P-value: <0.0001; 10.000 perm.) showed a covariation pattern that was more evenly distributed in several PLS axes, as compared to the oblique crania (PLS1: 75.47% of explained covariance; P-value: <0.0001; 10.000 perm.), which showed a higher percentage of covariance explained by PLS1 (Table 6 and Fig 6). In addition, the results for PLS2 are visually presented in S2 Fig.

Table 4

PLS results of the complete dataset.

Type of shape component	PLS	Singular value	P-value (perm.)	% total covariance	Correlation	P-value (perm.)
Symmetric component	PLS1	0.00033482	<0.0001	73.55	0.85774	<0.0001
	PLS2	0.0001568	<0.0001	16.131	0.93347	<0.0001
	PLS3	0.00008937	<0.0001	5.24	0.55733	0.0018
	PLS4	0.00005396	<0.0001	1.91	0.45688	0.0006
	PLS5	0.00003998	<0.0001	1.049	0.38626	0.0027
	PLS6	0.00003208	<0.0001	0.675	0.4065	<0.0001
	PLS7	0.00002996	<0.0001	0.589	0.37694	0.0007
	PLS8	0.00001957	0.0054	0.251	0.34792	0.0018
	PLS9	0.00001893	<0.0001	0.235	0.31456	0.007
	PLS10	0.00001563	0.0002	0.16	0.34212	0.0003
	PLS11	0.00001043	0.2282	0.071	0.28316	0.0108
	PLS12	0.00000843	0.4457	0.047	0.22377	0.1982
	PLS13	0.00000753	0.2289	0.037	0.17127	0.719
	PLS14	0.00000578	0.5397	0.022	0.18892	0.2476
	PLS15	0.00000519	0.2268	0.018	0.20055	0.0414
	PLS16	0.00000437	0.1146	0.013	0.16144	0.1452
	PLS17	0.00000214	0.699	0.003	0.17011	0.0831
Asymmetric component	PLS1	0.00004092	< .0001	76.517	0.85205	< .0001
	PLS2	0.0000178	< .0001	14.478	0.64862	< .0001
	PLS3	0.00000822	0.0079	3.089	0.51247	< .0001
	PLS4	0.00000745	0.0001	2.537	0.28785	0.0584
	PLS5	0.0000058	0.0007	1.537	0.30112	0.0065
	PLS6	0.00000392	0.0667	0.703	0.23064	0.2107
	PLS7	0.00000316	0.0671	0.456	0.21981	0.205
	PLS8	0.00000237	0.1713	0.257	0.24704	0.0407
	PLS9	0.00000219	0.0137	0.218	0.21242	0.1552
	PLS10	0.00000144	0.3728	0.094	0.2071	0.0821
	PLS11	0.00000116	0.3329	0.061	0.1865	0.085
	PLS12	0.00000083	0.4822	0.032	0.10635	0.8422
	PLS13	0.00000069	0.0815	0.022	0.13983	0.2208

Singular values and pairwise correlations of the PLS scores between corresponding blocks (i.e. neurocranium and viscerocranium) for the complete dataset (10.000 perm.).

Table 5

Singular values and pairwise correlations of the PLS scores between blocks for the different cranial categories.

Cranial category	PLS	Singular value	P-value (perm.)	% total covariance	Correlation	P-value (perm.)
Antero-posterior	PLS1	0.00022702	<0.0001	52.957	0.84244	<0.0001
	PLS2	0.00014615	<0.0001	21.948	0.80669	<0.0001
	PLS3	0.00009995	<0.0001	10.266	0.59511	0.2154
	PLS4	0.00007074	0.005	5.142	0.50818	0.4959
	PLS5	0.00005088	0.1631	2.66	0.59985	0.0259
	PLS6	0.00004936	0.0045	2.504	0.55449	0.0557
	PLS7	0.00003657	0.1146	1.374	0.44052	0.5675
	PLS8	0.00003181	0.0728	1.039	0.51616	0.0648
	PLS9	0.00002688	0.0906	0.742	0.47659	0.147
	PLS10	0.00002189	0.2138	0.492	0.42679	0.3276
	PLS11	0.00001675	0.6639	0.288	0.42802	0.202
	PLS12	0.00001426	0.6423	0.209	0.35027	0.6087
	PLS13	0.00001172	0.713	0.141	0.3684	0.3146
	PLS14	0.00001034	0.4734	0.11	0.30987	0.5///849
	PLS15	0.00000846	0.4329	0.074	0.28601	0.5716
	PLS16	0.00000594	0.6884	0.036	0.31792	0.1609
	PLS17	0.00000413	0.6133	0.018	0.20155	0.8294
Non-deformed	PLS1	0.00015927	<0.0001	51.948	0.91448	<0.0001
	PLS2	0.00010278	<0.0001	21.632	0.84507	<0.0001
	PLS3	0.00007465	<0.0001	11.413	0.621	0.0335
	PLS4	0.00004542	0.0205	4.225	0.51909	0.0647
	PLS5	0.00003983	0.0042	3.249	0.52602	0.0098
	PLS6	0.00003433	0.0016	2.413	0.43808	0.1617
	PLS7	0.00003128	0.0001	2.003	0.47722	0.0111
	PLS8	0.00002324	0.0168	1.106	0.38172	0.31
	PLS9	0.00001893	0.0671	0.734	0.35914	0.3731
	PLS10	0.000014	0.5889	0.402	0.37862	0.1328
	PLS11	0.00001231	0.4508	0.311	0.32843	0.3824
	PLS12	0.00001062	0.3832	0.231	0.33221	0.2094
	PLS13	0.0000077	0.8936	0.122	0.3121	0.2286
	PLS14	0.00000725	0.5126	0.107	0.2775	0.3353
	PLS15	0.00000498	0.9274	0.051	0.17241	0.9729
	PLS16	0.00000431	0.6768	0.038	0.18569	0.7411
	PLS17	0.0000028	0.5845	0.016	0.25283	0.1212
Oblique	PLS1	0.00036424	<0.0001	75.467	0.87452	<0.0001
	PLS2	0.00014999	<0.0001	12.797	0.92793	<0.0001
	PLS3	0.00007936	0.0258	3.582	0.67348	0.0249
	PLS4	0.0000714	0.0009	2.9	0.57496	0.1397
	PLS5	0.00005595	0.0068	1.781	0.54332	0.1515
	PLS6	0.00004043	0.2076	0.93	0.55662	0.0409
	PLS7	0.00003977	0.0046	0.9	0.46563	0.3604
	PLS8	0.00003243	0.0112	0.598	0.50844	0.0653
	PLS9	0.00002395	0.2544	0.326	0.35908	0.8733
	PLS10	0.00002257	0.0313	0.29	0.49321	0.0333
	PLS11	0.00001596	0.5815	0.145	0.43106	0.1573
	PLS12	0.0000148	0.195	0.125	0.53239	0.0016
	PLS13	0.00001192	0.3287	0.081	0.38458	0.2015
	PLS14	0.00000834	0.8654	0.04	0.37618	0.1294
	PLS15	0.00000618	0.9523	0.022	0.28098	0.6061
	PLS16	0.00000473	0.9139	0.013	0.26312	0.4833
	PLS17	0.00000305	0.7927	0.005	0.1995	0.8789

Symmetric component.

Fig 5

PLS1 correlation plot of the complete dataset (symmetric component) for the neurocranium and viscerocranium.

The ellipses represent the 90% confidence intervals.

Table 6

Singular values and pairwise correlations of the PLS scores between blocks for the different cranial categories: asymmetric component.

Cranial category	PLS	Singular value	P-value (perm.)	% total covariance	Correlation	P-value (perm.)
Antero-posterior	PLS1	0.00006503	<0.0001	74.488	0.89509	<0.0001
	PLS2	0.00003158	<0.0001	17.57	0.69941	0.0005
	PLS3	0.00001375	0.3215	3.332	0.47458	0.2284
	PLS4	0.00001025	0.45	1.85	0.48418	0.0982
	PLS5	0.00000868	0.2264	1.327	0.53345	0.0086
	PLS6	0.00000586	0.7514	0.606	0.42417	0.184
	PLS7	0.00000414	0.9481	0.303	0.35894	0.4781
	PLS8	0.00000377	0.6266	0.25	0.31224	0.7608
	PLS9	0.00000258	0.948	0.117	0.40304	0.1411
	PLS10	0.00000199	0.9701	0.07	0.2755	0.803
	PLS11	0.00000157	0.9498	0.044	0.21466	0.9477
	PLS12	0.00000145	0.5391	0.037	0.18816	0.9472
	PLS13	0.00000062	0.9791	0.007	0.24791	0.3926
Non-deformed	PLS1	0.00003401	<0.0001	70.459	0.83665	<0.0001
	PLS2	0.00001627	<0.0001	16.118	0.68219	<0.0001
	PLS3	0.00000932	0.1208	5.287	0.43128	0.1304
	PLS4	0.00000723	0.2256	3.18	0.42484	0.0507
	PLS5	0.00000574	0.3121	2.005	0.35483	0.2513
	PLS6	0.00000421	0.6866	1.08	0.32847	0.2929
	PLS7	0.00000375	0.3475	0.858	0.20194	0.989
	PLS8	0.00000251	0.8762	0.383	0.3035	0.2943
	PLS9	0.00000205	0.8267	0.257	0.3019	0.2034
	PLS10	0.00000173	0.6797	0.183	0.23182	0.6427
	PLS11	0.00000129	0.7887	0.102	0.19566	0.7694
	PLS12	0.00000105	0.5824	0.067	0.15091	0.9182
	PLS13	0.00000059	0.8039	0.021	0.12506	0.9397
Oblique	PLS1	0.00003953	<0.0001	67.032	0.85317	<0.0001
	PLS2	0.0000171	0.1233	12.547	0.58171	0.0527
	PLS3	0.00001355	0.1152	7.871	0.52199	0.0853
	PLS4	0.00001037	0.2228	4.615	0.49875	0.0577
	PLS5	0.00000889	0.0825	3.391	0.36861	0.6503
	PLS6	0.00000649	0.3315	1.807	0.35462	0.6034
	PLS7	0.00000511	0.4062	1.121	0.31515	0.7753
	PLS8	0.00000426	0.2876	0.779	0.30781	0.7411
	PLS9	0.00000266	0.9539	0.303	0.34745	0.3425
	PLS10	0.00000244	0.6727	0.256	0.2689	0.7634
	PLS11	0.00000199	0.5659	0.17	0.27819	0.4924
	PLS12	0.00000133	0.7774	0.076	0.13281	0.9981
	PLS13	0.00000087	0.6378	0.032	0.18318	0.8439

Asymmetric component

Fig 6

PLS1 correlation plots for the different cranial categories here analyzed (symmetric component.

a) Antero-posterior; b) Non-deformed and; c) Oblique.

The overall integration of shape variation component measured using the r-PLS was 0.826 (P-value: 0.001; 10.000 perm.), whilst the PLS1 axes explained 76.52% of the total covariance of the sample (P-value: <0.001; 10.000 perm.). This value was relatively similar to the ones obtained for the antero-posterior (r-PLS: 0.879; P-value: 0.001; 10.000 perm.) and oblique (r-PLS: 0.868; P-value:0.001; 10.000 perm.) crania, while the non-deformed skulls showed a lower integration value (r-PLS: 0.735; P-value:0.001; 10.000 perm.). Both oblique (PLS1: 67.03% of explained covariance; P-value: <0.001; 10.000 perm.) and non-deformed skulls (PLS1: 70.46% of explained covariance; P-value: <0.001; 10.000 perm.) exhibited a similar distribution of explained covariance by their PLS axes. However, the antero-posterior deformed crania a higher level of integration in their first PLS axes (PLS1: 74.49% of explained covariance; P-value: <0.001; 10.000 perm.) (Table 6).

Angular comparison between PLS axes

From all the possible comparisons between the paired singular axes between non-deformed skulls and the antero-posterior and oblique crania, only a few of them were significant. Below are presented the comparisons that were significant for the two PLS axes that accounted for the majority of the total covariance for the symmetric component.

Antero-posterior PLS2 v/s non-deformed PLS2 (θ = 62.762°; P-value = 0.00094)

Antero-posterior PLS1 v/s oblique PLS1 (θ = 61.069°; P-value = 0.00043)

Antero-posterior PLS1 v/s oblique PLS2 (θ = 64.017°; P-value = 0.00164)

Non-deformed PLS1 v/s oblique PLS2 (θ = 14.579°; P-value <0.00001)

Non-deformed PLS2 v/s oblique PLS1 (θ = 41.117°; P-value <0.00001)

Asymmetric component

From all the possible comparisons between the paired singular axes between non-deformed skulls and the antero-posterior deformed and oblique crania, just few of them were significant. Below are presented the comparisons that were significant for the two PLS axes that accounted for the majority of the total covariance for the asymmetric component.

Antero-posterior PLS1 v/s non-deformed PLS1 (θ = 39.064°; P-value <0.00001)

Antero-posterior PLS2 v/s non-deformed PLS2 (θ = 42.940°; P-value <0.00001)

Antero-posterior PLS1 v/s oblique PLS1 (θ = 28.325°; P-value <0.00001)

Non-deformed PLS1 v/s oblique PLS1 (θ = 30.414°; P-value <0.00001)

Angular comparison between principal components and PLS axes

The PC’s and PLS axes of the complete dataset were compared using the same procedure outlined above in order to test whether the observed pattern of variation (PC’s) followed a similar trend as the observed patterns of covariation between viscero- and neurocranium (PLS axes). In a similar fashion, from all the possible comparisons only a few of them were significant (Symmetric component: θ PC1 v/s PLS1 = 13.361°, P-value <0.00001; Asymmetric component: θ PC1 v/s PLS1 = 34.190°, P-value <0.00001; θ PC2 v/s PLS2 = 49.545°, P-value <0.00001).

Matrix correlations

The matrix correlation between the covariance matrices for the symmetric and asymmetric components of variation for the complete dataset was 0.638 (P-value: < 0.0001; 10.000 perm.). Below are presented the pairwise matrix correlation results for the different cranial categories for both the symmetric and asymmetric component of variation.

Antero-posterior v/s non-deformed: 0.618 (P-value: < 0.0001; 10.000 perm.).

Antero-posterior v/s oblique: 0.592 (P-value: < 0.0001; 10.000 perm.).

Non-deformed v/s oblique: 0.750 (P-value: < 0.0001; 10.000 perm.).

Antero-posterior v/s non-deformed: 0.754 (P-value: < 0.0001; 10.000 perm.).

Antero-posterior v/s oblique: 0.745 (P-value: < 0.0001; 10.000 perm.).

Non-deformed v/s oblique: 0.873 (P-value: < 0.0001; 10.000 perm.).

Discussion

The PCA showed that even though there is morphological continuum ranging from one deforming type to the other (with non-deformed individuals relatively in between), it is possible to notice differences between the analyzed cranial categories. This was confirmed when running the pairwise PERMANOVAs and CVA analysis since highly significant differences between the two deforming styles and the non-deformed crania were observed. There was a slight overlap between the groups that was expected and coherent with the morphological continuum outlined above. It is interesting that in spite of the several proposed classifications to define deforming styles [72,93,111-115], the simple binomial sorting applied here followed by a numerical shape analysis, allowed a well-defined distinction between dissimilar deforming traditions. This result has taxonomic implications related to the a priori nature and the qualitative character of the classifications usually applied in bioanthropoloy and archaeology to give diagnostics about specimens’ provenance and grouping [116]. An almost inevitable consequence of this typological approach is the over-representation of categories and subcategories defining each group of objects [117]. Future studies should test whether more complex classification schemes have an enhanced classificatory power when compared to the simplified grouping used in this research. In addition, it is interesting that the two deformed cranial categories display more variation along the two first PCs as compared to the non-deformed crania. Although these differences in variation are not completely unexpected, they represent a good example how developmental plasticity acts on the human skull. Environmental perturbations such as the ACD influence growth and development. Hence, the application of these deforming devices during early years of post-natal life is likely related to an increase in morphological variance in the two deformed groups, explained in turn, as has been observed by several authors, by the increase of the vertical development of the maxilla [118], by changes in basicranial shape and in cranial base angles, as well as in the relative position of the mandibular fossae [92], significantly greater frequencies of lambdoid ossicles, and more lambdoid wormian bones [119], and a significant greater level of bilateral asymmetry [120,121], compared to the non-deformed group.

The null hypothesis stating that there was no modular behavior on the skull based on two different developmental origins was rejected. For the symmetric component of variation, the CR coefficient was significantly lower than one, thus suggesting that there was a degree of independence between the two modules. The asymmetric component also showed a lower and significant CR value, which suggest a strong degree of independence between the two modules. Thus, there is support for the hypothesis that the human cranium displays significant modularity when compared to the null hypothesis of no modular structure. It is relevant to keep in mind that modularity and integration are not two ends of the same continuum, thus it is perfectly possible to have both modularity and integration in the same dataset as in the present case. Modularity implies greater covariation of variables within modules than between them, whereas integration means that two modules are more correlated than it would be expected from chance by random arrangements of pairs of observations from each module. Hence, it is possible to have both integration between modules and modularity within them.

Regarding the static morphological integration of the human cranium, it seems that the different categories analyzed here show relatively similar overall levels of integration as observed in their r-PLS values. Nevertheless, the oblique deforming style appears to constrain the strength of static covariation between the viscero and the neurocranium in a more defined direction as compared with both the antero-posterior crania and the non-deformed skulls (Table 4). Oblique skulls exhibit a covariance greater than observed on the antero-posterior and the non-deformed individuals for PLS1. This could mean that the normal growing vectors of the skull are at least partially canalized towards the inner portion of the facial skeleton due to the application of the deforming device, consequently increasing the covariation with the calvarium. Interestingly the PLS1 for the oblique skulls (symmetric component) accounted for 75.467% (P-value: <0.0001) of the covariance of the sample, while the PLS1 for non-deformed skulls accounted for only the 51.948% (P-value: <0.0001) and 52.957% (P-value: <0.0001) for the antero-posterior deformed crania. These latter categories exhibit a covariation pattern that is more distributed in other directions (i.e. axes) as compared with oblique skulls. One measure of the morphological integration level is the correlation between the first PLS scores of the analyzed blocks, which was higher for the oblique skulls as compared to the other crania (Table 4) [20,122]. It is evident from these results that the oblique deforming style produces a more directed static covariation pattern between the viscero and neurocranium (Fig 6C). The skulls at the top of the correlation plot have more elongated and oblique neurocrania associated to more prognathic faces, while the skulls at the bottom have more rounded cranial vaults. These results are in agreement with previous studies that have found that the ACD effect is not confined only to the braincase but also affects the skeletal face [81,83]. On the other hand, the correlation plot of the PLS1 for non-deformed skulls showed that crania at the top of the plot exhibited more flattened neurocrania, while at the bottom the skulls showed a more rounded cranial vault. In the case of the antero-posterior deformed skulls, those individuals with more flattened occipitals were at the bottom of the plot, while those at the top showed more inclined cranial vaults. PLS2 also exhibits some interesting patterns (S2 Fig). Whereas the PLS1 shows, in a similar fashion as the PCA, a morphological continuum ranging from one deforming type to the other (with non-deformed individuals in between) (Fig 5), the PLS2 correlation plot of the complete dataset (symmetric component) shows a total overlap between the categories (S4A Fig). From morphological perspective, the crania at the top of this correlation plot display more elongated and laterally expanded neurocrania with a depressed area between the frontal and the parietals, whilst the skulls at the bottom show more rounded and vertically expanded vaults and flatter faces. PLS2 correlation plots for each one the cranial categories (symmetric component) can be found in S2B, S2C and S2D Fig.

Conversely, the developmental integration analysis showed that non-deformed skulls exhibit an overall lower level of integration as compared to the deformed skulls as shown by their r-PLS values. This might mean that the deforming styles increase the covariation between the asymmetric components of the neuro and viscerocranium. This means that an increased asymmetry in one of these cranial parts would generate a stronger increase of the asymmetry of the other skull compartment as compared to the non-deformed skulls, or a decreased asymmetry in one block would produce a stronger decrease of the asymmetry of the other cranial portion. In addition, the antero-posterior deformed crania showed a higher level of developmental integration in their first PLS axes (PLS1: 74.488% of explained covariance; P-value: <0.001; 1,000 perm.), which suggests that this particular deforming style directs the covariation between the two cranial blocks in a more directed way (i.e. it is less distributed in other directions of covariation).

ACD affects the morphological integration of the skull, showing a particular covariation trend depending on the deforming style [92]. It seems that the different deforming types produce different covariation patterns between the viscero and neurocranium as observed by the total covariance explained by each PLS axes for each one of the shape variation components under analysis. One plausible explanation is that the different deforming devices exerted different force vectors (both in direction and magnitude) on the neurocranium. These differences in force vectors were generated probably by the different materials of the deforming devices (i.e. wood, fibers, ropes, pads combined with flatted stones, etc.) and by the exact anatomical location where the deforming device was placed.

Another perhaps more important explanation to the different degree of integration between modules in deformed and non-deformed crania, is the change in the normal loads acting on cranial bones during ACD. There is an important role of mechanical loads in bone growth and development [123-128]. In the human cranium, the most noticeable source of loads is mastication, and its effect has been shown to be more important in the lower face [16]. The remaining cranial structures including the neurocranium reflect other non-load-bearing aspects such as population history and climate [129,130] in a stronger manner compared to the mandible [131,132]. ACD could act by affecting this modular pattern of load distribution, increasing the strength of the covariation between the neurocranium and the face in a more constrained direction by re-distributing intracranial pressure during development. This new, altered pattern of load likely impacts on a developing face that is now not only under the effect of the mastication and the functional matrices, thus allowing the deformation to reflect at a facial level.

When quantitatively testing if the observed covariation patterns were similar or not depending on the analyzed cranial category, the angular comparison between the PLS vectors shown that the majority of them were dissimilar (i.e. their angular difference was near 90°). However, there were some significant similarities between some PLS axes that could be related to the fact that despite the differences in the intensity of the covariation, the overall covariation pattern could be relatively similar. The angular comparison analyses also showed that there was a significant relationship between the PC’s and the PLS axes, thus indicating that the patterns of observed variation are in concordance with the observed covariation patterns between viscero and neurocranium. This could mean that the observed variation in cranial shape is due to the morphological integration of the skeletal face and cranial vault.

Finally, the matrix correlation analysis between the symmetric and asymmetric components of variation showed a moderate-to-high value. When carried out the matrix correlations by separating the different cranial categories it was observed that for both the symmetric and asymmetric components of variation, non-deformed skulls are more similar to the oblique crania. This seems to indicate that in spite of the increased static covariation observed in PLS1 of the oblique crania, they are more similar to the non-deformed condition. These results are concordant with the CVA results that showed that antero-posterior deformed skulls are the most different of the three cranial categories.

Conclusions

The results from this research show that there is a modular organization of the human skull (i.e. neuro and viscerocranium). Furthermore, the present results show that the strength of the morphological integration between the neurocranium and viscerocranium is differentially augmented depending on the applied force vectors on the skull (i.e. oblique deforming style). Compressive forces onto the parietal bones (i.e. oblique ACD) increases the static morphological integration between these two anatomical regions, while compressive forces onto the occipital and frontal bones (i.e. antero-posterior ACD), increases the developmental integration of the skull. Although the underlying cause of this phenomenon is still unknown, it could be related with the specific mechanisms constraining the normal expansion of the brain and how this affects the normal growth and development of the skull. Further analyses are required to get a better insight of the possible effects of ACD on human biology. One interesting approach would be to use the present results to carefully design a biomechanical simulation of the growing skull while simulating compressive forces as proxies for the different deforming devices.

Sample used in this study.

Individual numbers and complete repository information, including museum name and geographic location.

(CSV)

Click here for additional data file.

Intra-observer error Procrustes ANOVA.

Mean squares values, degrees of freedom (df), F statistic, and statistical significance (P value) of intra-observer error.

(CSV)

Click here for additional data file.

Posterior probabilities from the cross-validated CVA.

Linear Discriminant Analysis using the Procrustes coordinates of the first 11 PCs after a broken-stick model for 3 classes /'Antero-Posterior', 'Non.deformed', 'Oblique'/.

(CSV)

Click here for additional data file.

CR coefficients.

CR coefficients obtained from permutation tests (999 rounds) of alternative partitions of a) the complete dataset, b) antero-posterior deformed skulls, c) non-deformed crania and d) the oblique sample, with the observed CR coefficients designated by a red arrow. Both the symmetric and asymmetric components of shape variation were analyzed.

(PDF)

Click here for additional data file.

PLS2 correlation plots (symmetric component).

a) the Complete dataset, as welll as for the different cranial categories under analysis b) Antero-posterior; c) Non-deformed and; d) Oblique.

(PDF)

Click here for additional data file.