Experimental evidence that physical activity affects the multivariate associations among muscle attachments (entheses)

Muscle attachment scars (entheses) comprise the areas of the bone where muscles or ligaments attach (Benjamin et al., 1986). In paleontological and bioarchaeological sciences, they are frequently utilized as a basis for reconstructing habitual physical activity patterns among past populations (Foster et al., 2014; Schrader, 2019). Previous research has focused either on morphological changes within entheseal surface areas (osteophytic or osteolytic traits; e.g. Mariotti et al., 2004; Villotte, 2006; Villotte et al., 2010; Henderson et al., 2017) or calculations of entire entheseal three-dimensional (3D) morphology (e.g. Zumwalt, 2005; 2006; Karakostis and Lorenzo, 2016; Wallace et al., 2017; Karakostis et al., 2017, 2018a,b). However, most of the previous methods for analyzing entheses have been associated with either low or untested measuring precision (for entheseal changes, see Davis et al., 2013; Wilczak et al., 2017; Villotte et al., 2016; for overall 3D morphology, see Noldner and Edgar, 2013), absence of multivariate statistical procedures (see Milella et al., 2015; Karakostis and Lorenzo, 2016; Karakostis et al., 2017, 2018a) and inadequate consideration of other factors affecting entheseal form such as age and body size (e.g. Williams-Hatala et al., 2016; see also Foster et al., 2014; Henderson et al., 2017; Schrader, 2019).

Recently, we proposed a new approach for analyzing entheses that is based on precise 3D measuring protocols (Karakostis and Lorenzo, 2016; Karakostis et al., 2017). This approach strives to control for the complex factors of inter-individual entheseal variation (see Foster et al., 2014; Schrader, 2019) by focusing on the multivariate proportions among attachments within each individual, rather than comparing single entheseal structures across individuals with distinct characteristics. It was recently applied on a unique human skeletal sample (the Basel–Spitalfriedhof collection) with detailed lifelong occupational documentation (exact jobs, positions and durations documented for each individual) as well as life histories (genealogical data, official medical records, socioeconomic details and more) (Hotz and Steinke, 2012), and a clear separation was found across individuals with distinct occupational characteristics (Karakostis et al., 2017). Specifically, that analysis revealed that long-term construction workers exclusively presented a power-grasping hand entheseal pattern, whereas lifelong precision workers only showed a precision-grasping pattern involving the thumb and index finger (Karakostis et al., 2017). Ultimately, based on the biomechanical literature, co-ordination of the particular muscle synergy group represented by each of the two observed entheseal patterns could be associated only with power or with precision grasping movements (e.g. Marzke et al., 1998; Clarkson, 2000; Huesler et al., 2000).

Furthermore, we recently validated the capacity of our entheseal method to detect repetitive muscle contraction in life, relying on a previously published experimental model involving electrically induced muscle stimulation in Wistar rats (Vickerton et al., 2014; Karakostis et al., 2019a). The results of that experimental analysis demonstrated that our multivariate 3D approach can accurately identify which individuals and muscles were systematically stimulated in the experiment, under blinded study conditions involving different research institutions. Nevertheless, the experimental model used in that previous analysis did not involve real-life activity patterns that rely on the coordination of certain muscle synergies. Therefore, despite its compelling results, direct experimental evidence of a causal relationship between physical activity patterns and entheseal form is still lacking in the literature. In fact, previous attempts to experimentally test for such relationships have failed to support the assumption of such a direct link (Zumwalt, 2006; Rabey et al., 2015; Wallace et al., 2017). Particularly, these previous works have reported the absence of a bivariate association between activity patterns and entheseal variation. However, owing to the methodological limitations explained above, we do not consider the conclusions of those previous experimental studies as definitive (for a more extensive discussion, see Karakostis et al., 2018a).

