The southern Andes of Argentina and Chile is a mountain belt formed at the convergent margin between the Nazca and South American plates. Our study area encompasses Central Argentina and Chile and is located in the transition zone between the Pampean flat slab, characterized by a shallow subduction angle resulting in a volcanic gap, and the normal subduction zone of central Chile and Argentina (Ramos and Folguera, 2009). This segment corresponds to the northernmost active volcanism of the Southern Volcanic Zone of the Andean belt. In this latitudinal block, the Andes are composed of seven major morpho-structural provinces that provide our framework for the sampling and analysis of bioavailable strontium. Since it represents a particular strontium biozone in terms of human diet and mobility (Banner, 2004; Standen et al., 2018), we also include the Pacific Ocean as the first unit along the west to east SAST, as follows: (1) Pacific Ocean; (2) Coastal Cordillera; (3) Central Depression-Central Chilean Valley; (4) Western Principal Cordillera; (5) Eastern Principal Cordillera; (6) Frontal Cordillera; (7) Precordillera; and (8) Active Foreland or Quaternary lowlands (Figure 2).

The first unit is the (1) Pacific Ocean that was sampled in coastal settings between the modern towns of Los Vilos in the north (31°55′S) and San Antonio in the south (33°40′S) (Figure 1). Eastwards from the Pacific coast follows the Coastal Cordillera (2), which has a Paleozoic-Triassic basement in the west and Jurassic-Cretaceous clastic and volcanoclastic intra-arc sequences in the east. These units are intruded by Upper Cretaceous igneous rocks of dioritic to monzodioritic composition. The Central Depression (3) is an extensional basin developed since late Eocene-Oligocene times and filled by volcanic and volcanoclastic rocks (Charrier et al., 2002). During the early to mid-Miocene, calc-alkaline andesitic lava, and acid pyroclastic flows of the Farellones Formation were deposited in the central part of the basin (Kurtz et al., 1997). From the eastern flank of the Central Depression develops the Cordillera Principal, divided into two geological provinces with different rock ages and compositions: the Western Cordillera Principal (4), formed by Oligocene and Miocene volcanic and intrusive rocks, and the Eastern Cordillera Principal (5), constituted by sedimentary Mesozoic rocks deformed into the Aconcagua fold and thrust belt (Giambiagi et al., 2009). The Jurassic to middle Cretaceous sedimentary rocks of the Eastern Cordillera Principal include marine evaporites, sandstones, clay-stones, and limestones. Immediately eastwards lies the Cordillera Frontal (6), containing a Paleozoic basement constituted by sedimentary, metamorphic and igneous rocks intruded by Upper Paleozoic granitoids unconformably covered by Permo-Triassic igneous rocks. The Precordillera unit (7) is emplaced eastwards of the Cordillera Frontal and has a southern limit at ∼33°S. It is composed of Lower and Upper Paleozoic (Cambrian to Permian) sedimentary rocks intruded by Paleozoic calc-alkaline granitoids (Furque and Cuerda, 1979; Astini et al., 1995; Giambiagi et al., 2003). Finally, late Quaternary clastic deposits occur in the Active Foreland (8), to the east of the Frontal Cordillera and Precordillera. This foreland basin results from the progressive erosion of the rocks forming the different morpho-structural Andean units. Synthesizing these trends, rock age increases eastwards from the Pacific coast, hence leading to low ⁸⁷Sr/⁸⁶Sr in the west and increasingly higher values eastwards (Barberena et al., 2017).

This latitudinal band of the Andes presents the highest average altitude separating the Pacific and Atlantic slopes of South America (Clapperton, 1993). The interaction between the Andean topography and the prevailing westerly storm tracks bringing humidity to the continent from the Pacific Ocean produces marked spatial variation in temperature, precipitation, and ecology (Garreaud, 2009; Oyarzabal et al., 2018; Masiokas et al., 2020). These communities are dominated by shrubs in the lowlands and by an increase in the herbaceous stratum up to 3,300 masl (Muñoz-Schick et al., 2000). There is a succession of vegetation communities by altitude (from east to west): open scrubland (1,000–1,500 masl), sub-Andean and Andean scrubland (1,500–2,700 masl), Andean steppe (2,700–3,300 masl), and Andean desert (>3,300 masl). Andean steppe and scrubland associated with the Patagonia phytogeographical province (3,300–1,400 masl) are followed to the east by the dry scrublands of the Monte phytogeographical province (Abraham et al., 2009).

