Southern Ethiopia is a region characterized by a complex and diverse set of rock formation histories (Abbate et al., 2015; Figure 1). The area encompasses part of the East African Rift System (EARS), mainly the Main Ethiopian Rift (MER), with significant varations in the underlying geology, ranging from pre-rifting Precambrian basement rock to the more recently dated Quaternary formations (Abbate et al., 2015). The Precambrian basement is the oldest rock outcrop in the region with lithologies consisting of deformed gneisses, amphibolite, and granulite (Yibas et al., 2002; Emishaw et al., 2017). The East African Orogen with the Arabian-Nubian Shield and the Pan-African Mozambique Belt is represented by the Precambrian of southern Ethiopia. These Precambrian rocks are exposed in the Moyale, Adola and Megado, Yabelo, Negelo, Agere Mariam, and Konso areas dominated by granite-gneiss, ophiolitic fold and thrust-out crops dated between 880 Ma and 526 Ma (Yibas et al., 2002). The rocks are interlayered with Archean to Proterozoic plutonic rocks (Philippon et al., 2014; Abbate et al., 2015). Tholeiitic basalt flows and subordinate tuffs dated between 48 Ma and 30 Ma are exposed above the basement rock in the Gamo and Amaro area (Philippon et al., 2014; Emishaw et al., 2017).

Geomorphological, geochemical and structural investigations have shown the western and eastern MER are dominated by several fault lines. Border faults across the Chew Bahir to Gamo Gofa areas were formed through earlier rifting processes during the Miocene (WoldeGabriel et al., 2000). Volcanism progressively stretched southward from the Afar Depression between Somali and the southern Ethiopian plateau. This continued in correspondence with the formation of the EARS rifting which propagated southward. Rifting in the Broadly Rifted Zones (BRZ) of southern Ethiopia is a northward propagation of the Kenya-Turkana Rift of the southern MER which followed a different trend from the MER (Emishaw et al., 2017). In the northern and central MER, volcanism resumed around 10 Ma with exposed outcrops consisting of late Pliocene volcanic basalts and lacustrine sediments in the areas of Nech Sar, Segen, and Mursi (Teklay et al., 1998; Emishaw et al., 2017).

In the MER, tectonism, volcanism, and rifting have formed diverse sets of lithological features (WoldeGabriel et al., 2000). Inside the rift and along its shoulders, Quaternary formations corresponds to alluvial, fluvial and lake deposits. The MER divides the southern sector into the Chamo and Galana basins to the west and east, respectively, by the Amaro Horst. Outside the Ethiopian rift, Precambrian basement rock and banded gneiss are distributed in the Borena, Yabelo, and Konso areas (WoldeGabriel et al., 2000). Finally, climatic conditions triggered intense weathering of outcropping rock formations and deposits thus promoting the formation of extensive pedological cover.

Southern Ethiopia is home to many palaeoanthropological and archaeological sites dating from the Plio-Pleistocene to the Holocene period, offering immense insights into hominin bio-cultural trajectories. Over the last four decades, research in southern Ethiopia has revealed several fossiliferous and artifact-bearing sites. Early palaeontological and palaeoanthropological expeditions took place in the 1970s across the lower Omo, Kibish and Fejej areas in the southernmost Ethio-Kenya border, yielding rich evidence of human bio-cultural evolution and adaptation dating back to the Plio-Pleistocene (Wood and Leakey, 2011; Shea, 2008; Asfaw et al., 1991). The Konso-Gardula Formation in southern Ethiopia has yielded the earliest Acheulean hand axes dated to the early Pleistocene (Asfaw et al., 1992). The region is also known for yielding one of the earliest records of MSA assemblages from the site of Gademotta in the MER (Sahle et al., 2014).

Owing to the rich MSA and Late Stone Age (LSA) cultural and occupational sequence, the southwestern highlands of Ethiopia were hypothesized to have acted as refugia for hunter-gatherers in the backdrop of the Late Pleistocene climatic oscillation (Brandt et al., 2017, 2012; Ossendorf et al., 2019). Evidence of population contact, long-distance raw material procurement, and trade and exchange of artifacts from the area has been documented using geochemical techniques (Brandt et al., 2017). The region has also been considered a corridor for the dispersal of H. sapiens within and beyond Africa (Brandt et al., 2017). Some of the best-known MSA and LSA sites such as Goda Buticha (Pleurdeau et al., 2014), Mochena Borago (Brandt et al., 2012), and Fincha Habera (Ossendorf et al., 2019), are situated in southern, southeastern, and southwestern Ethiopia.

