In forensic anthropology and bioarchaeology, the accurate construction of an unknown individual’s biological profile is essential. Hence, a large number of relevant works focus on methods for the accurate assessment of key parameters, such as sex and age. The developmental stages of the human skeleton have been tightly linked to the age of an individual [1, 2], especially in the earlier ages. As such, methods that assess the long-bone length together with the epiphyseal fusion state have been developed [3, 4], as well as methods that rely on dental evidence [5-7]. Similarly, the macroscopic degeneration of bones has been used to estimate age in adults, with the majority of traditional methods focusing mainly on the pelvic bones [8-10], the sternal end of ribs [11-13], and more advanced methods employing transition analysis on various skeletal age markers [14, 15].

Research has also been carried out regarding the spine and its association with aging. The vertebral column undergoes major morphological changes due to age degeneration and activity strain, which can be used as general indicators of aging [2]. The majority of age-related studies on the spine has focused on the observed morphological changes on the vertebral bodies and, more specifically, on the degree of the formed osteophytes [16-21]. The latter have been strongly associated with aging, as they are the result of the intervertebral disc degeneration and the subsequent osseous changes [22]. Other studies have used specific measurements to examine the effect of age on additional parameters of the vertebral body or other relevant skeletal structures [23-25].

Estimating the age-at-death of unknown skeletal remains based on the degree of osteophyte formation of the vertebral column has been extensively used in bioarchaeological and forensic contexts, especially when other skeletal age markers are missing or severely damaged. However, there is a limited number of studies focusing on age estimation from morphological properties of the vertebral body recorded as continuous variables, since the majority are using categorical scoring techniques based on macroscopic observation. Thus, this study proposes a number of geometrical variables deriving from the vertebral endplates of the 12th thoracic vertebra (T12), which was chosen due to its easily distinguishable morphology and high likelihood of remaining intact. The aim of this study is to examine whether the variables show significant association with age, using samples from two intertemporal European populations, and, if so, how aging affects the morphology of the vertebral endplates. The effects of sexual dimorphism are also evaluated and taken into account. Additionally, the proposed variables are examined for potential applications in forensic and bioarchaelogical contexts.

Sample description

The working sample consists of two European populations: a modern Greek sample and two Danish archaeological assemblages, one unearthed in Aalborg (ÅHM6093 AS 40/14), a medieval urban population (AD 1100–AD 1600) and another from Holmens-Kirke (NM538/2013 AS 4/14), also urban dated back to the 18th century. The Greek sample is part of the Athens collection, which is housed in the Department of Animal and Human Physiology at the National and Kapodistrian University of Athens. It comprises 168 adult individuals (93 males and 75 females) born between 1879 and 1965, whose age-at-death ranges from 18 to 99 years old (mean: 55.7 years; standard deviation: 18.6 years; skewness: 0.02) as documented in their death certificates [26]. The Danish archaeological samples are housed in the Department of Forensic Medicine at the University of Copenhagen and are pooled for the purpose of this study. Pooling aims for adequate sample size and any population differences between the two chronologically distant archaeological assemblages do not pose any constrains in the experimental design (analysis) of this study. In total, they comprise 107 adult individuals (56 males and 51 females) with an estimated age-at-death mean range between 21.5 and 82.5 years old (mean: 46.5 years; standard deviation: 18.5 years; skewness: 0.42). The age estimation of the Danish samples has been previously assessed using standard methods [27]. The distributions regarding the age-at-death of each population sample are shown in Fig. 1. It should be noted that the histogram for the Danish sample is based on the midpoints of the estimated age-at-death range of each individual. The working sample explicitly concerns individuals with no apparent pathological deformities and at least one endplate of the T12 being intact.

Figure 1. Age and sex distributions of Greek and Danish samples.

Data collection

The intact upper and lower endplates of the T12 vertebrae were photographed with a Canon EOS 100D camera and a Canon EF 40 mm f/2.8 STM fixed focal length lens. All images were shot from a stationary position with the camera perpendicular to the endplates, which were coplanar with the reference scale as shown in Fig. 2. The photographed endplates were imported in Autodesk 2012 software (Autodesk Inc., California, USA) as scaled raster images, and the coordinates of two two-dimensional (2D) polylines were recorded: 1) one polyline of the external boundary (red) outlining the outside edge of the ring of solid epiphyseal rim including any osteophytes and 2) a second polyline of the internal boundary (blue) separating the epiphyseal rim from the more porous central endplate (Fig. 2). Nine morphometric variables describing the morphological variation of the endplate’s top view surface were calculated from the polyline coordinates. These variables are listed in Table 1 along with their definition and abbreviation.