The Andean highlands are characterized by large amounts of precipitation occurring mostly during winter as snow and, as altitude decreases, ecosystems become increasingly dry east and westwards. In this geographical context, the SAST between the Pacific Ocean and the eastern Andean lowlands connects a succession of tightly packed altitudinally arranged ecosystems (Figure 3). This is a highly seasonal landscape where large areas above 2,500 masl would only be available for systematic occupation and circulation during the summer months (Durán V. et al., 2018). This is a fundamental landscape feature that would have conditioned patterns of trans-Andean human movement and social interaction (Cornejo and Sanhueza, 2011; Cortegoso et al., 2016).

The archeological record from the southern Andes provides a window to study diachronic socio-demographic processes spanning the early human dispersion and colonization of the continent beginning in the late Pleistocene, the effective occupation of available spaces within socially defined territories, processes of economic intensification and reduction of mobility, development of agropastoral economies, and socio-political complexity ending with the successive annexation by the Inka and Spanish empires. The scale(s) of human paleomobility and the role of migration through time would have been key components of the trajectories of these and other historical processes.

Humans have occupied the coastal and mountain environments of the southern Andes of Argentina and Chile since the late Pleistocene (García, 2003; Méndez Melgar, 2013; Méndez et al., 2018). For most of the Holocene, mobile foraging was the dominant way of life, documented by a variety of hunting weapon systems, the consumption of diverse plant taxa such as algarrobo (Prosopis sp.) and molle (Schinus polygamus), and the wild camelid guanaco as some of the main staples (Bárcena, 2001; Llano, 2015; Cornejo et al., 2016). There are indications of a major regional shift in human adaptations and economic practices starting before 2,000 years BP. Domestic plants are recorded, including maize (Zea mays), quinoa (Chenopodium quinoa), beans (Phaseolus sp.), and squash (Cucurbita sp.) (Planella et al., 2015; Llano et al., 2017), as well as the earliest pottery (Sanhueza and Falabella, 1999; Marsh, 2017). Burial practices became more elaborate and some of the first formal cemeteries in central-western Argentina are dated to these centuries (Rusconi, 1962; Novellino et al., 2013; Menéndez et al., 2014). The llama domestic camelid (Lama glama) would have been present in areas within this region (Gasco, 2018). For at least a millennium, people had an agropastoral lifeway that would have combined hunting, gathering, farming, and herding in diverse proportions. The Inka Empire arrived in central Argentina and Chile immediately after AD 1400 (Cornejo, 2014; Marsh et al., 2017).

The SAST extends for over 350 km between the Pacific Ocean (Chile) and the eastern lowlands in Mendoza Province (Argentina) (Figure 1). In order to build an isoscape of bioavailable strontium across the southern Andes we sampled small rodents—with restricted home ranges—and plants (Price et al., 2002; Copeland et al., 2016; Snoeck et al., 2018; Scaffidi and Knudson, 2020). Our strategy was initially designed to target rodents as the main proxy for bioavailable strontium (Barberena et al., 2017), and we turned to plant samples for areas where rodent samples were not readily available. Accordingly, the sampling was not designed to formally compare these two substrates. However, since we detected important differences for rodents and plants coming from some of the geological regions, we develop two alternative Random Forest models, respectively, analyzing rodents, on the one hand, and rodents and plants, on the other hand. We then explore the robustness of the resulting models and the main biophysical variables that explain the observed isotopic variation.

