**3. Results**

The asymmetry in polar moment of inertia, J, calculated at 40% of the biomechanical length shows a general trend of lateralization ranging between 0.36% and 10.37%. In all but 4 individuals, J is larger on the right side (Table 2). There are no statistically significant differences between occupation and sex group means (Figure 3), as determined by two-way ANOVA of J among sexes or occupations.

The first two PCs of the symmetric component of cortical thicknesses account for 78% of the total variance (PC1 = 72.33%; PC2 = 5.79%) (Figure 4). On average, the two sexes are separated along PC1 with the females towards the positive limit, and males, the negative (t14.48 = 4.66, *p* < 0.01). The three occupation groups largely overlap but with the "building", "army" and "desk" groups approximately distributed in this order from negative to positive limits of PC1.

**Table 2.** Description of the sample and calculation of lateralization. For all the individuals, we report age, weight, height and occupation. We calculated the biomechanical length of the humerus and the polar moment of inertia, J, at 40% of the biomechanical length on both sides (JL and JR) and the degree of lateralization, JLAT%, between the two sides. Values of J are multiplied by 10<sup>3</sup> (Ruff 2000). D = working at desk job; B = working in building/mining companies; A = working in the armed forces.


**Figure 3.** Violin and box plots showing the percentage of right lateralization pooled by sex (**left**) and occupation (**right**, A = "army", B = "building", D = "desk"). Violin plots illustrate the density distribution of the data.

**Figure 4.** PCA of the symmetric component. Circular and square dots represent male and female individuals, respectively. Violet = "Army", orange = "Building" and green = "Desk". The double-sized circles and squares show the mean values of PC1 and PC2 pooled by sex and occupation. In the morphometric maps, the diaphysis is unrolled from the anterior border (on the left of the x-axis) and follows the medial, posterior, lateral and anterior borders. Violet and green respectively indicate greater or smaller values of cortical thickness. The tops and the bottoms of the morphometric maps correspond to the proximal (80% of the biomechanical length) and distal (20% of the biomechanical length) parts of the diaphysis.

The visualisations of the morphometric maps represented by the extremes of PCs 1 and 2 highlight a general increase in cortical thickness at negative values of PC1. PC2 represents a different pattern of thickening/thinning of the humeral diaphysis. With increasingly negative values of PC2, the cortical bone is thicker in the mid and distal portion of the medial and anterior margins and between the posterior and lateral margins. Conversely, with more positive values, the proximal portion of the diaphysis is thicker anteriorly (Figure 4).

The PCA of the asymmetric component (Figure 5) indicates how the cortical thickness of the entire diaphysis differs from symmetry. The distance of points from the origin indicates the degree of asymmetry represented by the first two PC scores. Following Mardia et al. [46], we found that directional (mean difference between sides) and fluctuating (average differences from the symmetric mean) components account for 16.00% and 84.00% of the total variance, respectively.

**Figure 5.** PCA of the asymmetric component shown in two separated plots for the left (**A**) and right (**B**) sides. Arrows connect points representing the left (**A**) and right (**B**) sides of each individual from the origin (i.e., zero asymmetry between left and right side). Violet = "Army", orange = "Building" and green = "Desk". In the morphometric maps (**C**), violet and green palettes respectively indicate larger and smaller values of cortical thickness. In these, the diaphysis is unrolled from the anterior border (on the left of the x-axis) and follow the medial, posterior, lateral and anterior borders. Violet and green respectively indicate greater or smaller values of cortical thickness. The tops and the bottoms of the morphometric maps correspond to the proximal (80% of the biomechanical length) and distal (20% of the biomechanical length) parts of the diaphysis.

PC1 explains 41.05% of the total variance. This axis describes generalised thickening or thinning of the cortex as seen in the morphometric maps representing the extremes of this PC. Scores on PC1 indicate that the right-sided cortex tends to be thicker than the left (with a few exceptions, plausibly explained by handedness). PC2 (7.47% of the total variance) shows a different pattern of asymmetry. The morphometric maps indicate that, from positive to negative limits, this PC represents posterior and anterior thinning of the diaphysis (Figure 5).