Table 1 . List of the morphometric variables evaluated on the lower and upper endplates of T12

Variable	Abbreviationa	Definition
Internal area	LiArea/UiArea	The area (in mm2) of the more porous central endplate as circumscribed by the blue outline in Fig. 2.
Internal perimeter	LiPer/UiPer	The length (in mm) of the boundary of the more porous central endplate as traced by the blue outline in Fig. 2.
External area	LeArea/UeArea	The area (in mm2) of the external boundary of the endplate including any osteophytes as circumscribed by the red outline in Fig. 2.
External perimeter	LePer/UePer	The length (in mm) of the external boundary of the endplate including any osteophytes as traced by the red outline in Fig. 2.
Internal-external area	LieArea/UieArea	The area (in mm2) of the epiphyseal rim of the endplate as circumscribed by the two outlines in Fig. 2.
Internal-external area ratio	LieAreaRatio/UieAreaRatio	The ratio of the area of the epiphyseal rim over the external area of the endplate.
Internal area perimeter index	LiArPerIndex/UiArPerIndex	The ratio {x∈R|0<x≤1} of the internal area to the squared internal perimeter multiplied by 4π. It reaches maximum when the circumscribed area is perfectly circular. This is a measure of smoothness of the perimeter of the more porous central endplate (blue outline).
External area perimeter index	LeArPerIndex/UeArPerIndex	The ratio {x∈R|0<x≤1} of the external area to the squared external perimeter multiplied by 4π. It reaches maximum when the circumscribed area is perfectly circular. This is a measure of smoothness of the perimeter of the external boundary of the endplate (red outline).
Internal-external area perimeter index	LieArPerIndex/UieArPerIndex	The ratio of the area of the epiphyseal rim to the squared external perimeter multiplied by 4π. This is a complementary measure of smoothness of the external boundary of the endplate (red outline).

T12, 12th thoracic vertebra. ^aThe initials L and U indicate the lower (caudal) or upper (cranial) endplate, which the variable corresponds to.

Figure 2. Polyline tracing of external (red) and internal (blue) endplate boundaries on the endplate’s surface.

Repeated measurements on both upper and lower endplates of 10 randomly selected vertebrae from the Greek population were performed for evaluating the intra- and inter-observer error by the first author (NG) and the second author (AB).

Data analysis

The intra- and inter-observer errors were evaluated in terms of the relative technical error of measurement (TEM%) for each morphometric variable, as well as its corresponding reliability coefficient R [28]. The expression of sexual dimorphism was also examined for each morphometric variable in both populations. This was done by means of one-way ANOVA or Kruskal-Wallis test according to the assumption of normality, which was tested with Shapiro-Wilk’s test, and effect size estimates ω² and ε² were calculated respectively. The correlation between the morphometric variables and age-at-death was explored with generalized additive models (GAMs) separately for males and females due to sexual dimorphism revealed from the previous analysis. All GAMs herein employ a mixed model approach via restricted maximum likelihood (REML) for estimating the smoothing parameter from the available data. Univariate GAMs were developed for every morphometric variable on the Greek population data, initially using the age-at-death as a predictor for the value of individual morphometric variables, hence particular traits of the endplate’s shape, and subsequently using the morphometric variables for predicting age-at-death. Multivariate GAMs were also constructed with selected morphometric variables, separately for the lower and upper endplates. This selection was through trial and error among the best performing univariate GAMs in terms of the higher deviance explained. The best performing multivariate GAMs were finally applied on the Danish population for predicting the age-at-death of each individual, which was cross-validated against the previously assessed age-at-death range estimates.

The statistical analysis including plots and results presented in this work has been implemented in GNU Octave and R programming languages and the entire code base including the raw data (original 2D polyline coordinates and replicates for measurement error) is available at https://doi.org/10.5281/zenodo.5776793. Use the ‘run_analysis.m’ script to reproduce this work from within the GNU Octave environment.

