*2.5. Morphometric Mapping*

The SMA asymmetry values were obtained using the dominant asymmetry equation for all 21,960 (360 × 61) landmarks, and the results for each paired humeri were deposited in a matrix with 61 rows (sorted by the order of cross-sections) and 360 columns (sorted by the order of directions). These matrices were then visualized as morphometric maps to display the distribution characteristics of bending rigidity asymmetry along the proximodistal humeral diaphysis (Figures 1 and A1). The asymmetry values of J35 and J50 for all individuals were also calculated using the same equation.

**Figure 1.** The positional and directional correspondence between humeral external structure, diaphyseal cross-sections, and morphometric map exhibiting bending rigidity asymmetry. Abbreviations for anatomical terms are as follows: prox: proximal; mid: middle; dist: distal; lat: lateral; post: posterior; med: medial; ant: anterior.

#### *2.6. Methods to Estimate the Variation of Humeral Biomechanical Asymmetry*

To explore the variation in humeral asymmetry patterns in modern humans, 40 individuals were divided into sub-groups defined by sex and population. The three populations, which varied in geographic location, chronological age, and subsistence pattern, were supposed to vary in their habitual behaviors, so population was set as one variable. Sexual division of labor is an important issue when discussing historical populations, and the sexual dimorphism of humeral asymmetry can be affected by nonbehavioral factors such as genetics or hormones [27]. Therefore, sex was set as another variable. Mean morphometric maps exhibiting SMA asymmetry values for each sub-group were qualitatively compared. Additionally, a two-way multivariate analysis of variance (MANOVA) was conducted to quantitatively test whether sex and/or population were significant sources of variation. When fitting the regression model for MANOVA, SMA asymmetry values at all landmarks were set as the dependent variables, while sex and population were set as the independent variables with interaction. Customized in-house scripts, mainly sourced from R package 'geomorph' and 'RRPP', were utilized to conduct MANOVA [53,54]. In addition, the coefficients of variation (CV) for SMA asymmetry values at all landmarks were calculated in sub-groups and visualized by morphometric maps to exhibit intra-group variation characteristics. Only sub-groups defined by sex or by population were included in this analysis to reduce the impact of outliers.

#### *2.7. Methods to Test the Representativeness of J35 or J50 Asymmetry*

The reliability of using J35 or J50 asymmetry to represent the overall humeral asymmetry was tested using several statistical methods. First, a multivariate regression model was built on all specimens to statistically test the degree of correlation between overall SMA asymmetry and J asymmetry. When fitting the model, the SMA asymmetry values at all landmarks were set as the dependent variables, and the J35 or J50 asymmetry value as the independent variable. Customized in-house scripts, mainly sourced from R package 'geomorph' and 'RRPP', were utilized to carry out this fitting [53,54]. Second, to investigate the association of every SMA asymmetry value and J asymmetry value across the entire humeral diaphysis, the correlation coefficients between each SMA asymmetry value and J35 or J50 asymmetry value (CC35 and CC50) were calculated within sub-groups. The same protocols for visualizing SMA asymmetry values were applied to CC results to generate morphometric maps. The CC morphometric maps of sub-groups were qualitatively compared to reveal inter-group variations.