Here, in order to assess the validity of our new method for investigating variation in entheseal morphology directly related to physical activity patterns (Karakostis and Lorenzo, 2016; Karakostis et al., 2017, 2018a), we applied our method to the skeletons of animals used in a previous experiment that did not detect a link between entheseal changes and physical activity (Wallace et al., 2017). That study was co-authored by two of us (I.J.W. and N.K.) and found that the applied activity regime (running) had no effect on the morphology of the lateral epicondyle of the femur (i.e. attachment area of the lateral gastrocnemius, an ankle extensor), based on an experimental model involving 10 turkeys (five downhill runners and five controls). Turkeys were chosen for that experiment because bipedal birds are considered to be appropriate models for human locomotion, as they exhibit relatively similar patterns of limb movement (Gatesy and Biewener, 1991), basic body center-of-mass dynamics and pendular exchange (Biewener and Daley, 2007). Although Wallace et al. (2017) found that the employed activity regime led to significant changes in diaphyseal and trabecular bone morphology, their analysis did not detect significant alterations in entheseal topography. The variables representing entheseal morphology described surface curvature (‘Dirichlet normal energy’), surface relief (‘relief index’), and surface complexity (‘orientation patch count rotated’). In contrast to those results, here we show that our newly developed multivariate method is actually able to detect differences in entheseal morphology due to physical activity. Our study thus provides experimental validation of our method and the first direct evidence of a causal relationship between physical activity and entheseal morphology.

All experimental procedures (originally published in Wallace et al., 2017) were reviewed and approved by the Institutional Animal Care and Use Committee of Brown University.

For a detailed description of the experimental design, see Wallace et al. (2017). Briefly, fifteen 1-year-old female Eastern wild turkeys (Meleagris gallopavo Linnaeus 1758) were purchased from a licensed breeder, group-housed in a climate-controlled room with rubber flooring (6×10×3 m; width×length×height) on a 12 h:12 h light:dark cycle with food and water freely available. Animals were randomly assigned to one of three activity groups: uphill runners, downhill runners and controls (N=5/group). Running groups were exercised on either an inclined or declined motor-driven treadmill (uphill and downhill runners, respectively) at a speed of 2.5 m s⁻¹ for 30 min a day, 4 days a week, over a period of 10 weeks. The treadmill slopes were gradually adjusted from 3 to 8 deg during the first week of training to 20 deg during the last weeks of training, as per the capability of each individual. The number of strides was 137,087±9929 (mean±s.d.). Controls ran only for the initial 4 days of the experiment and only on a level surface. At the end of the experiment, animals were anesthetized using isoflurane (4%; 2 l min⁻¹ O₂) and, once in a deep anesthetic plane, euthanized with a sodium pentobarbital overdose (i.v.). Subsequently, their femora were extracted and scanned using a micro-computed tomography (microCT) scanner (Nikon XT H 225 ST) and a voxel size of 45 μm³. 3D iso-surfaces were then extracted from microCT images and saved in ‘.stl’ format.

The distinctiveness of bony muscle attachments presents high variability both among and within individuals of a species, with the areas of some entheses being more distinguishable than others (Foster et al., 2014; Schrader, 2019). In this study, we selected three entheses which were clearly defined and directly homologous across the femora of all 15 subjects (see examples in Fig. 1). These involved the fossa of the medial epicondyle of the femur (origin attachment of the muscle medial gastrocnemius), the insertion point of the gluteus primus (located in the proximal femoral shaft) and the medial supracondylar line (partial attachment of vastus medialis and adductor magnus) (Hudson et al., 1959; Hudson and Lanzillotti, 1964; Koenig et al., 2017). The muscles associated with these structures are coordinated during running locomotion in galliform birds (Hudson and Lanzillotti, 1964; Koenig et al., 2017). Therefore, if associations among muscle attachments reflect habitual muscle synergies, we expect the entheses of these muscles in the running animals to present a multivariate pattern that would distinguish them from those of the control animals.