Recent archaeological studies have also documented early and late Holocene occupations in the region (Assefa et al., 2014). Evidence of microlithic artifacts, plant remains and other evidence documented from sites such as Mochena Borago, Dendi and Nech Sar Park provide insights into the subsistence and technologies of hunter-gatherers (Schepers et al., 2020). Several rock art sites dated to the Holocene period are distributed across the southern Ethiopian highlands and lowlands, providing valuable insights into the symbolic and artistic expressions of past societies (Gallinaro and Zerboni, 2021; Assefa et al., 2014; Solomon, 2020). Furthermore, abundant megalithic monuments, stelae and rock art sites are found in the Gurghe, Sidama, and Gedio areas, mainly dating to the medieval period (Duff et al., 2018; Roger and Jean-Paul, 2020).

In this study, a total of (n = 48) samples including plants (n = 15), water (n = 24), soil (n = 7), and tooth enamel of micromammals (n = 2) were used for building the ⁸⁷Sr/⁸⁶Sr isoscape. Meanwhile, a further 35 archaeological fauna specimens were sampled for Sr isotope analysis (Supplementary Table). Plant, soil, water, and enamel samples were collected from southern Ethiopia during 2022 and 2023 field. Sampling surveys were subdivided based on geological formations key to the area of interest, with locations randomly selected across the different lithologies. Soil samples were collected from 15 to 20 cm below the surface. For all samples, we attempted to avoid areas of anthropogenic disturbance and human settlements. For plant samples, grasses and leaves of trees and bushes were targeted. A Garmin ETREX 22x GPS device was used to record points set to the WGS84 reference system for the geographical coordinates. Photographic documentation was produced for each sample and the surrounding area, with additional information such as the type of geological formation and a description of the sample was also being recorded with an iPad. Plant samples were dried in dry air in the Ethiopian Heritage Authority (EHA) lab.

For water sampling, 250 ml LDPE plastic bottles with narrow necks and screw caps were used. Water samples were collected primarily from lakes located in the MER. Additional Sr isotope data (n = 19) of water ⁸⁷Sr/⁸⁶Sr values (Rango et al., 2010; Markowska et al., 2022) and (n = 2) plant ⁸⁷Sr/⁸⁶Sr values (Ghouri, 2024) were extracted from previous studies.

Stable isotope sampling of fossil tooth enamel powder was conducted at EHA. Preparation was performed following the internal isotope protocols of the Laboratory of Paleoanthropology and Bioarchaeology of the Sapienza University of Rome. Prior to this, each tooth surface was carefully cleaned by gently abrading it with a diamond drill bit (>1.0 mm) to eliminate potential contaminants. Enamel from mammalian tooth was then sampled using a Dremel tool fitted with a diamond-tipped bit (>2.0 mm). Between 7 and 15 mg of enamel powder was extracted for Sr isotope analysis by drilling from the root to the crown of the tooth, ensuring an average dietary and environmental record for the period of enamel formation.

Different sample preparation protocols were involved for soil leachate, plants, water and tooth enamel samples prior to analysis in the column chemistry. Soil samples were ground into powder using a mortar and pestle, weighed, and placed into test tubes. Between 1 and 10 g of samples were weighed and placed in a centrifuge tube. Samples were washed with Milli-Q ultrapure water to remove external dust and placed in crucibles of 50 ml quartz. The samples were subsequently added to 8 mL of 1M acetic acid using a pipette. The samples underwent sonication, followed by 12 h resting period, another round of sonication, centrifugation, and transfer into a clean 15 ml Teflon beaker using a sieve to filter out organic materials. Subsequently, 2 ml of Milli-Q Water was added, followed by further sonication and centrifugation. Finally, the substrate was transferred back into a beaker after another round of sieving. The beaker was then dried on a hotplate at approximately 125–130°C, and 1 ml of 12% hydrogen peroxide (H₂O₂) was added to eliminate organic materials, followed by another drying cycle. The dry sample was transferred onto weighing paper to determine the weight of the fraction leached from the original sample.

Plant samples were prepared by breaking them into small pieces using pliers over a large clean Kimwipe and placed into 15 mL centrifuge tubes. The centrifuge tubes were labeled with sample names and were placed in an oven at approximately 50°C for 12–24 h. The samples were placed in a crucible, sonicated in Milli-Q water, followed by 3 times rinse with Milli-Q water, and dried for more than 12 h. The samples were then ashed in a Thermo Lindberg/blue M^TM Furnace at a different temperature starting at 100°C with gradually increasing temperature until reaching 650°C. The samples were then baked for 2 h to ensure complete combustion. The ashed samples were transferred to 1.5 mL centrifuge tubes and dissolved.