Intra- and inter-observer error results showed the relative TEM is below 1% and the reliability coefficient R is above 0.98 for every variable included in this analysis. This indicates that the measurements extracted from the photographs are very reliable. The shape variables (AreaRatio and ArPerIndex) performed slightly worse than the size variables (perimeter and area), since the former are direct derivatives of the latter and consequently the measurement error accumulates. Detailed values of TEM, TEM%, and R coefficient for intra- and inter-observer repetitive measures per variable are given in supplementary Table 1.

All variables from lower and upper endplates and from each population plotted against age-at-death and colored by sex (blue: males; red: females) are provided in supplementary Figs. 1-4, while variable statistics (mean, standard deviation) per sex are also available in Table 2. The results for sexual dimorphism are summarized in Table 3. All size-related variables, i.e. perimeters and areas, in both populations were found to be sexually dimorphic with male vertebral endplate variables being larger than their female counterparts. It is remarkable, though, that the lower and upper internal ArPerIndex shape variables also exhibited sexual dimorphism in the Danish populations, although this is not true for the Greek population. The two statistically significant cases of LieAreaRatio (P-value: 0.007; ω²: 0.041) and LieArPerIndex (P-value: 0.017; ω²: 0.030) in the Greek sample should not be considered as strong evidence of sexual dimorphism, since their respective effect size is too small and most likely their statistical significance is sample related rather than population representative.

Table 2 . Descriptive statistics for each variable for lower and upper endplates

Measurement	Greek population					Danish population
	Male μ	Male σ	Female μ	Female σ		Male μ	Male σ	Female μ	Female σ
LiAreaa	591.7	150.7	549.8	155.0		647.7	212.5	646.4	177.44
LiPerb	90.9	11.2	87.1	15.2		94.7	14.6	93.7	17.46
LeAreaa	1204.9	233.1	1088.5	227.1		1296.6	311.9	1263.4	319.84
LePerb	132.1	13.5	124.8	19.2		138.2	17.1	133.9	24.13
LieAreaa	613.1	136.6	538.8	135.1		648.9	161.3	617.0	184.15
LieAreaRatio	0.510	0.068	0.490	0.100		0.506	0.080	0.477	0.0947
LiArPerIndex	0.886	0.034	0.870	0.115		0.883	0.035	0.877	0.1319
LeArPerIndex	0.861	0.033	0.848	0.111		0.843	0.035	0.841	0.1273
LieArPerIndex	0.438	0.054	0.420	0.082		0.426	0.067	0.410	0.0814
UiAreaa	577.5	134.3	537.4	138.3		662.3	180.7	611.2	161.48
UiPerb	90.0	9.9	86.6	14.5		95.9	12.9	91.0	16.91
UeAreaa	1126.2	222.8	1018.4	221.3		1231.8	278.9	1187.0	303.35
UePerb	127.4	13.4	120.6	19.5		135.0	16.3	130.8	24.65
UieAreaa	548.7	136.9	481.0	138.1		569.5	144.0	575.8	178.24
UieAreaRatio	0.486	0.067	0.464	0.096		0.465	0.064	0.471	0.0965
UiArPerIndex	0.884	0.048	0.863	0.118		0.889	0.034	0.875	0.141
UeArPerIndex	0.865	0.037	0.849	0.116		0.841	0.041	0.826	0.1356
UieArPerIndex	0.420	0.055	0.399	0.079		0.391	0.054	0.397	0.0804

μ, mean; σ, standard deviation. ^aMeasured in millimeters, ^bMeasured in squared millimetres.

Table 3 . Expression of sexual dimorphism of the lower and upper endplates of T12

Measurement	Greek population						Danish population
	Male normality	Female normality	F/Chi-sq	P-value	ω2/ε2		Male normality	Female normality	F/Chi-sq	P-value	ω2/ε2