After extraction of the 15 iso-surfaces, all 3D areas of entheses were measured without considering the activity group of each turkey subject. Each individual's 3D model was imported into the software package Meshlab (CNR-INC, Rome, Italy), where the three selected entheses were delineated on the bone surface using the highly repeatable imaging techniques introduced in our previous research (Karakostis and Lorenzo, 2016; Karakostis et al., 2017, 2018a,b). In brief, this methodology involves the use of various software-provided 3D imaging filters for delimiting muscle attachments on the basis of bone elevation and surface complexity (i.e. geodesic distances within entheses as well as calculation and color-mapping of the surrounding bone area's principal directions of curvature). After delineation of their borders on the bones, entheseal surfaces were isolated and measured in mm² (Karakostis and Lorenzo, 2016; Karakostis et al., 2018a). According to our previous intra-observer and inter-observer repeatability tests, the maximum mean precision error of this new measuring technique was 0.62% (Karakostis and Lorenzo, 2016; Karakostis et al., 2018a), while entheseal measurements are shown not to vary significantly among different 3D scanning technologies (computed tomography, laser-scanning and structured-light scanning) (see Karakostis et al., 2018a).

The resulting 3D surface measurements were subjected to principal component analysis (PCA), an exploratory technique that reveals multivariate patterns of variation in a sample, without using a priori group classification (Field, 2013). The correlation (rather than covariance) matrix was used, because the measurement scales of the three variables varied substantially (Table 1). Prior to the analysis, we verified that the dataset met all assumptions needed for PCA (see Field, 2013), which include normal distribution of the variables (based on Shapiro–Wilk tests, P-values ranged from 0.13 to 0.28), absence of significant outliers (based on the z-scores approach; see Field, 2013), minimum sample size requirements (five cases per variable) and linearity among variables. All three PCs extracted (summarizing 100% of the variance in the sample) were plotted, and individual cases were labeled by activity group.

In order to investigate whether body size or bone length was associated with the observed multivariate patterns, we applied a series of Kendall's tau-b tests to evaluate the strength of correlation between these parameters and each of the three PCs (i.e. a total of six comparisons performed). This non-parametric test, which is advisable for relatively small sample sizes (Field, 2013), was preferred because of slight violation of the assumption of approximate normality in PC2 scores (Shapiro–Wilk test, P-value=0.02). All statistical analyses were performed in the IBM SPSS software package (IBM Inc., Armonk, NY, USA; v. 24 for Windows).

The descriptive statistics for all three variables used are summarized in Table 1; Fig. 2 presents the values from each animal. These 3D surface measurements were obtained following protocols with thoroughly tested intra- and inter-observer repeatability (Karakostis and Lorenzo, 2016; Karakostis et al., 2018a; also see Karakostis et al., 2018b). Within each separate variable there is extensive overlapping among uphill/downhill runners and controls. This demonstrates that statistical attempts to compare single entheseal structures would likely not identify any marked difference among groups in entheseal 3D surface size.

By contrast, the resulting plot of the multivariate analysis (correlation PCA) on the same dataset showed a distinctive separation between controls and the animals exercised for 10 weeks (uphill runners and downhill runners), which broadly overlapped (see Figs 3 and 4 for 3D and 2D plots of the three PCs, respectively). The multivariate patterns responsible for this differentiation were represented by PC2 and PC3, which explained 30.64% and 16.82% of the total sample variance, respectively (Table 2). On PC2, uphill and downhill runners presented higher PC scores, which were associated with proportionally higher values for all three entheseal variables. However, the entheses of the gluteus primus and medial gastrocnemius muscles presented substantially higher PC2 factor loadings (0.73 and 0.62, respectively), whereas the correlation between this component and the medial supracondylar line was minimal (0.02, an almost ‘zero coefficient’; see Field, 2013). As a consequence, PC2 does not reflect overall entheseal size, which would require considerable factor loadings in all three variables analyzed (over 0.30; see Field, 2013). Similarly, PC3 (16.82%) explains variation in the relative size of the medial epicondyle (attachment area of the medial gastrocnemius) compared with the other two entheses. On this axis, three of the five control animals presented lower values than all the rest, showing a proportionally larger medial epicondyle of the femur (Fig. 4B). The minimal effect of individual size (body and bone dimensions) on all observed multivariate patterns is directly supported by the Kendall's tau-b correlation tests, which found no significant correlation between PC scores and the individuals' body mass or bone length (Table 3).