The rodent samples include modern and archeological specimens (Barberena et al., 2019, 2020) and the plant samples are all modern. Two strategies for sampling plants were applied: the plant samples from Argentina combine five near-located individuals from the same species in each sampling locality and the samples from Chile correspond to only one plant in each case. It is not expected that this may cause any significant difference in the results obtained. Whenever possible, native taxa were selected (Table 2). Considering the wide altitudinal (∼0–3,400 masl), climatic, and ecological variation in the SAST, three different plant types were sampled: in the Monte of the dry lowlands of eastern Argentina we sampled shrubs such as Larrea, Prosopis, Schinus, and Adesmia, with deep rooting systems reaching the groundwater at beyond 20 m-depth; in the dry highlands of the Patagonian phytogeographical Province we selected grasses (Poaceae) with low-depth root systems; and finally in the more mesic lowland settings west of the Andes we sampled evergreen trees (Quillaja saponaria, Peumus boldus, Eucalyptus robusta), where a large part of the roots providing the soil nutrients are located near the surface (Díaz et al., 1999; Jobbágy et al., 2011). Considering that plant physiology conditions the sources of strontium uptake within the same landscape (Capo et al., 1998; English et al., 2001; Reynolds et al., 2012; Grimstead et al., 2017), we will augment our current sampling in order to account for possible inter-specific differences within geological regions.

The geological provinces described above were used as the units of analysis for the bioavailable strontium sampling. However, the samples are not randomly distributed within these regions since in many cases work focused in areas of archeological interest, such as Los Vilos, Combarbalá, and the Maipo River in Chile and the Uspallata Valley and Laguna del Diamante in Argentina. However, these are representative of the geologic areas from which the samples were gathered. The plant and rodent samples are distributed in the spatial units of analysis as follows: Pacific coast, n = 15; Coastal Cordillera, n = 10; Western Principal Cordillera, n = 10; Eastern Principal Cordillera, n = 10; Frontal Cordillera n = 8; Precordillera, n = 24; and Active Foreland, n = 14.

Strontium isotope analysis was performed in the Department of Geological Sciences, University of Cape Town, following routine chemical and MC-ICP-MS methods (Copeland et al., 2010, 2016; Scott et al., 2020). Plant sample material was cut into <1 cm pieces and subsamples of plants from each sampling locality where placed in a pure quartz silica crucible. The crucibles were placed in a muffle furnace set at 300°C, and the temperature ramped by 100°C/h to 650°C and held at that temperature overnight. Once cooled, the ashed plant material was transferred to a 7 ml Savillex Teflon beaker. For digestion, 2 ml of a 4:1 mixture of 48% HF: 65% HNO₃ was added, the beaker closed and placed overnight on a hotplate set at 140°C. Following two stages of dry down and re-dissolving in 1 ml 65% HNO₃, 1.5 ml 2 M HNO₃ was added and allowed to settle overnight to be ready for strontium separation chemistry. Powdered bone and enamel samples were weighed into 7 ml Savillex Teflon beakers, 2–3 ml of 65% HNO₃ added and the closed beakers placed on a hotplate at 140°C for an hour. Following complete sample dissolution, the beakers were opened, and the samples dried. Once dry, the samples were redissolved in 1.5 ml 2 M HNO₃ to be for strontium separation chemistry.

Strontium separation chemistry for all samples followed the same method (Pin et al., 1994) and used Triskem Sr. Spec resin. The separated strontium fraction for each sample was dried down, dissolved in 2 ml 0.2% HNO₃ and diluted to 200 ppb Sr concentrations for Sr isotope ratio analysis using a Nu Instruments NuPlasma HR MC-ICP-MS. Analyses were referenced to bracketing analyses of NIST SRM987, using a ⁸⁷Sr/⁸⁶Sr reference value of 0.710255. All strontium isotope data are corrected for isobaric rubidium interference at 87 amu using the measured signal for ⁸⁵Rb and the natural ⁸⁵Rb/⁸⁷Rb ratio. Instrumental mass fractionation was corrected using the measured ⁸⁶Sr/⁸⁸Sr ratio, the exponential law, and an accepted ⁸⁶Sr/⁸⁸Sr value of 0.1194. An in-house carbonate reference material processed and measured with the batches of unknown samples in this study yielded results (⁸⁷Sr/⁸⁶Sr 0.708914 ± 0.000045 2σ; n = 12) in agreement with long-term results for this reference material from this facility (⁸⁷Sr/⁸⁶Sr 0.708911 ± 0.000040 2σ; n = 414). Total procedural blanks processed with these samples yielded background Sr levels <250 pg, and therefore negligible.

