## 2 Materials and Methods
### Principal Component Analysis Applied to Asymmetry  
To compare the asymmetry of the specimens, a mirrored and re-labelled copy of each configuration was made. From these two copies, the individual asymmetric component was found by computing the differences between each landmark configuration and the reflected version. A modified version of Principal Component Analysis (PCA), defined as asymmetric PCA, was used to explore the asymmetries. PCA is a method that looks at complex multi-variate data and can be used to find patterns between variables and quantify this [41]. In an asymmetric PCA, the \"origin of the axes represents perfect symmetry\", and the specimen positions were represented by vectors departing from the origin, with their length and orientation indicating respectively the magnitude and direction of each specimen's asymmetry [23; 32]. Therefore, this PCA analyses only the asymmetric component of shape variation. This method allows for conclusions to be drawn about the nature of the asymmetry and the asymmetric variation between specimens. The proportions of directional and fluctuating asymmetry were then estimated using a sum of squares decomposition of the total asymmetry.  
Furthermore, a PCA of the symmetric component was performed to visualise the symmetric shape variations within the sample. The symmetric shape was obtained by averaging each landmark configuration and its superimposed relabelled counterpart [38; 40].  
Figure 3: Process of Procrustes superimposition, modified and adapted from [PERSON] and Gunz [39]. 2009, [PERSON] and [PERSON]. This proceeds by standardising the size, position, and orientation of the landmark configurations through registration. Squared distances between equivalent landmarks obtained from all specimens are minimised by centering (translating), scaling, and rotating (aligning) the individual’s landmark configurations using their centroid (the arithmetic mean of all landmark coordinates). Scaling is carried out in order that the centroid size (CS or Csize—the square root of the sum of squared distances of a set of landmarks from their centroid) of each configuration is scaled to 1 by dividing the raw landmarks by their centroid size. Configurations are then superimposed at their centroids and iteratively rotated with respect to each other to minimize the sum of squared distances of the specimens from the mean shape. In the first iteration, the specimens are aligned to an arbitrarily chosen one and once all configurations are fitted to this, the rotation minimises distances between equivalent landmarks. Once all specimens are fitted, the sample mean of each coordinate is computed. In subsequent iterations, the configurations are fitted to the mean, which is recomputed and used in the next iteration. The algorithm stops when the sum of residuals from the mean reaches a minimum, usually after 3-5 iterations. The resulting registered landmark coordinates are known as “shape coordinates” and can be submitted to statistical analysis.  
Shape variations of the first and second Principal Component axes (PC1 and PC2) for both the asymmetric and symmetric PCA were visualised using Thin-Plate-Splines (TPS) [42, 43]. First, the 3D mesh of the cranium showing the smallest Procrustes distance from the mean shape was chosen for the visualization and warped onto the mean shape. To visualise the shape variations, the mesh was warped using TPS onto the maximum and minimum scores of the PC1 and PC2 for both the symmetric and asymmetric PCAs [23, 38].