Finally, on PC1 (52.54%), there was extensive overlapping among the three cohorts. Nevertheless, it is worth mentioning that this component revealed certain differences between uphill and downhill runners (for biomechanical interpretations, see Discussion). In particular, four of the five downhill runners showed slightly higher scores than both uphill runners and controls. Based on the factor loadings (Fig. 3, Fig. 4A, Table 2), the higher PC1 scores of downhill runners correlated with a proportionally larger 3D surface area of the medial supracondylar line and the medial epicondyle (Table 2). At the same time, the mean 3D surface size of these two entheses – and especially of the medial supracondylar line – was also considerably larger in downhill runners (Table 1).

The results of our multivariate analysis (Fig. 3) showed that certain entheseal patterns (PCs) provided a clear separation between all 10 turkeys exercised for 10 weeks and the five controls. This is consistent with the results of our previous research on human individuals with detailed lifelong occupational documentation and life history, which found differences between long-term heavy manual laborers and precision workers of lower intensity (Karakostis et al., 2017; see also Karakostis et al., 2018a). It also agrees with our recent experimental results of controlled muscle stimulation in Wistar rats, which showed correlations with predicted changes in entheseal multivariate patterns (Karakostis et al., 2019a). Moreover, it accords with the conclusions of a recent experiment on adult mice demonstrating that artificial muscle paralysis can lead to significant bone loss (at the millimeter scale) in the associated entheseal surfaces (Deymier et al., 2019).

The fact that the three turkey groups overlapped extensively when the surface area of each enthesis was considered separately (Fig. 2) is also consistent with the results of past experimental research, which reported no significant association between the morphology of single entheseal structures and physical activity (Zumwalt, 2006; Rabey et al., 2015; Wallace et al., 2017). Rather than performing multivariate analyses to explore the underlying patterns among different entheses for each individual, those previous studies – including the experiment from which our animal sample derives (Wallace et al., 2017) – focused on comparing single entheseal structures directly across individuals and groups. It is the application of our multivariate method to the same animals that allowed us to identify linear combinations (PCs) of the original variables (surface areas) separating uphill/downhill runners and controls without any prior group classification (Fig. 3). Our novel approach therefore reconciles the conflicting results of previous studies.

Even though multivariate statistical procedures are considered to be a standard step in many biological sciences (e.g. genetics, paleontology, forensic anthropological methods for human identification, and others; Karakostis et al., 2014; Posth et al., 2017; Ioannidou et al., 2019; Harvati et al., 2019) as well as some skeletal markers of physical activity patterns (e.g. trabecular morphology; Scherf et al., 2013, 2016; Stephens et al., 2018), their application on entheses remains scarce in the literature. This is in spite of a few previous anthropological studies reporting intriguing results after multivariate analyses, either involving overall 3D morphology of entheseal surfaces (Karakostis and Lorenzo, 2016; Karakostis et al., 2017, 2018a) or ordinal scores describing entheseal changes (Milella et al., 2015). It must be emphasized that our method, which focuses on the underlying shared correlations among different entheses in each individual, aims at identifying muscle synergy groups closely related to certain physical tasks (e.g. running locomotion). Therefore, its biomechanical interpretations do not focus on the function of each muscle separately, as this would be potentially associable with a diverse range of physical movements. By performing multivariate statistics, we attempt to control as much as possible for the numerous factors affecting direct interindividual variation in entheses, as we rely on the proportions among muscle attachments within each individual rather than the morphology of each single entheseal structure. By contrast, the widely applied direct statistical comparisons (i.e. bivariate correlations or mean comparisons) across distinct individuals are likely more prone to the numerous and complex factors affecting direct interindividual entheseal differences. These may involve age, body size, bone dimensions, sexual dimorphism, systemic factors (hormones and nutrition), pathology, genetic variability, microscopic and macroscopic taphonomical processes, and complex interactions between these factors (Rauch, 2005; Foster et al., 2014; Henderson et al., 2017; Karakostis et al., 2017, 2018c).