The Random Forest approach is based in a bootstrapping procedure that generates an ensemble (or forest) of trees resulting in different combinations of the target variable (⁸⁷Sr/⁸⁶Sr) and its biophysical predictors (Breiman, 2001; Bataille et al., 2018). Compared to other classification or regression methods, the average of the Random Forest trees has smaller variance and a higher predictive performance. In addition, compared with other machine learning methods it is less prone to overfitting (Kuhn, 2013). As has been robustly demonstrated (Bataille et al., 2018, 2020), this method allows combining quantitative and categorical variables and has a remarkable versatility to generate ⁸⁷Sr/⁸⁶Sr isoscapes.

In order to determine the most parsimonious subset of variables with the highest predictive power, the Variable Selection using Random Forest (VSURF) package was used (Genuer et al., 2019). We incorporate spatial, geological, and bioclimatic variables in the model (Table 1). We also utilize as a predictor variable the “bedrock model” (Bataille et al., 2018, 2020) that predicts the median, Quartile 1 and Quartile 3 of ⁸⁷Sr/⁸⁶Sr values for bedrock in a global scale. This will expand the recent application developed by Serna et al. (2020) for northern Patagonia. The variables are introduced in a stepwise manner in order to detect the combination of predictors that reduces the mean “Out-of-Bag Error” (see details in Supplementary Data).

The Random Forest approach divides the total isotopic sample into a training phase of the model, encompassing 80% (n = 72) of the results, while the remaining 20% (n = 19) was used in a testing phase to assess the predictive capacity of new data of the model produced (see script in Supplementary Material). The full dataset is first randomized to avoid any structure resulting from record or sampling order. The Random Forest was first passed through a tuning phase in order to find the optimal number of variable combinations with respect to the Out-of-Bag error estimate, by using only the subset of variables automatically selected by VSURF. Afterward, we selected the model showing the smaller Mean Squared Error (MSE) between predicted and observed ⁸⁷Sr/⁸⁶Sr values. The relative weight of each variable in reducing the overall error in the Random Forest model was measured by the “variable importance impurity measure,” graphically represented by the partial dependence of ⁸⁷Sr/⁸⁶Sr values on those variables. The ⁸⁷Sr/⁸⁶Sr values predicted by the Random Forest model in relation to the predictor variables were then used to produce the isoscapes, while the residuals between the predicted and observed values provide the uncertainty associated with the isoscape. In addition, the uncertainty in the prediction for each case was measured using the standard error generated by the data used to generate the model. Finally, we used a Mantel correlogram (Legendre and Fortin, 1989) to assess the existence of remaining spatial structuration in the ⁸⁷Sr/⁸⁶Sr values unaccounted for by the covariables of the model, that would indicate that unknown variables have a weight in structuring the observed variation. To accomplish this, we utilized the residuals of the Random Forest model with 10,000 permutations on two distance-matrices produced on the basis of the residuals and the spatial coordinates on the training model (80% of the sample).

The Random Forest model was adjusted with the software Ranger (Wright and Ziegler, 2017) and tuneRanger (Probst et al., 2019) that implement quantile regression for Random Forest, and with VSURF (Genuer et al., 2019). The graphs of partial dependence were produced with the pdp package (Greenwell, 2017) and the analysis of spatial structure in the residuals was conducted with the packages ncf (Bjornstad and Cai, 2020) and ecodist (Goslee and Urban, 2007). All the analyses were performed in R Core Team (2020). For the geostatistical analysis, we compiled a georeferenced database incorporating the results of the Random Forest model. This information was processed in the platform QGIS 3.14 in order to build the isoscapes of bioavailable strontium using the Inverse Distance Weighing—IDW—function (Kootker et al., 2016; Bataille et al., 2018).

We present new ⁸⁷Sr/⁸⁶Sr results for modern plant samples (n = 26) (Table 2), which are combined with previously published results for rodent samples (n = 65) (Table 3). While this sample can be considered as preliminary given the large size of the study-area, the results offer a first characterization of the geological provinces across the southern Andes.

