Age estimation for two Mediterranean populations: rib histomorphometry applied to forensic identification and bone remodelling research

Numerous intrinsic and extrinsic factors influence bone remodelling rates and have shown to affect the accuracy of histological aging methods. The present study investigates the rib cortex from two Mediterranean skeletal collections exploring the development of population-specific standards for histomorphometric age-at-death estimation. Eighty-eight standard ribs from two samples, Cretans and Greek-Cypriots, were processed histologically. Thirteen raw and composite histomorphometric parameters were assessed and observer error tested. The correlation between age and the parameters and the differences between sex and population subsamples were explored through group comparisons and analysis of covariance. General linear models assessed through data fit indicators and cross-validation were generated from the total dataset, and by sex and population subsamples. Most of the histological variables showed a statistically significant correlation with age with some differences observed by sex and by sample. From the twelve models generated, the optimal model for the whole sample included osteon population density (OPD), osteon perimeter, and osteon circularity producing an error of 10.71 years. When sex and samples were separated, the best model selected included OPD and osteon perimeter producing an error of 8.07 years for Greek-Cypriots. This research demonstrates the feasibility of quantitative bone histology to estimate age, obtaining errors rates in accordance with macroscopic ageing techniques. Sex and sample population differences need further investigation and inter-population variation in remodelling rates is suggested. Moreover, this study contributes to the creation of population-specific standards for Cretans and Greek-Cypriots. Supplementary Information The online version contains supplementary material available at 10.1007/s00414-022-02812-2.


Introduction
Micro-anatomical features of bone used in forensic age estimation methods were first developed using long bones such as femur, tibia, and fibula [1]. Nearly 30 years later, non-weight bearing bones such as the rib and clavicle were also used for age-at-death (ADD) estimation [2]. Since then, costal elements have been extensively used for histological age estimation [3,4].
A variety of bone histological features (e.g. secondary osteons number, osteon metrics) applied on age estimation methods have demonstrated different correlation to age [3][4][5]. Moreover, the accuracy of the age estimates obtained have shown to be further dependent on intrinsic and extrinsic factors such as historic versus modern samples, pathology and physical activity [6,7]. Additionally, methodological approaches and observer's experience can also influence the results [8,9].
Controversial conclusions have been drawn regarding inter-population variation age-related changes observed in bone microstructure [3,10,11]. Some studies have demonstrated that sample demographics and intrinsic sample characteristics do not have an impact on the reliability of the methods, reporting similar accuracy rates for age estimation population-specific methods as for non-related population formulae [10]. Other researchers have shown that the application of histological population-specific methods provides more accurate age estimates as existing methods produced age estimates not falling within the reported error rates [11][12][13]. Thus, exploring the development of populationspecific histological age estimation formulae has been the target of many studies [4,14].
This paper is a continuation of previously reported results specific to the validation study performed by testing four existing rib histological methods on Cretan and Greek-Cypriot samples [12]. A systematic underestimation of three of the aging methods with an overall increase in errors as age increases has been demonstrated [2,3,15]. Furthermore, higher accuracy of one of the methods was noticed although just for individuals over 60 years of age [5]. In view of these findings, rib histomorphometric parameters are further investigated in the current study to understand the implications of bone remodelling and age-related changes observed on rib cortical bone on the Mediterranean sample (Cretan and Greek-Cypriots). Accordingly, this paper proposes revised population-specific standard methods to estimate AAD through rib histomorphometry for Modern Mediterranean populations.

Materials
Two Mediterranean populations of known AAD and sex were used for this study (Table 1) [12,16]. Ethical approval for this research was given by the University Hospital of Heraklion (Crete) and the University of Edinburgh (UK) ethics committees. Cretan sample consisted of 41 individuals, 34 individuals from the Cretan Osteological collection [17], and seven individuals collected from routine autopsy examination at the Forensic Medicine Unit of the University of Crete. A total of 41 individuals (23 males and 18 females) with a mean age of 57.48 (SD = 21.17) were included. The Greek-Cypriot sample was collected from individuals from the Cypriot Osteological Collection (Limassol, Republic of Cyprus) [18]. A total of 47 individuals (17 males and 30 females) with a mean age of 62.81 (SD = 14.20) were included. The combination of both samples resulted in a total number of 88 individuals with a mean age of 60.33 (SD = 17.89) ( Table 1). The age distribution of the total sample as well as the sample divided by populations can be seen in Fig. 1.