The water sample was filtered and a volume of 10–50 mL of filtered water (0.45 μm filter) was evaporated to dryness in a Teflon beaker at approximately 80°C. The residue was dissolved in 3 M nitric acid (HNO₃) to remove any remaining organic material and salts. The solution was then centrifuged to remove any insoluble material prior to the transferring the supernatant to a clean Teflon vial. Strontium was isolated and purified Sr fraction was then dried, redissolved in a small volume of 0.1 M HNO₃, and prepared for isotopic analysis.

Approximately 5–10 mg of enamel samples were pretreated with 1 mL of 0.1 M acetic acid for 30 min and rinsed multiple times with Milli-Q water to remove residues, carbonates and other contaminants. Samples were then dried again at 50°C and subsequently dissolved in HNO₃.

Samples were processed in the Department of Geological Sciences Isotope Geochemistry Clean Laboratory at the University of North Carolina at Chapel Hill. Approximately 7–10 mg of each sample was dissolved in 550 μL 3.5 M HNO₃. Sr was isolated through ion exchange chromatography on a column with EiChrom Sr-Spec^TM resin using 3.5 M HNO₃ to wash the columns and H₂O to elute Sr. Approximately 20 μL of 0.03 N H₃PO₄ was added to the sample and evaporated to dryness. Sr samples were loaded onto single Re filaments with TaF₅ and analyzed in triple-dynamic multi-collector mode with ⁸⁸Sr = 3 V (resistor of 10¹¹ Ω) on the VG Sector-54 thermal ionization mass spectrometer. All data are normalized to ⁸⁶Sr/⁸⁸Sr = 0.1194 assuming exponential fractionation behavior. Data are reported relative to a value for NBS-987 of ⁸⁷Sr/⁸⁶Sr = 0.710250 ± 0.000020 (2s; n = 10). The internal run precision of all analyses (<0.0018%) is far smaller than the long-term external reproducibility of standards (0.0028%). Therefore, we report all uncertainties using the larger value (± 0.000020, 2s, absolute; Supplementary Information).

The Ordinary Kriging approach was chosen to model a bioavailable isoscape of southern Ethiopia. Ordinary Kriging offers the spatial interpolation approach with statistical autocorrelation between two pairs of data points (in this case environmental ⁸⁷Sr/⁸⁶Sr values) based on a variogram model. This approach has been proven to be the effective way of producing Sr isoscapes as it considers distance and degree of variation among sampled points to estimate unknown areas (Britton et al., 2020; Lugli et al., 2022). As mentioned earlier, for the implementation of Ordinary Kriging, it is necessary to develop a variogram with a specific model, which in our interpolation case is Gaussian, as it appears to best fit the shape of the point distribution. Subsequently, the raster obtained was imported into QGIS 3.34.6 to clip the interpolation through a polygon representing the area of interest in our study. We used QGIS 3.34.6 and SAGA 9.3.1 to produce the bioavailable ⁸⁷Sr/⁸⁶Sr isoscape model in the area of interest for the Ordinary Kriging following the import of the shapefile (shp) containing the sample points in the WGS 84 UTM 37N reference system (EPSG: 32637). To measure and evaluate model performance, we employed validation techniques and calculated statistical metrics including the Mean Square Error (MSE), Root Mean Square Error (RMSE), and R² values. A cell size of 1,000 m was chosen, and the program was instructed to perform cross-validation, specifically Leave-One-Out Cross-Validation tests (LOOCV).

All statistical analysis, tests, and visualization were performed in the R Studio (R Core Team, 2023). Descriptive and inferential statistics were used in this study. Summary statistical, such as the mean, median, standard deviation, interquartile range, and range are provided. Inferential statistics were used to draw conclusions about the characteristics of our sampled datasets and determine whether our data followed any observed discrepancies or relationships that were statistically significant for the fauna data. A Shapiro-Wilk test confirmed the enamel data were not normally distributed while parametric tests were applied for the environmental data which met the assumption of normality. A p > 0.05 does not provide enough evidence to reject the null hypothesis. Additionally, linear regression analysis was performed to evaluate the relationships between ⁸⁷Sr/⁸⁶Sr values and important variables such as elevation. Spearman's correlation coefficient was performed mainly to know the relationship between two variables. The Mann-Whitney U test was used for comparing two variables, while the Wilcoxon Signed-Rank Test was used for paired independent variables. An ANOVA test was used to determine the differences among the means of multiple variables.