The ⁸⁷Sr/⁸⁶Sr values are normally distributed (W = 0.98094, p = 0.2039) ranging between 0.70393 and 0.71184 with an average of 0.70723 ± 0.00177. When considering the plant and rodent samples separately, the 26 plant samples have a mean ⁸⁷Sr/⁸⁶Sr value of 0.70689 ± 0.00149, a maximum value of 0.71184 (Frontal Cordillera), and a minimum value of 0.7042 (Western Principal Cordillera). The 65 rodent samples have a mean ⁸⁷Sr/⁸⁶Sr value of 0.70736 ± 0.00187, a maximum value of 0.7109 (Precordillera), and a minimum value of 0.70393 (Western Principal Cordillera).

Particularly for the Precordillera geological province, that includes the key archeological locality of the Uspallata Valley discussed here (Figure 1), we have adequate sample sizes for plant and rodent samples to allow a first comparison. Unexpectedly, the respective isotopic ranges do not overlap. Plants have a mean value of 0.70747 ± 0.00047 and rodents of 0.70964 ± 0.00082. Interestingly, not only rodents but also wild herbivore camelids from Uspallata-Precordillera fail to reproduce the signal for modern plants (Barberena et al., 2019). While this sampling was not designed to compare the performance of different substrates to build strontium isoscapes, these differences are considered large and they would have implications for the reconstruction of human paleomobility and migration. Accordingly, we explore next if there are differences in the main predictor variables of ⁸⁷Sr/⁸⁶Sr values in two alternative Random Forest models: one combining values for plants and rodents (n = 91) and another one including only rodents (n = 65). Based on this comparison we assess the predictive capacity of the two alternative models.

In the “Plants + Rodents” dataset, the selected Random Forest model explains a 76% of the variation in the observed OOB values with an MSE of 0.00000008 (⁸⁷Sr/⁸⁶Sr RMSE = 0.0003), suggesting a good power of the predictive variables utilized and a low error. The variables with the highest predictive power are—in decreasing order of importance—Geological Province, Bio15 (temperature seasonality), and Mean-age, with a subordinate role by six other bioclimatic variables (Figure 4A). In the case of the “Rodents” sample, the selected Random Forest model accounts for 85% of the observed values with an OOB MSE = 0.0000007 (⁸⁷Sr/⁸⁶Sr RMSE = 0.0007). The variables with a high predictive power are only three: Geological Province, Mean-age, and Distance to coast (Figure 4B). In sum, the “Rodents” model has a higher predictive capacity than the “Plants + Rodents” model (85 vs. 76%). Additionally, in the “Rodent” model this variation is accounted for by only three predictor variables, two of which are tightly connected with the geological substrate (Geological Province and Mean-age), while the third (Distance to coast) would be related with the decreasing incidence of marine strontium—0.709202 ± 0.000003—(Kuznetsov et al., 2012) with increasing distance from the coast (Alonzi et al., 2020). Importantly, in the “Plants + Rodents” model several bioclimatic variables contribute to structure variation with a substantial role of temperature seasonality (Bio15).

The analysis of the linear fit between the predicted and observed ⁸⁷Sr/⁸⁶Sr values shows that the plant samples have a larger estimation error than the rodent samples, likely affecting the amount of variation, respectively, explained by the two Random Forest models presented, as well as the number of variables selected by each model (9 in the “Plants + Rodents” model vs. 3 in the “Rodents” model, Figure 4). In the “Plants + Rodents” model the validation of the Random Forest model with the new cases from the test set (20% of the sample) accounts for 72% of the variation (R-squared = 0.72) with a RMSE = 0.0008. The Mantel test produces a value of r = −0.07, p > 0.05, suggesting that there is no significant remaining spatial structure in the residuals, although a low autocorrelation was detected in two distance intervals (Supplementary Figure 1A). While the linear regression between the values included in the training and testing sets for the “Plants + Rodents” samples shows a good overall fit, there is a larger dispersion of the residuals in the area of the most radiogenic values—between 0.709 and 0.710—(Figures 5A,B). In the SAST, these samples are located adjacent or within the Precordillera highlands (including the Uspallata Valley), formed by Lower Paleozoic rocks that are the oldest present (Figure 2). This increase in the magnitude of residual values in contexts with higher ⁸⁷Sr/⁸⁶Sr ratios has been shown to be a global phenomenon (Bataille et al., 2020) and explains the existence of considerably higher residuals in the “Plants + Rodents” model compared to the “Rodents” model (Figures 5B,D).