LiArea	0.283	0.406	77.015	<0.001	0.333		0.147	0.712	62.281	<0.001	0.373
LiPer	0.886	0.950	87.425	<0.001	0.362		0.745	0.851	51.632	<0.001	0.330
LeArea	0.190	0.043	58.698	<0.001	0.385		0.523	0.923	120.481	<0.001	0.537
LePer	0.580	0.012a	47.007	<0.001	0.307		0.033a	0.708	55.992	<0.001	0.544
LieArea	0.530	0.005a	21.295	<0.001	0.135		0.304	0.253	61.928	<0.001	0.372
LieAreaRatio	0.623	0.170	7.540	0.007	0.041		0.046a	0.523	0.137	0.712	–0.009
LiArPerIndex	0.002a	0.482	0.197	0.657	–0.005		0.612	0.114	34.596	<0.001	0.246
LeArPerIndex	0.502	<0.001a	0.264	0.607	–0.005		0.046a	0.573	3.012	0.083	0.020
LieArPerIndex	0.871	0.290	5.781	0.017	0.030		0.320	0.291	0.037	0.847	–0.009
UiArea	0.166	0.183	68.868	<0.001	0.314		0.729	0.973	47.926	<0.001	0.345
UiPer	0.458	0.199	77.112	<0.001	0.340		0.653	0.813	37.870	<0.001	0.293
UeArea	0.124	0.231	91.750	<0.001	0.380		0.512	0.625	72.804	<0.001	0.447
UePer	0.006a	0.073	52.428	<0.001	0.352		0.309	0.452	48.727	<0.001	0.349
UieArea	0.007a	0.007a	29.460	<0.001	0.195		0.035a	0.162	28.568	<0.001	0.317
UieAreaRatio	0.711	0.959	0.560	0.456	–0.003		0.230	0.222	0.069	0.793	–0.011
UiArPerIndex	0.137	0.001a	0.069	0.793	–0.006		0.044a	0.350	13.833	<0.001	0.148
UeArPerIndex	0.001a	0.001a	0.012	0.914	–0.007		0.635	0.115	3.021	0.086	0.022
UieArPerIndex	0.624	0.614	0.612	0.435	–0.003		0.391	0.975	0.109	0.742	–0.010

T12, 12th thoracic vertebra. ^aChi-square statistic and ε² effect size correspond to variables with at least one group statistically significant on Saphiro-Wilk’s test.

Univariate regression of age-at-death against the morphometric variables of the endplates on the Greek population revealed that the total area of the endplate, the area of the epiphyseal rim, and the amount of shape irregularities of the endplate’s external boundary are mostly affected by age, with the lower and upper endplates as well as males and females following somewhat different patterns. Our results, summarized in Table 4, showed that the morphology of the upper endplate is more correlated with age than the lower endplate, as expressed by the explained deviance and the F-statistic of the respective GAMs. Regarding sex differences, the total area and the perimeter of the endplate (*eArea, *ePer) as well as the shape of the external boundary (*eArPerIndex) are more prominently affected by age in females, whereas the area of the epiphyseal rim and its ratio over the external area, especially on the upper endplate (UieArea, UieAreaRatio), are heavily affected by age in males. All regression models for predicting variables from age-at-death are plotted in supplementary Figs. 5–8. Regressions of selected variables with statistically significant GAMs are shown on Fig. 3 with male and female groups plotted side-by-side. As shown in figure 3, all size-related variables increase with age including the ratio between the endplate’s total area and the central endplate, whereas the external Area-Perimeter Index decreases, manifesting a more irregular shape of the external boundary at older ages. Most notably, the external area and perimeter of both endplates increase continuously in females, whereas in males they remain constant at younger ages before they start increasing after the age of 40. The opposite was observed in the external Area-Perimeter Index, especially of the upper endplate, which appears relatively unaffected by age in young adult females. According to our results, the most consistent pattern observed in both sexes concerned the area of the epiphyseal rim (*ieArea), which gradually increases with age at a similar rate. On the other hand, the ratio of the epiphyseal rim (*ieAreaRatio) is prominently affected by age only in the upper endplate of males (deviance explained 32.3%), whereas in male lower endplate as well as in females the explained deviance was significantly lower (6.7%–12.8%).

Table 4 . Univariate regression with generalized additive models in the Greek population

Morphometric variable	Predicting variables in males					Predicting variables in females					Predicting age in males					Predicting age in females
	R2 adjusted.	Deviance explained	F statistic	P-value		R2 adjusted.	Deviance explained	F statistic	P-value		R2 adjusted.	Deviance explained	F statistic	P-value		R2 adjusted.	Deviance explained	F statistic	P-value