Importantly, our results also showed partial separation between the two experimental runner groups along PC1: four out of five downhill runners showed marginally higher PC1 scores than uphill runners and controls (Fig. 3, Fig. 4A). Higher PC1 values reflect proportionally greater 3D areas of the medial supracondylar ridge (partial attachment area of the vastus medialis and adductor magnus muscles) and the medial epicondyle (enthesis of the medial gastrocnemius muscle) compared with the enthesis of gluteus primus (Table 2). As a consequence, the convex hulls of the two activity groups do not tend to be aligned (Fig. 4A). Interestingly, this tendency is also reflected on the raw 3D area measurements of downhill runners in the same entheses (see Fig. 2 and Table 1). The most notable difference is found in the medial supracondylar line, where downhill runners present an average surface area that is 10 mm² larger than in the other two groups (Table 1). This translates to a mean increase of approximately 40% for downhill runners, which does not seem to be associated with a considerably different standard deviation or error in the estimation of the mean (Table 1). It is likely that these differences between the two activity groups are due to muscle physiology. In downhill locomotion, the gravity-resisting muscles (i.e. the leg extensors; e.g. vastus medialis that inserts along the medial supracondylar line) are subjected to stretch when activated in order to dissipate gravitational potential energy. However, in uphill locomotion, the same muscles shorten while activated in order to elevate the center of mass. As a consequence, the leg muscles of downhill runners are expected to generate higher muscle forces, and therefore the entheses would experience greater peak stresses per stride compared with uphill runners. In this context, one might be tempted to apply statistical probability tests (e.g. a general linear model) for directly comparing PC scores between the uphill and downhill runners (Konow and Sanford, 2008). Nevertheless, in the dataset of the present study, as a result of violation of certain statistical assumptions (i.e. homogeneity of variances, normal distribution or linearity among variables; Field, 2013) as well as a relatively small sample size per group, we do not consider statistical significance testing to be a fully reliable avenue (for more discussion on the potential risks of statistical significance testing in certain cases, see Amrhein et al., 2019). By contrast, our dataset met all statistical assumptions required for performing a correlation PCA on the raw entheseal 3D surface calculations (see Materials and Methods). That observational approach revealed morphological variation between uphill/downhill runners and controls without taking into consideration any a priori group classification.

In spite of the above observations, there is substantial overlap between uphill and downhill runners in the PCA (Figs 3 and 4). This could be related to the selection of entheses, given that all three correspond to muscles recruited for all running activities in galliform birds (Hudson et al., 1959; Hudson and Lanzillotti, 1964). Alternatively, one could wonder whether the resolution of biomechanical information provided by entheseal patterns can allow for the identification of such fine details regarding the nature of physical activity. To further investigate the limits of our method, as well as of entheseal variability in general, future experimental research could compare individuals subjected to more distinct activity regimes (e.g. climbing versus running tasks; see Rabey et al., 2015) and/or focus on diverse muscle synergy groups.

In conclusion, the results of our experimental study on turkeys demonstrate that it is possible to detect the performance of physical activity based on entheseal morphology, as long as entheseal analysis relies on precise 3D quantification and multivariate statistical procedures focusing on the proportions among entheses rather than the analysis of single muscle attachments. Because muscle attachments are the only bone areas directly associated with the musculotendinous unit, the positive results of this study are relevant for scientific fields focusing on reconstructing activity patterns based primarily on skeletal remains (e.g. paleoanthropological, paleontological, bioarchaeological or forensic anthropological studies). In particular, given that galliform birds have been proposed as appropriate locomotor models for human-like locomotion (Gatesy and Biewener, 1991; Biewener and Daley, 2007; Wallace et al., 2017), and based on the previously successful application of our method on uniquely documented human individuals (Karakostis et al., 2017, 2018a), we strongly believe that biomechanical interpretations based on entheseal multivariate patterns can be informative for human evolutionary studies.

We are grateful to Professor Thomas Roberts (Department of Ecology and Evolutionary Biology, Brown University) for donating the turkey bones that form the basis of this study. We would also like to thank Doug Boyer for microCT scanning the bones analyzed in this study. Also, for helpful discussions, we thank the participants of a symposium entitled ‘Evidence and methods surrounding the use of entheses to reconstruct muscle anatomy and activity’ held at the 2019 Annual Meeting of the American Association of Physical Anthropologists in Cleveland, OH, United States. We also thank two anonymous reviewers for their helpful comments and suggestions.