The regression between the predicted and observed results for the “Rodents” sample shows a high fit between predicted and observed ⁸⁷Sr/⁸⁶Sr values (Figure 5C) and the validation of the Random Forest model with the new cases from the test set (20% of the sample) accounts for 95% of the variation (R-squared = 0.95, RMSE = 0.003) (Figures 5C,D). The Mantel test produces a value of r = 0.005, p > 0.05, suggesting that there is no remaining spatial structure in the residuals (Supplementary Figure 1B). In Supplementary Figure 1 we present the isoscapes with the Standard Error, Quartile 1 (25%), Quartile 3 (median), and Quartile 3 (75%). In the Supplementary Materials file we also present the results for all the covariates analyzed in this model.

The analysis of partial dependence between the predictor variables and ⁸⁷Sr/⁸⁶Sr values for the “Plants + Rodents” model suggests an eastwards exponential increase of the values beyond 150 km from the Pacific coast, reaching the highest values at a distance of 250 km from the coast (Figure 6). This is compatible with the underlying trends in bedrock age from west to east, characterized by the succession of increasingly older rocks: Western Principal Cordillera (Oligocene-Pliocene), Eastern Principal Cordillera (Cretaceous-Jurassic), Frontal Cordillera (Permian-Triassic), and Precordillera (Lower Paleozoic) (Figure 2). Eastwards from Precordillera, there is a decrease in ⁸⁷Sr/⁸⁶Sr values in the Active Foreland, a Quaternary basin that averages sediments from Precordillera and younger formations such as the Cordillera Frontal and Cordillera Principal (Figure 6A). While sampling is small in these Quaternary deposits, it is expected that spatial resolution will be lower compared to the Andes (see below). The most important bioclimatic variable identified is Bio15, a coefficient of Seasonality based on temperature, followed by Tmin07—Minimum July temperature.

The “Rodents” model, on the other hand, is accounted for in a much more simplified fashion by only three predictor variables: Geology, Mean-age, and Distance to coast (Figure 7), indicating that the rodent samples track the geological variation more closely than the “Plants + Rodents” model, where several bioclimatic variables contribute significantly.

We have presented the results for bioavailable strontium across the SAST connecting the Pacific Ocean in Chile and the eastern lowlands in Argentina. The application of the Random Forest approach pioneered by Bataille et al. (2018) has allowed comparing two substrates that can be used to build frames of reference for the interpretation of ⁸⁷Sr/⁸⁶Sr values in archeological remains. The “Rodents” model presents the best fit between predicted and observed values (85%) in association with a small error (RMSE = 0.0007). This large amount of variation can be accounted for with a simpler model where only the predictors Geological province, Mean rock age, and Distance to coast explain the observed variation, whereas in the “Plants + Rodents” model several bioclimatic variables—particularly seasonality of precipitation (Bio15)—contribute to explain variation in addition to geological province and Mean-age (Figures 6, 7). Considering that the plant samples analyzed are modern, an issue that will need to be explored in the future is the possible incidence of modern mining activities in the ⁸⁷Sr/⁸⁶Sr values in plants (Burger and Lichtscheidl, 2019; Thomsen and Andreasen, 2019).

While this project was not formally designed to compare alternative substrates to build strontium isoscapes, we suggest that by its higher fidelity to bedrock geology, rock age, and distance to coast, rodents provide a more accurate proxy for the interpretation of values in human remains in our study area (Price et al., 2002; Bentley, 2006). In Figure 8A we present the strontium isoscape and the associated residual values (Figure 8B) to assess the uncertainty of the isoscape across the SAST.