LiArea	0.0137	3.6	1.8493	0.4250		0.0145	4.0	1.5581	0.4508		–0.0118	0.0	0.0233	0.8797		–0.003	1.2	0.8024	0.3738
LiPer	0.0186	4.1	2.2534	0.3612		0.0081	3.2	1.3757	0.5216		–0.0096	0.2	0.1975	0.6579		–0.0056	0.1	0.6343	0.4289
LeArea	0.1345	15.7	13.5473	0.0044		0.2143	22.6	19.0073	<0.0001		0.0933	10.4	9.6460	0.0026		0.2143	22.6	19.0043	<0.0001
LePer	0.1968	21.6	21.4629	0.0002		0.2980	30.9	29.0392	<0.0001		0.1708	18.1	18.3015	0.0001		0.2980	30.9	29.0173	<0.0001
LieArea	0.2127	22.5	25.5148	<0.0001		0.2647	27.6	24.7600	<0.0001		0.2083	21.8	23.107	<0.0001		0.2755	29.3	24.4858	<0.0001
LieAreaRatio	0.0755	8.7	7.8608	0.0063		0.1069	12.8	7.7663	0.0172		0.0755	8.7	7.8627	0.0063		0.0965	11.0	8.0514	0.0061
LiArPerIndex	–0.0008	1.1	0.9311	0.3374		–0.0058	0.9	0.6197	0.4341		0.0310	6.0	3.7325	0.3350		–0.0058	0.9	0.6197	0.4342
LeArPerIndex	0.2539	26.3	29.5794	<0.0001		0.3127	33.2	31.2599	<0.0001		0.2539	26.3	29.58	<0.0001		0.2940	30.5	28.4789	<0.0001
LieArPerIndex	0.0244	3.6	3.1050	0.0817		0.0533	8.4	3.8480	0.1962		0.0244	3.6	3.1053	0.0817		0.0103	2.5	1.6827	0.1993
UiArea	0.0416	6.0	3.6402	0.1132		0.0857	11.9	7.0206	0.0869		0.0311	4.3	3.5968	0.0616		0.0237	3.9	2.5745	0.1136
UiPer	0.0380	5.8	4.2327	0.1435		0.0859	12.0	6.9904	0.0944		0.0208	3.3	2.7228	0.1030		0.0224	4.0	2.7106	0.1815
UeArea	0.0660	8.6	6.6975	0.0465		0.2706	29.3	24.6259	0.0001		0.0490	6.1	5.1704	0.0257		0.2438	25.5	21.9678	<0.0001
UePer	0.1382	15.7	13.9764	0.0024		0.3377	36.0	33.7354	<0.0001		0.1190	13.0	11.9427	0.0009		0.2975	30.8	28.5353	<0.0001
UieArea	0.2991	30.8	36.2725	<0.0001		0.2124	22.5	18.5347	0.0001		0.2984	30.7	35.4489	<0.0001		0.2124	22.5	18.5325	0.0001
UieAreaRatio	0.3143	32.3	38.1285	<0.0001		0.0528	6.7	4.6243	0.0353		0.3143	32.3	38.1295	<0.0001		0.0529	6.8	4.6335	0.0355
UiArPerIndex	–0.0014	1.1	0.8833	0.3502		–0.0090	0.7	0.4188	0.5199		0.0395	7.3	4.0652	0.3420		–0.0042	1.5	0.3091	0.6526
UeArPerIndex	0.2001	21.0	21.2698	<0.0001		0.2299	25.6	20.5795	0.0005		0.2087	22.4	20.5262	0.0001		0.1957	21.5	17.7321	0.0006
UieArPerIndex	0.1889	19.9	19.8570	<0.0001		0.0207	4.7	2.1276	0.3607		0.1889	19.9	19.8627	<0.0001		0.0000	1.5	1.0012	0.3214

Deviance explained in percentage (%).

Figure 3. Regression of age against selected morphometric variables with statistically significant models.

The results of the univariate regression of the variables for predicting age with GAMs are also provided in Table 4, while the corresponding regression plots are available in supplementary Figs. 9–12. As expected, they are in accordance with their ‘predict variables from age’ counterparts although to some lesser degree as expressed by the lower explained deviance of these models. These dissimilarities are expected since the majority of the produced GAMs are non-linear. It is important to note, though, that with respect to age prediction the analyzed variables in females can account up to about 31% of the explained deviance regarding the external perimeter of both endplates as well as the external Area-Perimeter Index of the lower endplate. In males, on the other hand, the best age-predicting variables concern the area and ratio of the upper endplate’s epiphyseal rim (UieArea: 30.7% deviance explained; UieAreaRatio: 32.3% deviance explained) followed by the external Area-Perimeter Index of the lower and upper endplates with explained deviance of 26.3% and 22.4% respectively.