Methods
The target costal element for this study is the 6th left rib; however, left or right standard ribs (4th-8th) were collected when the 6th rib was not available, as no bias is reported by the histomorphometric examination of other standard ribs [19]. The sampling area selected was the rib midshaft with the inclusion criteria considering no damage of the periosteal surface and no trauma or evident pathology affecting bone metabolism.
Thin-sections were processed following standard procedures for bone histological analysis [20,21]. A Leica DM 750P research microscope fitted with a Leica MC 170 HD camera (Leica Application Suite V4 software) was used for capturing 4 × and 10 × magnification microphotographs employed for data collection. Table 2 includes the definition and related data acquisition procedure for all raw and composite histomorphometric parameters explored in this study. Cross-section area measurements were taken on stitched microphotographs, and osteon counting was performed using a standard research light microscope (with polarised filter). The microphotographs were also used to record a count of micro-structures. Single osteon metric parameters were collected through 10 × magnification semi-polarised microphotographs using the segmented tool available in Image J 1.48 software platform [22]. Examples of different age cohort individuals with examples of the parameters collected are presented in Fig. 2. Intra-and inter-observer errors have been already reported for the present histological data [12]. Considering both Technical error of measurement (TEM) and Intraclass Correlation Coefficient (ICC) results, it was demonstrated that overall within and between observer errors fell within the acceptable levels of error, except for N.On. Fg, osteon population density(F) (OPD(F)) and On.Cr that showed that 15-30% of the variance can be attributed to measurement error and moderate to good agreement was reported.
Osteon perimeter observer errors were not tested previously, and the results will be reported here. A total of 22 thin sections were randomly selected from the total sample with intra-observer error being performed on the slides after 3 months from first data collection, while inter-observer was performed by one of the co-authors who has basic training in histomorphometric analysis. Observer error was assessed through TEM reporting TEM values, relative TEM (rTEM), and R as well as through ICC [24,25]. ICC was performed using a two-way mixed effects model under absolute agreement with over 0.80 and 0.90 ICC and 95% confidence intervals considered as good and excellent agreement, respectively [26]. The histomorphometric data was assessed and analyzed through different statistical tests. First, descriptive statistics were calculated to present simple descriptors. Normal distribution of the data was explored through histograms, Q-Q plots, and Shapiro-Wilk test [27]. Depending on the distribution of the data, the relationship between the parameters and age for the total sample was assessed using Pearson's correlation or Spearman's Rank correlation coefficients [28]. The age distribution between sex groups and populations was compared using Welch's t-test [29]. Depending on whether the assumption of normality was confirmed or not, Welch´s t-test or Mann-Whitney U were performed to compare the differences on the histomorphometric parameters and age in the two groups, sex and population samples, separately. One-way analysis of covariance (ANCOVA) was performed to further explore the relationship between age, samples (sex and population) and the parameters. The dependent variable was each histological parameter, and group was set as the factor (either sex or population sample). This analysis allowed to determine whether group membership had a significant effect on the dependent variables and whether age was a significant covariate [30]. If the variables were not normally distributed, the natural log transformation was performed. Finally, general linear regression models (GLM) were generated through linear and multiple regression analyses to produce the age estimation equations for the total sample and sub-samples (sex and population, separately). Residuals were inspected for the assumptions of linearity, normality, independence, and homogeneity of variance [31,32]. For multiple linear regression analysis, independence of the parameters was examined through bivariate correlation, tolerance, and the variance inflation factor statistics (VIF). The GLM models were performed using Gaussian link function and maximum likelihood fitting assessing model selection according to significance levels, ANOVA x 2 test as well as Akaike information criterion (AIC) and Bayesian information criterion (BIC) [33]. Moreover, goodness of fit indicators such as R 2 and standard error of the estimate (SEE) were considered as well as the number, repeatability, and reproducibility of the parameters for the final selection of the optimal models. Leave-one-out (LOO) cross-validation was performed for all regression models to avoid splitting the sample (https:// github. com/ topepo/ caret). Cross-validated results were compared to the initial regression results by means of adjusted R 2 , cross-validated SEE, and mean absolute error (MAE). Statistical analysis was performed using IBM SPSS 27 and R version 4.1.0.