While the differences between sexes (F = 18.40, Df = 1, *p* < 0.01) in asymmetry are significant, as indicated by two-way ANOVA (Figure 6), there are no statistically significant differences in asymmetry between the occupation groups.

**Figure 6.** Violin plots of the total length of the displacement vectors of asymmetry in relation to sex (**left**) and occupation (**right**, A = "army", B = "building", D = "desk"). Violin plots illustrate the density distribution of the data.

The lengths of the vectors in Figure 5 indicate the magnitude of asymmetry of cortical thickness, represented by the first two PC scores. The sum of the vectors from the entire matrix of PC scores represents the overall magnitude of asymmetry. In this sample, it is correlated with the index of lateralization (JLAT%), calculated using the polar moment of inertia, J, (correlation = 0.58, *p* < 0.01).

In multivariate regressions, body weight (R<sup>2</sup> = 0.19, b = 0.22, *p* < 0.01), height (R<sup>2</sup> = 0.22, b = 0.31, *p* < 0.01) and biomechanical length (R<sup>2</sup> = 0.16, b = 0.14, *p* < 0.01) significantly predicted asymmetry, while age and occupation score are not statistically significant predictors (Figure 7).

A variance partitioning analysis was performed to evaluate the percentage of variance in asymmetry associated with weight, height, age and occupation score. The combination of all tested variables explained 24.72% of the variance in asymmetry calculated from the PC scores of asymmetric component. Weight (2.09%) and height (2.39%), and their interaction (16.31%) explain the largest portion of asymmetry (Figure 8).

**Figure 7.** Map of R<sup>2</sup> and beta coefficients calculated from multivariate regression of the asymmetric component on independent variables of interest. Cortical thickness values were rank-transformed. In the first row, maps indicate which regions of the diaphysis show asymmetric variation in thickness with age, weight, height and biomechanical length. The R<sup>2</sup> range is reported using a rainbow palette. White cells indicate statistically insignificant relationships. In the second and third rows, the beta coefficients from the multivariate regressions are mapped on the left and right sides, respectively. Warm colours describe an increase in cortical thickness, and cold colours indicate a reduction.

**Figure 8.** Variation partitioning Venn diagram. Variation in asymmetry expressed as adjusted R<sup>2</sup> explained by weight, height, occupation and age and their intersections. Significance codes: \*\* = 0.001. Outside the Venn diagram are reported the total explained variances by each variable, taking into account the interactions between them. Note, the sum of the explained variance is not 100% (artefact of the variance partitioning algorithm due to the calculation of negative adjusted R2).

> Multivariate regression was used to assess the relationship between asymmetry in cortical thickness and the independent variable of interest (i.e., weight, age, height, biomechanical length and occupation scores) at each cell of the morphometric maps (i.e., the cortical thickness measured at each semilandmark). At each cell, the explained variance (R2) and slope (Beta coefficient) can be mapped to visualise the relationship between cortical map asymmetry and the independent variables. Maps can be drawn to represent the (exactly opposite) effects of these relationships on the right or left sides. Such maps are presented in Figure 8. The maps of R<sup>2</sup> indicate that each independent variable is associated with a different pattern of asymmetry, with different localised regions showing an association with each variable. Unsurprisingly, regressions of height and biomechanical length produce the most similar diagrams. On average, the slopes (beta coefficients) of asymmetry of cortical thickness on the independent variables indicate that cortical thicknesses tend to be greater in the right arm (likely this is a mostly right-handed sample). The maps of Figure 8 show values only for those regions where the regression is significant. Weight, height and mechanical length are associated with a larger rate of increase in cortical thickness in the right (usually the dominant) compared to the left humerus. In contrast, with increasing age, cortical thickness decreases more slowly in the dominant arm.

## **4. Discussion and Conclusions**