The summary results of the bioavailable Sr samples for plants, archaeological fauna, soil leachate, and water are presented in Table 1. The ⁸⁷Sr/⁸⁶Sr isotope values for the environmental samples ranged from 0.70396 to 0.71183 (mean, 0.70673 ± 0.00188). Table 2 presents the values for different geological and lithological groups. Plant values display the largest variations (range = 0.70396–0.71183 ± 0.00254) in the overall range compared to other sample types such as soil and water. This is expected as the sampling covers a wider geological rock outcrops and landscapes. The soil and water ⁸⁷Sr/⁸⁶Sr values also vary significantly with ranges from 0.70368–0.70898 ± 0.00189 to 0.70445–0.70887 ± 0.00119, respectively. The variation for water is interesting despite most water samples being taken from similar sources of water bodies mainly from Rift Valley lakes. Minimal isotopic offsets are observed among different substrates of soil leachate, plants and waters within similar geological formations. In all Sr sample types, specimens collected within small distances and the same geology generally show ⁸⁷Sr/⁸⁶Sr isotope values of the Gotera fauna are properly estimating the local Sr distributions among different lithological units and geological outcrops. The Sr isotope data for the environmental variables are normally distributed (P = 0.00070. P < 0.05).

Analysis based on local geological formations also show large variance between them in terms of bioavailable ⁸⁷Sr/⁸⁶Sr values. The most radiogenic values are displayed from the Precambrian rock unit of the Yabelo and Jinka provinces with an empirical ⁸⁷Sr/⁸⁶Sr ratio >0.707, whereas limited values exhibited across the rift valley lakes and Quaternary deposits and soils. Quaternary sediments and Tertiary volcanic rocks have almost overlapping ⁸⁷Sr/⁸⁶Sr values (Figure 2) being deposits the result of the reworking of volcanoclastic material after erosion. Samples from the rift valley lakes of Langano, Shala, Abijata, Awasa, Abaya, and Chamo have generally higher Sr isotopic ratio. Lake Chew Bahir and Abaya have ⁸⁷Sr/⁸⁶Sr mean value of ~0.706, while Lake Langano, Awasa, Shala, and Zeway have a mean value of 0.707. Despite its geographical proximity to Lake Abaya, Lake Chamo has the lowest (~0.704) ⁸⁷Sr/⁸⁶Sr values from the lake samples. Within different lithological units, recent unconsolidated Quaternary sediments have the lowest ⁸⁷Sr/⁸⁶Sr values followed by basic volcanic rocks, while carbonate rocks of the study area show highest ⁸⁷Sr/⁸⁶Sr values.

A linear regression analysis (Figure 3) of the Sr isotope values by elevation shows that isotopic values increase as elevation increases. This is confirmed by a statistically significant correlation coefficient analysis (P < 0.05, P = 0.024). Overall, ⁸⁷Sr/⁸⁶Sr values show large variabilities across our selected areas of interest in southern Ethiopia. These variabilities are correlated with geological and lithological units whereby younger lithologies display lower ⁸⁷Sr/⁸⁶Sr values compared to older Precambrian units with higher ⁸⁷Sr/⁸⁶Sr values.

The ⁸⁷Sr/⁸⁶Sr range of the predicted values in the study area of southern Ethiopia is < 0.705 to > 0.709 (Figure 4). The Sr isoscape reflects the spatial heterogeneity of our area of interest (Figure 5A), with the lowest values displayed to the direction transecting the Gotera area to the northeast and southwest (0.704–0.705). The bioavailable isoscape shows more radiogenic ⁸⁷Sr/⁸⁶Sr values (>0.708) predicted for the eastern and western uplands of Jinka and Yabello areas. Relatively lower values were associated with Quaternary sediments and Tertiary volcanics of the rift floor around the Arba Minch, Lake Chamo, and Gotera areas. Most of the lake water bodies in the MER such as Lake Langano, Zeway, Abijata, Shala, and Abaya displayed predicted ⁸⁷Sr/⁸⁶Sr values between 0.706 and 0.707 (Figure 4). Lake Chamo area displays the lowest predicted ⁸⁷Sr/⁸⁶Sr values. Areas around Lake Langano with environmental samples from hot springs in the central MER show relatively higher ⁸⁷Sr/⁸⁶Sr values compared to other areas in the rift setting.

The RMSE value of 0.00132 indicates a high model accuracy while the R² of 99.206862 reflects over 99% of the variance in the dataset can be explained by the model, implying an excellent model fit. This study models the prediction standard error of the isoscape model to evaluate the prediction performance (Figure 5). The predicted standard error from the model is relatively low (SE = 0.0000000032) and evenly distributed in the sampled locations. We found that the prediction error (Figure 5A) becomes larger in the areas with fewer or no samples. This is expected because the Ordinary Kriging overlooks the consideration of auxiliary environmental predictors such as local geologies, climates, and lithological factors. Overall, our model demonstrates strong predictive performance with a high R² (99.206862) and low MSE (0.000002) and RMSE (0.00132), representing accurate and precise prediction.