Intra-and inter-observer error-On.Pm
On.Pm intra-and inter-observer error values as obtained through TEM analysis fell within the limits of agreement with rTEM and R reporting values under 5 and above 0.95, respectively. ICC analysis reports excellent agreement for both intra-and inter-observer error with ICC estimates over 0.98 and 95% confidence intervals falling between 0.96 and 0.99 indicating excellent agreement. Table 3 presents data for the total sample including descriptive statistics for the raw and composite histomorphometric parameters. Pearson's correlation coefficient or Spearman's rank correlation coefficient for each parameter are reported (see supplementary material for sex and population subgroups). Among all the variables, the highest correlation with age was reported for OPD(F), OPD, On.Pm, and On.Cr. Figure 3 presents information about the relationship between age and the parameters. Density plots show the frequency distribution of each parameter. The relationship between the parameters is presented using Pearson's

Differences between population samples and sexes
The histomorphometric parameters were explored using Pearson´s correlation coefficients or Spearman's Rank correlation for each subgroup (sex and population) to assess the strength and direction of the relationship between known age and each variable (Table 4). A similar pattern in the direction of the association (positive or negative) as seen for the   Sex and population sample effect was further assessed using one-way ANCOVA. For all the dependent variables tested, known age was reported to be a statistically significant covariate. Only those parameters that do not violate the assumption of homogeneity of regression slopes are presented here. Regarding sex differences and the parameters explored, the only variable indicating a statistically significant sex effect was Ct.Ar (LnCort. Ar: F(1, 85)