Studies of patterns of lateralization in archaeological populations often suffer from the lack of a reference sample with known loading history to contextualise the findings. Additionally, the assessment of lateralization is commonly limited to the analysis of one or just a few levels along the diaphysis of paired long bones. The publication of the NMDID [42] offers the prospect of directly relating skeletal form (total body CT-scan) with known biological profile and loading history (metadata with 60 variables). The time and effort in gathering skeletal data is much reduced by the functions available in the R

package *morphomap* [41], a recently published toolkit providing functions to extract from CT data, the segmented long bone of interest and based on that. Here, we further extend *morphomap* to visualize and analyze asymmetry in paired long bones. Specifically, the new implementation: (i) performs a PCA on the symmetric and asymmetric component of form variation; (ii) creates morphometric maps of symmetric and asymmetric variation on single individuals or on entire samples from PC scores; (iii) calculates the total magnitude of asymmetry of cortical bone distribution, quantifying deviations from symmetry.

PCAs of both symmetric and asymmetric components indicate that cortical bone distribution differs between sexes, but not between the occupation groups considered in our analyses. On average, male individuals possess thicker and more asymmetric humeri than females. All our measurements of magnitude of asymmetry (cross-sectional geometry and vector lengths from PC scores) present a general pattern of right lateralization. This finding is consistent with previous studies indicating 90% preference for right-handedness in modern humans [47–49]. Further, the analyses of morphometric maps indicate that males are more asymmetric than females. However, since males are larger on average, their more asymmetric cortical thickness might be guided by allometric effects. In contrast, the index of lateralization based on cross-sectional geometry (JLAT%) does not significantly differ between sexes. This contrast may be due to the fact that JLAT% is calculated at a single level along the diaphysis (40%v of bone biomechanical length), whereas the calculation of absolute lateralization from PC scores takes into account the entire diaphysis. In fact, JLAT% calculated at 70%, 75%, 76%, 78% and 79% of bone length is statistically different between sexes.

Regression analyses on morphometric maps show that, as body weight, height and longbone biomechanical length increase, so does asymmetry. In contrast, with increasing of age, asymmetry decreases (i.e., the oldest individuals are less asymmetric than the younges<sup>t</sup> individuals). Interestingly, loading history (occupation scores) does not affect the pattern of asymmetry. The main effect of body proportions (weight and height) as a key factor in determining the degree of asymmetry is also confirmed by the partitioning of variance analysis performed on values of total asymmetry calculated from the PCA.

This study was able to test the hypothesised relationships between loading history and cortical bone distribution in the humeral shaft using the unique and extensive collection curated in the NMDID [42]. Despite the quality of these data, our analyses do not show a significant association between occupation and asymmetry. This analysis was confined to occupations that would be expected to lead to extremely different occupational loading histories (army, building vs. desk) in order to emphasise any effect of occupation. This suggests either that there is no strong difference in the effect of these occupations on loading history, or that occupational history does not reflect the full loading history, and that our categorisation by occupation is inadequate to describe individual loading history. Further studies are required to clarify this finding, which is potentially of grea<sup>t</sup> importance in archaeological and forensic contexts.

**Author Contributions:** Conceptualization, A.P., D.M., G.M. and P.O.; data curation, A.P. and C.Z.; formal analysis, A.P., C.Z., I.M. and A.M.; funding acquisition, A.P. and P.O.; investigation, A.P. and C.Z.; methodology, A.P. and P.O.; project administration, A.P. and P.O.; resources, A.P. and C.Z.; software, A.P., P.R. and P.O.; validation, C.Z.; visualization, A.P. and A.M.; writing—original draft, A.P., C.Z., P.R. and P.O.; writing—review and editing, C.Z., I.M., A.M., P.R., D.M., G.M. and P.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from the European Union's Horizon 2020 research and innovation program (H2020 Marie Skłodowska-Curie Actions) under gran<sup>t</sup> agreemen<sup>t</sup> H2020-MSCA-IF-2018 No. 835571 to Antonio Profico and Paul O'Higgins.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of the University of York, Department of Archaeology.

**Informed Consent Statement:** Not applicable. **Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors are grateful to two anonymous reviewers for helpful comments.

**Conflicts of Interest:** The authors declare no conflict of interest.