The ⁸⁷Sr/⁸⁶Sr isotope range of archaeological fauna displays a smaller range (range = 0.704565–0.706145) (Table 3) compared to the different environmental archives utilized in the current study. The mean value for the Gotera fauna of all analyzed taxa appears to be identical (mean ~0.705). The ⁸⁷Sr/⁸⁶Sr values of the Alecelaphini (range = 0.70465–0.70587), Antilopini (range = 0.70470–0.70566), Bovidae (range = 0.70479–0.70615), Equidae (range = 0.70471–0.70585), and Suidae (range = 0.70462–0.70592) are not significantly different (Figure 5). While relatively larger variations are observed for Suidae taxa (sd = 0.00072), limited variation is documented for Antilopini (sd = 0.00035) (Table 3). Other individuals including a rodent (0.70485), a tragelaphini (0.70457), and a Viverridae (0.70464) show more or less similar ⁸⁷Sr/⁸⁶Sr ranges (Figure 6) corresponding to the GOT10 site plant and soil ⁸⁷Sr/⁸⁶Sr values.

Bioavailable Sr isoscape is key to being able to properly interpreting ⁸⁷Sr/⁸⁶Sr ratios across different geographic areas and its relevance for past mobility (Maurer et al., 2012; Bataille et al., 2018; Wang et al., 2023). Over the past decades, studies have increasingly employed different modeling approaches to construct a bioavailable ⁸⁷Sr/⁸⁶Sr isoscape and Sr isotope studies to understand the origin and mobility patterns of archaeological fauna and human groups in the past as well as for ecological studies (Price et al., 2002; Alexander Bentley, 2006; Maurer et al., 2012; Crowley et al., 2017; Bataille et al., 2018; Britton et al., 2020; Gigante et al., 2023; Wang et al., 2023). Such studies interpret ⁸⁷Sr/⁸⁶Sr values by comparing them with the predicted bioavailable ⁸⁷Sr/⁸⁶Sr isoscape (Bataille et al., 2018).

This study demonstrates that the bioavailable isoscape produced for southern Ethiopia here appears to correspond closely to the underlying local lithology. Though limited, previous Sr isotope researches in Ethiopia on major lithological units (e.g., Pryor et al., 2020) has revealed ⁸⁷Sr/⁸⁶Sr values across the Miocene/Pliocene and Oligocene basalts (0.703–0.706), Cretaceous terrestrial sandstones (0.7068–0.7078), Jurassic limestones (0-709-0.710), and Precambrian basement geology (0.7098–0.7105). Such measured variations in the ⁸⁷Sr/⁸⁶Sr values across the different topographic features reflect the contribution and interplay between various geological formation histories (Crowley et al., 2017). For instance, volcanic rocks generally display low ⁸⁷Sr/⁸⁶Sr values while plutonic rocks display high ⁸⁷Sr/⁸⁶Sr values. In the present study, the geolithologies of the area, according to the GLiM database (Hartmann and Moosdorf, 2012) covers the Precambrian basement rock outcrops prominently exposed across both sides of the rift valley, while the rift floor is comprised of Quaternary sediments and Tertiary volcanic outcrops. In areas to the west and southeast of Gotera, the Precambrian deposit seems to display similar and higher ⁸⁷Sr/⁸⁶Sr values, whereas the Pleistocene unconsolidated sediments show lower ⁸⁷Sr/⁸⁶Sr values. Both the measured and predicted ⁸⁷Sr/⁸⁶Sr values agree with the ⁸⁷Sr/⁸⁶Sr values documented in previous studies across Ethiopia where plutonic rocks display > 0.709 (Teklay et al., 1998). Similarly, the Miocene-Pliocene volcanic basalts from the MER produced very low ⁸⁷Sr/⁸⁶Sr values of 0.70354–0.70431. This is consistent with our measured ⁸⁷Sr/⁸⁶Sr values and the predicted model for the Tertiary Miocene volcanic-dominated areas of Segen and part of Konso.

The observed minimal offsets in the empirical ⁸⁷Sr/⁸⁶Sr values in some areas, for instance, at Lake Abaya and Lake Chamo can perhaps be explained by variations in sediment load and runoff sources. It has been documented that Lake Abaya is known for its high turbidity due to the large amount of radiogenic sediment flow and riverine Sr flux from its catchment areas (Teffera et al., 2019; Markowska et al., 2022). Similarly, the high radiogenic ⁸⁷Sr/⁸⁶Sr isotopic composition of Lake Chew Bahir, is due to the substantial contribution of the Precambrian basement catchment area (Markowska et al., 2022). The low ⁸⁷Sr/⁸⁶Sr values (~0.704) recorded from Lake Chamo is probably due to the Neogene/Paleogene volcanic contributions and limited sediment flow from the high radiogenic catchment areas. Likewise, the relatively high ⁸⁷Sr/⁸⁶Sr values of hot spring waters documented across the MER are likely due to strong influences from bedrock geology. Boschetti et al. (2018) demonstrate that groundwater passes through mineral rich bedrock substrates tends to display high radiogenic Sr values compared to lake waters formed as a result of surface runoffs.