The combination of variables for multivariate regression has been explored among the variables with statistically significant univariate regression models and the resulted optimal multivariate GAMs are summarized in Table 5. The best performing multivariate GAM concerns the upper endplate of male individuals explaining more than 46% of the observed variation, whereas the same endplate side albeit with different variables in females yielded the poorest performance with only 32.5% explained deviance. On the lower endplate, females’ GAM outperforms males’ GAM with explained deviance of 41.6% and 38.6% respectively. Conforming with the observations of the univariate results, the area of the epiphyseal rim (*ieArea) is included in all four models, the external Area-Perimeter Index is the next most used predictor, the external perimeter plays an important role only in females, whereas the ratio of the epiphyseal rim’s area (UieAreaRatio) is an effective predictor only in males’ upper endplate.

Table 5 . Multivariate Generalized Additive Models and their application on Danish sample

Group	Variable	Greek population			Danish population
		R2 adjusted	Deviance explained		Under-estimated	Within-range	Over-estimated
Males–lower endplate	LeArPerIndex, LieArea	0.3678	38.3		0	18	35
Females–lower endplate	LePer, LeArPerIndex, LieArea	0.3777	41.6		3	15	32
Males–upper endplate	UieAreaRatio, UieArea, UeArPerIndex	0.4296	46.1		1	14	31
Females–upper endplate	UePer, UeArea, UieArea	0.2925	32.5		1	11	31

Deviance explained in percentage (%).

The application of the multivariate GAMs on the archaeological Danish material yielded a correct age prediction rate ranging from 25.6% (female upper endplate) up to 34% (male lower endplate). These results are included in Table 5 as underestimated, within range (correctly predicted), and overestimated individuals for each corresponding multivariate GAM. The majority of the misclassified individuals were overestimated in age compared with the estimated age-at-death ranges, a comparison illustrated in the histograms of Fig. 4. It should be noted that for the correctly estimated (within range) individuals, the colored range maps of the ‘within range’ bar represent the ratio of the distance of predicted age over the mean age range previously estimated osteologically multiplied by a fifth of the estimated age range, whereas the colored range maps of the bars for the under- and over-estimated individuals represent the distance in years between predicted age and the closest margin of the previously estimated age range.

Figure 4. Evaluation of multivariate models on age prediction of the archaeological sample.

The present study introduces a robust and highly reproducible method for quantifying the morphology of the vertebral endplates for exploring correlations with age. Despite certain amount of variation being unaccounted for due to the 2D nature of the captured images and their subsequent post-processing for variable extraction, our method presents significant advantages over previous studies using categorical scoring schemes for vertebral osteophyte formation. Beyond the reliability of the adopted method, as shown by the measurement error results, a major advantage is that the extracted variables are continuous, hence offering a greater flexibility for statistical analysis, especially for exploring correlations with age. This is most prominent with the LeArPerIndex/UeArPerIndex variables that consistently capture and quantify the endplate’s external boundary’s irregularities, which are directly related to the formation of osteophytes. Another advantage is related to the acquired images, which bear the potential for a fully automated variable extraction in future works, although in this study the manual extraction of the polyline coordinates was favored for its simplicity despite its significant labor.

The presence of sexual dimorphism regarding the size-related variables is in accordance with previous findings on the Greek population [29] as well as the well-established fact that males are larger than females, in general. On the other hand, the shape of the endplate and especially its external boundary, which relates to osteophyte formation, does not exhibit differences between males and females. Hence, our results suggest that the processes involved in the formation of osteophytes at the boundaries of the vertebral endplates share the same underlying causes between sexes, such as trauma, congenital diseases, or aging [18]. The only exception is observed in the Danish population and concerns the shape of the boundary of the more porous central endplate (*iArPerIndex), which is sexually dimorphic with considerable estimated effect sizes, although not as much as those observed in the size-related variables such as the area or perimeter of the endplate. Since age-at-death does not have any significant effect on the respective variables (LiArPerIndex, UiArPerIndex) according to our regression results, this exception might be related to differential mechanical loading or growth patterns between males and females of the Danish archaeological population, which in turn may be related to sex division in labor or differences in food consumption and/or health status [30, 31].