GLM: population-specific standards for age estimation
Simple and multiple GLM were generated using the entire sample through the inclusion of those variables that demonstrated a statistically significant relationship with age (Table 3). Only those models meeting the assumptions for linear regression, producing a SEE lower than 15 years, the lowest AIC and BIC values as well as the optimal LOO results are reported. Other considerations related to the general applicability of the models such as repeatability and reproducibility of the histomorphometric parameters were also taken into account. Therefore, only those GLM meeting the above criteria will be fully presented and described (Table 5). For sex and population sample-specific models, only those models showing an improvement in accuracy and fitness indicators as compared to the total sample models will be provided. The total sample dataset (N = 88) was used to conduct simple linear regression analysis through the inclusion of single raw or composite parameters producing models providing a SEE ranging from a maximum of 17.49 years for N.On and a minimum of 11.32 for OPD(F). Although the latter model provided the lowest SEE after cross-validation, this model was not selected as optimal due to the interobserver error reported for this specific parameter somewhere else [9,12,30]. The analysis results in M1 as the optimal univariate model, which includes OPD as single parameter producing a SEE of 12.78 years and cross-validated SEE of 12.88 years (Adj. R 2 = 0.48). For the application of the model for age estimation on an unknown individual, the practitioner needs to apply the values reported in Table 5. For example, the value collected for OPD needs to be multiplied by its unstandardised coefficient (2.893) and the constant added (15.660). The resulting value will be the age estimate (age estimate = (2.893*OPD value) + 15.660). The next best model (M2) includes On.Pm producing a slightly higher SEE and cross-validated SEE (13.36 and 13.53, respectively) ( Table 5). Multiple linear regression analysis was performed exploring all possible combination of parameters having into account the considerations presented above. Three models are provided which include a combination of OPD and osteon metrics parameters (M3 and M4 Total) as well as a model that does not require the inclusion of OPD, but instead selects On.Cr and Ct.Ar (M5 Total). Among these three models, the lowest SEE (10.71 years) and cross-validated SEE (11.01 years) are reported for M4 Total (Adj. R 2 = 0.64).
Sex-specific models were generated by dividing the dataset into males and females, separately, and running simple and multiple linear regression analysis ( None of the population-specific models for Cretans showed an improvement when compared to the total sample  Figure 4 presents the diagnostic plots for the best models for the total sample, and for the sex and population-specific subsamples.

Discussion
It is known that the relationship between micromorphological features and age is not fully consistent, following the same pattern as other biological age markers [34]. Thus, intra-and inter-variability due to sex, pathology, nutrition, physical activity, and genetics -among other factors -may consequently affect not only the application of the methods but also the expression of the microstructural features [35]. Information regarding life style or pathology is not available for the sample under study, and it will be recommended for future research performed on Mediterranean populations to reach a full understanding of remodelling rates in relation to biomechanics and the impact of specific disorders [36,37]. There is still an open debate regarding sex differences in histomorphometric variables, with some studies reporting sex differences in osteonal variables [6,38] while others did not report any [1,39]. Whether the earlier completion of the cortex in females would have an impact on the histological variables and on the age estimates due to differences in the mean tissue age needs further investigation [35]. In the early years of adulthood, males present around 40% larger bone area than females. This observation is associated with the larger body size of males. Additionally, a decrease in cortical area produced at the cross-section periosteal and endosteal regions observed for both sexes around middle age would potentially contribute to these differences [40]. As humans age, males and females experience a general decline in bone formation at the periosteal surface with a reduction in bone formation and continued resorption at the basic multicellular unit level of bone, resulting in a negative remodelling balance [41]. This imbalance is triggered by an increase in cortical remodelling which accelerates cortical and trabecular thinning with interstitial bone becoming highly mineralised in postmenopausal women [41]. Considering that 73% of the females in this sample are over 50 years old, the differences observed on the histomorphometric parameters may result from postmenopausal associated osteoporosis and related bone loss. The sex differences reported here are in agreement with other studies [38,42,43], although caution with microstructures, definition and skeletal element used by the other authors must be considered. In our results, Ct.Ar showed statistically significant differences with age as a covariate, and Tt.Ar and Es.Ar reported differences when the two sexes were directly compared, in agreement with other studies [42,43]. Perhaps, the effect of periosteal apposition seen in postmenopausal women as an adaptive response to the decrease of cortical bone strength and to the increase in bone fragility is responsible for the observed sex differences. However, this can be simply a reflection of sexual dimorphism in rib size as other studies stated substantial sexual differences for Cretans and Greek-Cypriots [17,18]. The key parameter used for age estimation, OPD, did not show any statistically significant differences between the sexes in disagreement with other studies [42]. Further research should be performed on a larger sample with a higher representation of individuals from different age ranges, overcoming the limitations of the current sample and further confirming our results.
Differences between the two population samples, Cretans and Greek-Cypriots, should be minimal. Both share dietary and cultural habits, along with similar climate [18]. However, the populations overall health and habits could have been impacted by differences in historical events. For example, the invasion of Cyprus in 1975 increased the poverty levels for Greek-Cypriots [44] and both islands experienced shifts in economic activities [45], (https:// doi. org/ 10. 1080/ 13683 50030 86679 43). Frequency number of osteons and OPD(F) showed statistically significant differences with higher values for Greek-Cypriots in comparison to Cretans. The same pattern is seen for Ct.Ar, On.Ar, and On.Pm with age as a covariate. This could mean that the packing factor effect relating osteons number, osteon size, and cortex dimensions may be reflected in the trend observed for OPD values between the samples (mean OPDs: Cretans = 14.50, Greek-Cypriots = 16.30) [46]. Cortical bone phenotypic traits -e.g. cross-sectional area and density -are determined in the early years of life. However, during growth and development, these traits are also modulated by environmental factors, with disease and life-style influencing the amount of bone mass reached later in life [47]. Moreover, polygenetic interactions and aging effects on the molecular and cellular processes, decrease of muscle mass, and life style -among other factors -are responsible of inter-individual and inter-population bone modelling and remodelling variability [48]. This topic deserves further investigation, and the inclusion of other parameters such as bone mineral density or other bone surfaces such as trabecular tissue might assist in future research related to the populations under study.
The ultimate goal of this research was the generation of rib histological age estimation standards, as it was demonstrated that existing methods presented several drawbacks for the sample under study [12]. Twelve GLM were reported as the optimal ones considering prediction accuracy, goodness of fit indicators, and cross-validation (Table 5), as well as observer errors [9,30]. If the sex of the individual is known, the use of sex-specific prediction equations is recommended, since overall, they provide slightly more accurate results than the general equations. Regarding the application of the models and considering the potential fragmentation of the remains, sex might be unknown implying that equations for the total sample should be applied. Within those, M3 to M5 provided similar accuracy rates, with the optimal model being M4 included OPD, On.Pm, and On.Cr. Although all these variables reported overall acceptable observer errors [5,12,49], the practitioner must have a suitable training on histomorphometric assessment to ensure reliable and accurate results, especially when OPD is assessed. Metric measurements on secondary osteons might be more feasible for experts with basic training in histology, and thus, M2 and M5 including just On.Pm and On.Pm and Ct.Ar, respectively, might be a better option. Note that On.Cr observer error agreement ranged from poor to good in several studies [12,13], and standardisation for this parameter might need further attention in the future [49]. Regarding sample-specific models, the total GLM should be applied if the unknown individual is suspected to be from Crete as the Cretan dataset did not produce lower errors than the total sample dataset. Now, for individuals of Greek-Cypriot origin, the SEE ranged from 8 to 9.75 years accounting for the lowest errors reported for all the tested models. The same principle applies here for the use of M2 and M4 to avoid errors in the assessment of secondary osteon densities recorded by inexperienced anthropologist, as they apply histological methods for assessing AAD on human skeletal remains.
In general, most of the predicting models included osteon density and/or osteon metrics as age predictors confirming their use for histological aging formulae. Even if OPD is still one of the best predictors, the inclusion of osteon size and shape descriptors has been recommended by other authors; osteon metric parameters are not directly subject to the asymptotic effect, and thus, their use compensates the inconsistencies of OPD for advanced age samples [5]. It is expected that the older the sample, the higher the error rates. This is due to a phenomenon known as the "trajectory effect" which implies that biomechanical and physiological changes in the aging indicators undergo more alterations as the later years in the life span are approaching [50]. Although individuals from younger cohorts should be included in future research in order to fully represent the general population, life expectancy in these countries is increasing steadily, and thus, development of anthropological methods for the elderly are currently required and methods validation needed [51,52]. Our method falls within the expectations reported for histological methods, demonstrating their suitability for age estimation in adults from modern human remains.

Conclusions
The present research is a continuation of the validation study conducted on the Mediterranean sample under study that highlighted the potential need to explore remodelling rates and histomorphometric data on Cretans and Greek-Cypriots. Our results demonstrate that several parameters differ between the sexes and the populations, possibly accounting for variation related to aging as a natural process and age-related changes such as hormonal alterations. Other considerations such as life history and clinical data would be needed to provide further conclusions on the sample under study. The proposed aging equations can be applied on unknown individuals in future forensic cases in Crete and in the Republic of Cyprus. Especially for Greek-Cypriots, this method could potentially assist in the identification of more than a thousand individuals that have yet to be identified, according to the Committee of Missing Persons (www. cmp-cyprus. org). Further research will focus on examining a larger sample size including a detailed representation of adult age cohorts, as well as validating the presented age estimation prediction equations on individuals or collections of the same populations.