Our bioavailable isoscape (Figure 7) suggests Ordinary Kriging appears to be effective in the areas of interest as demonstrated by the excellent model performance fit and low prediction standard errors (Figures 5B, C). Furthermore, the prediction model works well with higher precision and accuracy particularly in areas with the availability of dense samples. This shows that Ordinary Kriging approach produces a high-resolution bioavailable isoscape which is useful for modeling variations in ⁸⁷Sr/⁸⁶Sr across the landscape. Furthermore, the spatial variations in the predicted ⁸⁷Sr/⁸⁶Sr isoscape correspond well with the geological map of the study area (Figure 1), an alternative map traditionally utilized for Sr isotope studies. Previous studies (Lugli et al., 2019; Wang et al., 2023) have demonstrated the potential of Ordinary Kriging to generate detailed bioavailable isoscapes by integrating different Sr pools including soil, water, plants, and biological samples. Willmes et al. (2021) has also highlighted that Ordinary Kriging works best for small-scale analysis like the present study where with significant localized variations in ⁸⁷Sr/⁸⁶Sr values.

It is, however, important to note that while our Ordinary Kriging model generally works well to produce ⁸⁷Sr/⁸⁶Sr isotope map, it may overlook integrating Sr isotope influences, such as catchment area contributions and other factors highlighted above. This is reflected in the predicted and observed high ⁸⁷Sr/⁸⁶Sr values, while still geologically situated within the Quaternary sediment, for instance, some lakes and hot springs in the central MER. In this study, plant and soil samples perform well and generally reflect the local geology, while water samples display limited discrepancies which is also captured in the predicted ⁸⁷Sr/⁸⁶Sr isoscape.

An Sr isoscape for Ethiopia is available as part of the global isoscape modeled in Bataille et al. (2020) produced using the random forest (RF) machine learning technique, it has been observed that this isoscape shows a poor modeling performance in the areas of Ethiopia compared to the bioavailable isoscape in the part of Ethiopia studied here based on our new data. Bataille et al. (2020) has shown that the global isoscape expected to display low accuracies in the regions with fewer environmental samples regardless of the RF modeling. This suggests that a large sample size remains a fundamental prerequisite coupled with an appropriate choice of modeling approaches. Recent studies demonstrate the better modeling performance of the RF technique for constructing bioavailable isoscapes with the potential of exploring environmental covariates to predict values for less-sampled areas.

In this study, we explored factors including the role of elevation in influencing the bioavailable ⁸⁷Sr/⁸⁶Sr isotopes. Regression analysis (P < 0.05, P = 0.024) indicates that elevation appears to have a positive influence for the high ⁸⁷Sr/⁸⁶Sr values across the high-altitude areas situated just outside the EARS. This can be expected as younger lithological units are across the rift floor with lower elevation while old bedrock lies outside the rift or across the rift shoulders.

Apart from the good model performance, limited offsets have been observed between the measured and predicted ⁸⁷Sr/⁸⁶Sr ratios. Possible variations would be due to several factors such as atmospheric deposition, pedogenic, and geological factors (Maurer et al., 2012; Crowley et al., 2017). Similarly, given that Ethiopia's position to the coast is greater than 100 km, the influence of sea salt aerosols on the isotopic values of southern Ethiopia is expected to be minimal and, therefore, less of a concern in this study. Copeland et al. (2016) suggests areas with over 100 km distance from the coastal lines are less likely to be influenced by the sea spray effect. Similarly, anthropogenic factors including the use of fertilizer and human settlement may have minimal effects as most sampling areas are farther away from human settlements. Therefore, variations and/or similarities in the predicted isoscape may be explained by several factors including atmospheric dust deposition, lithological variations, climate and precipitation patterns.

This is reflected in the overlapping ⁸⁷Sr/⁸⁶Sr values of the Quaternary and Tertiary geological formations which have been shown in the predicted bioavailable isoscape of southern Ethiopia making our interpretation somewhat challenging particularly in areas with Miocene and Pliocene volcanic rocks. Furthermore, as Ordinary Kriging inherently overlooks the influence of geolithological factors, this may have precluded defining better resolution isoscape in the geological boundaries between two or more units (Zaugg et al., 2023). In this regard, complementing and utilizing alternative modeling approaches such as random forest regression (Bataille et al., 2020; Janzen et al., 2020) may overcome caveats as the technique considers environmental, climatic, lithological and geological factors.