The rationale behind using GAMs for investigating the effects of age on the vertebral endplates, including osteophyte formation, is that the age-related processes “driving” the observed changes are not necessarily linear, a fact also corroborated by our results. Furthermore, imposing certain assumptions regarding the expected fitted models would be unreliable, since not much is known on how the age affects the morphology of the vertebral endplates in order to select specific non-linear kernels for our regression models. GAMs, however, allow the computation of the best fit regression model on a dataset solely based on the available data and, most importantly, a mixed model approach via REML for estimating the smoothing parameter provides the appealing balance between flexibility and interpretability [32]. Hence, we argue that GAMs provide the optimal regression analysis tool for modeling the effects of aging on morphometric variables of the vertebral endplates, despite linear regression being regularly employed in relevant works to date [16, 21, 24, 25].

Our results on age-related changes of the vertebral endplates show a positive correlation of osteophyte formation with age-at-death, which is in accordance with previous studies [16, 19]. Interestingly, osteophyte formation in males increases consistently throughout adult life, while in females their development appears to start at a later stage in adult life and increases more rapidly thereafter, an observation more prominent on the upper endplate (Fig. 3). Despite this observation might being sample specific and our results certainly not allowing for any generalization, this differentiation could be attributed to sex differences in labor or habitual activity patterns [33, 34]. Nevertheless, our findings indicate that the effect of aging on osteophyte formation is more prominent in older females rather than males. Accordingly, a more pronounced effect of aging in females is also seen in the size-related variables of external area and perimeter of the endplates, which grow larger as age increases. However, it should be noted that the increase in these variables is most likely related to the formation of osteophytes around the endplate rather than the actual size of the vertebra itself.

The multivariate results confirm a strong correlation between age-at-death and the morphology of the vertebral endplates, as described by the selected morphometric variables and at least for the Greek population, whose sample was employed in this analysis. This correlation is stronger in males than females, which is in accordance with previous works [16, 19]. Most notably, not only the explained deviances of the GAMs differentiate between males and females as well as between upper and lower vertebral endplates, but there are also sex differences in the optimal, i.e. most predictive, variables involved in each GAM. Although there is no apparent explanation as to why females perform better in the lower endplate while males in the upper, it should be noted that osteophyte formation in this work is measured as a 2D planar projection of the endplates outline instead of the typical four-grade scoring system proposed by Watanabe and Terazawa [18]. Hence, our results are not directly comparable to earlier works and the different osteometric approach adopted in this work may well produce different patterns. Moreover, the observed sex differences in the multivariate GAMs could just be specific to the Greek population, in which aging has a more pronounced effect on the lower endplate of females rather than males, while, on the contrary, the upper endplate of males is more affected by age than in females. This distinct pattern is also corroborated by our univariate regression results (Table 3). Nevertheless, the here proposed method for quantifying the endplate morphology, including the degree of osteophyte formation, provides a consistent platform for inferring age-at-death albeit with limited predictive power as demonstrated by the application of multivariate GAMs on the archaeological Danish sample, in which about one third of the sample was accurately predicted within range, while the majority of the remaining sample was overestimated. Taking into account the possible population differences between the referenced modern Greeks and the archaeological medieval Danes and the fact that our results rely on osteological age estimates for the latter, we argue that the proposed method could provide an age-estimation framework for bioarchaeological settings, but only in the concept of broader age ranges, such as discriminating between young and old adults.

Our proposed method for reliably capturing the morphological variation of the vertebral endplates, including osteophytosis, in a continuous manner revealed a significant positive correlation between age-at-death and the morphology of the vertebral endplates of T12 tied to the formation of osteophytes. The observed correlation was stronger in males than females in the Greek population, with the latter developing osteophytes at later stages in adult life, possibly due to sex differences in labor or other habitual activity patterns. The developed multivariate GAMs showed differences in predicting age-at-death between the two sexes, as well as between the upper and lower endplates. However, from a forensic and bioarchaeological perspective, their maximum explained deviance on the Greek population as well as their application on the archaeological Danish sample showed low to moderate capacity in age-at-death estimation. Based on our results, we argue that the proposed method conveys certain merits in age-estimation frameworks, especially in distinguishing between younger and older adults, when more informative age markers are unavailable.

This work was partially supported by an Erasmus+ internship program (7034/2019) issued by the State Scholarships Foundation (Greece) granted to Nefeli Garoufi.

acb-55-4-441-supple.pdf

Conceptualization: AB. Data acquisition: NG, AB. Data analysis or interpretation: NG, AB. Drafting of the manuscript: NG. Critical revision of the manuscript: AB, MLSJ, CV, MEC. Approval of the final version of the manuscript: all authors.

No potential conflict of interest relevant to this article was reported.