Researchers have employed several approaches to define a local bioavailable baseline (Price et al., 2002; Maurer et al., 2012; Lugli et al., 2022). ⁸⁷Sr/⁸⁶Sr values that fall within the estimated local range are considered “local” while those which fall outside are considered “non-local” or migratory (Alexander Bentley, 2006). We employed a combination of different local vs. non-local provenance assignment approaches. The most common method of estimating the local range for a given archaeological material is to define the ± 2 sd of the mean of local ⁸⁷Sr/⁸⁶Sr ratios (Price et al., 2002; Alexander Bentley, 2006; Lugli et al., 2022).

To define the local baseline, we considered the range of microfauna, local plants, and the soil ⁸⁷Sr/⁸⁶Sr ratio of the Gotera area. We estimated the local range from soil (n = 3), small mammals (n = 2), and plants (n = 3) from the Gotera area. The local ± 2 sd range equivalent to the 95% confidence interval of the Gotera area falls between 0.70386 and 0.70568. This study compared the ⁸⁷Sr/⁸⁶Sr ratios of small mammals from the Gotera archaeological site to the local modern baseline (plants and soil). We found that the ⁸⁷Sr/⁸⁶Sr values are consistent over time. Unsurprisingly, the micro-fauna, plant, and soil values from our case study area show relatively limited baseline ⁸⁷Sr/⁸⁶Sr range implying lower variance in the local range and excellent agreement among the different Sr sources. Figure 7 shows the provenance assignment graph for the local archives, reflecting the fact that most of the faunal values fall within the local range. We also compared the corresponding faunal ⁸⁷Sr/⁸⁶Sr values against the empirical and predicted values for the Gotera area. Accordingly, the ⁸⁷Sr/⁸⁶Sr values of most of the Gotera faunal falls within the local range. Recent studies on the fossils and their modern counterparts has documented that species with savannah and open woodland habitats such as the Gotera fauna tend to have wider home ranges (Copeland et al., 2016; O'Brien et al., 2024).

The ⁸⁷Sr/⁸⁶Sr values of the Gotera fauna predominatly aligns well with the the predicted Sr isoscape within 30 km radii of the Gotera landscape, suggesting faunal exploitation most likely occurred in the vicinity of the site. Nearly all faunal ⁸⁷Sr/⁸⁶Sr values falls within the local range established by local fauna and bioavailable isoscape, suggesting limited mobility with local foraging radii of fauna and reliance on local resources. Alternatively, faunal movement could occur from regions or geographical areas with similar ⁸⁷Sr/⁸⁶Sr values. However, given the absence of high ⁸⁷Sr/⁸⁶Sr values exclusively for the archaeological faunal, which were observed in both measured and predicted bioavailable isoscapes of environmental samples, more plausible scenario is likely to be the limited mobility patterns. ⁸⁷Sr/⁸⁶Sr values similar to the fauna in the predicted isoscape is found mainly across the rift setting of the Arba Minch area. However, movement into the area would likely be hindered by the Konso ranges and uplands with Precambrian rock outcrops which have limited sign of isotopic similarity or mobility into the adjacent Precambrian basement areas exposed mainly in the high-altitude vicinities. Nevertheless, it is crucial to highlight the issue of equifinality which obfuscates the discussion of ⁸⁷Sr/⁸⁶Sr data.

The ⁸⁷Sr/⁸⁶Sr values of small fauna such as viverrids and hares in Gotera show values identical to those of modern plant and soil samples from the Gotera area. One individual Tragelaphus (ID, 532) shows similar Sr isotope ratios to the local plant and soil values. Except for micro-mammals, the relatively higher ⁸⁷Sr/⁸⁶Sr ranges of Alcelaphini, and Bovidae (indet.), Equidae and Suidae display almost similar variances. The sparse and limited distribution (Figure 4; Table 3) of the faunal ⁸⁷Sr/⁸⁶Sr isotope ratios is suggestive of the nature of different home ranges for each species. Higher ranges within the same taxon indicate different home-range behavior and mobility patterns both within and probably outside Gotera. Different faunas occupy varying degrees of home ranges from a couple of square kilo meters to hundreds of kilo meters depending on habitat preference, seasonality, and social structure (O'Brien et al., 2024). At Gotera, these animals could have had similar dietary and drinking sources and likely resided in a similar locality. Based on their ⁸⁷Sr/⁸⁶Sr values, none of the Gotera faunas were accessing the Precambrian basement formations across the uphill surrounding area.

Careful consideration, however, should be given to the small sample size (n=35) compared to the anthropogenically assembled large faunas from the GOT10 site. Similarly, integrating Sr isotopic composition against the bioavailable isoscape can be challenging due to issues of equifinality (Holt et al., 2021). Furthermore, it is worth bearing in mind that bulk samples used in this study have provided time-averaged ⁸⁷Sr/⁸⁶Sr values of faunal mobility history during tooth mineralization. One issue is that animals whose tooth ⁸⁷Sr/⁸⁶Sr matches the “local” Sr signature could be quite migratory, traveling to geologies with much higher and much lower ⁸⁷Sr/⁸⁶Sr values, but when homogenized during isotopic sampling, this will “average” to a value that appears identical to the local range. In this regard, sequentially sampled tooth enamel often provides better resolution, mainly to elucidate the seasonal mobilities of game animals. However, this study intends to test whether the faunas are broadly mobile across major geologies and topographies in search of resources, or whether the Gotera area may have offered adequate resources for faunas and humans to be residentially stable. While acknowledging the issue of equifinality, bulk samples have provided important clues about landscape use and mobility patterns particularly in areas such as the Gotera depression with diverse geological and topographic variabilities across the immediate vicinities.

The ⁸⁷Sr/⁸⁶Sr values for some taxa outside the shaded area in Figure 7, including Alcelaphini, and Equidae, could suggest that these individuals are “non-local”. These faunas may have been hunted from the Tertiary rock outcrops probably within the plain area southwest of Gotera toward Yabelo direction with relatively wider ⁸⁷Sr/⁸⁶Sr ranges of 0.70399 to 0.70820. The escarpments across the adjacent highlands are a less likely origin for the non-local fauna. This shows that either at least some individuals within the identified taxa could be migratory or hunter-gatherers may have traveled to hunt them.

The modern savannah ecosystem in eastern Africa hosts one of the largest seasonal wildlife migrations (O'Brien et al., 2023). Migration and mobility patterns in fossil records analogous to the modern fauna have not yet been properly researched except for an ongoing effort (O'Brien et al., 2023, 2024). Likewise, research on the impacts of faunal mobility behavior and seasonal migrations on hunter-gatherer lifestyle is minimal.

The faunal specimens that were exploited by MSA hunter-gatherers in Gotera did not show evidence of significant mobility, suggesting that faunal specimens were accessible resources throughout the year. We argue that the different taxa would have had access to resources year-round in the study area suggesting that MSA human groups may have had adequate resources across the Gotera landscape, which may explain the persistence of intense human occupation in the area. The availability of local resources mainly water and vegetation may also have controlled faunal home ranges by offering productive semi-arid lake-margin environments. The limited ⁸⁷Sr/⁸⁶Sr variations for faunas as shown in Figure 7 indicates MSA hunter-gatherers would have explored a wider yet still home-range-sized, area of the paleolandscape of the Gotera depression. Our current Sr isotope results reveal that the fauna may have roamed within the Gotera area possibly following the palaeobasin watercourses. The Gotera palaeolake offshore as interpreted by sedimentological analysis (Fusco et al., 2021) may have provided suitable palaeoecological niches for diverse faunal and human groups.

Tracking faunal mobility across nearby Tertiary rock outcrops, such as those in the Konso area, is hampered by the overlapping nature of Sr isotope values between the Quaternary sediments and Tertiary volcanics. The bioavailable isoscape also performs less well in this part of the study area. Based on the bioavailable Sr isoscape, low ⁸⁷Sr/⁸⁶Sr values are predicted in the refit floor from the Arba Minch in a north-to-south direction, transecting Gotera. Thus, it is likely that north-to-south small-scale faunal movement would have been possible following the rift basin within an isotopically homogeneous area.

O'Brien et al.'s (2024) recent Sr and δ¹³C stable isotope investigation on the extant fossils from four archaeological sites in Kenya suggests limited herbivore migration during the Late Pleistocene. Migratory and non-migratory behaviors of animals have been recognized from modern and archaeological fauna through intra-tooth sequential Sr isotope analysis (e.g., O'Brien et al., 2023). Fauna may migrate seasonally for different reasons, including searches for food, water, and mating (O'Brien et al., 2024). Future research on the intra-tooth analysis of ⁸⁷Sr/⁸⁶Sr ratios coupled with other stable isotopes will contribute additional insights into the nature of migration and seasonal site use of Late Pleistocene human groups and faunal communities in the region.

Our Sr isotope record indicates the low residential mobility of hunter-gatherers at Gotera. This can be related to the occurrence of rich material records at the site, indicating high occupational intensity and site-use of the Gotera area. Our interpretation corresponds to the exploitation of over 90% of locally sourced lithic raw materials at the GOT10 site (Fusco et al., 2024), suggesting that hunter-gatherers were predominantly relied on local resources. This study supports the ecological assumption that environmental setting with evenly distributed resources and diverse faunal availability characterize prehistoric hunter-gatherers residential stability and territorial lifestyles (Harpending and Davis, 1977). This study further supports Tryon and Faith's (2016) proposed model that suggests reduced residential mobility during the MSA to LSA transition period due to high population density and intensive exploitation of immediate local landscape.