A random forest isoscape model of bioavailable Sr for South America: a focus on southern Brazil

In recent years, advances in machine learning have greatly improved the generation of maps showing the geographic distribution of isotope ratios (isoscapes), which have become essential tools for environmental, mobility and provenance studies in both modern and archaeological contexts. Among the various isotopic systems employed, strontium (Sr) is particularly useful because its ⁸⁷Sr $/$ ⁸⁶Sr ratio in the environment is largely controlled by the underlying geology through the composition of local soils and rocks.

In this work, we present a new dataset of bioavailable ⁸⁷Sr $/$ ⁸⁶Sr ratios derived from n=233 plant samples collected across southern Brazil, covering the states of Santa Catarina and Rio Grande do Sul (ca. 370 000 km²). The measured ratios span from 0.70521 to 0.76039 and capture the bioavailable Sr isotope signatures over all major geological units in the region.

We combined these new data with an extensive compilation of published bioavailable Sr measurements from across South America (including plants, fauna, ancient human remains, shells, snails, lichens, water and soils) to construct three random forest Sr isoscapes using different subsets of the combined dataset at the regional and continental scales. The first model incorporates the entire dataset (“All” dataset, n= 883 sites), the second is based on plant + fauna + lichen + human (n= 661 sites) and the third is limited to plant + lichen samples (n= 531 sites). Among the three models, the full dataset model shows lower predictive power, while the plant + fauna + lichen + human and the plant + lichen models yield better results, with similar RMSE (0.0049 and 0.0054) and R² values (ca. 0.76). Compared to existing Sr isoscapes of South America, our models significantly enhance both spatial coverage and resolution of bioavailable Sr predictions, particularly in southern Brazil.

The new bioavailable Sr isotope dataset from Santa Catarina and Rio Grande do Sul states is available at https://doi.org/10.5281/zenodo.17988601 (Scaggion et al., 2025a) and the compiled literature dataset is reported as supplementary material.

1 Introduction

Maps of isotopic variability through the landscape are powerful tools used for answering questions related to migration, trade, diet, provenance, forensic investigations, and environmental changes in both modern and ancient contexts (e.g. Hobson, 1999; Benson et al., 2006; Makarewicz and Sealy, 2015; Bataille et al., 2018, 2020; Hoogerwerff et al., 2019; Lugli et al., 2022; Gigante et al., 2023; Armaroli et al., 2024, Asrat et al., 2025; Dosseto et al., 2025; Scaggion et al., 2025b).

Among the various isotopic systems, strontium (Sr) isotopes play a central role in geotracing and provenance studies of biological materials because biota obtains Sr primarily from soils derived through weathering and disaggregation of local rocks. As a result, the so-called bioavailable ⁸⁷Sr $/$ ⁸⁶Sr act as a proxy for the underlying geology of a region.

Radiogenic ⁸⁷Sr originates by β-decay of rubidium-87 (⁸⁷Rb), while ⁸⁶Sr is stable. Consequently, the ⁸⁷Sr $/$ ⁸⁶Sr ratio in rocks and minerals changes depending on their ⁸⁷Rb content, their ages and geological history providing characteristic isotopic signatures for different geological domains (Faure and Mensing, 2005; Dickin, 2018). Through alteration and pedogenetic processes, Sr ions are leached from rocks into soils and waters, where they mix with other Sr pools derived from rainfall, snow, groundwater, atmospheric deposition or anthropic activities (Vitousek et al., 1999; Voerkelius et al., 2010). These processes can shift the local Sr isotopic signature away from the original bedrock. From soil and water, Sr enters the local ecosystem (Bullen et al., 1996; Capo et al., 1998). This bioavailable Sr is passively absorbed by plant roots and transported to leaves, where Sr²⁺ substitutes for calcium²⁺ (Ca) in metabolic processes (Poszwa et al., 2000; White and Broadley, 2003). From plants and water, the bioavailable Sr enters the local trophic chain. In vertebrates, Sr replaces Ca in the crystalline structure of carbonated-hydroxyapatite (Pors Nielsen, 2004; Bentley, 2006), the main inorganic constituent of bones and teeth (Weiner and Wagner, 1998; Hughes and Rakovan, 2002). Thus, the ⁸⁷Sr $/$ ⁸⁶Sr ratios of bone and dental tissues reflect those of water, plants and animals consumed, which, in turn, mirror the bioavailable isotopic signature of the environment (Schwarcz et al., 2010; Crowley et al., 2017).

Although isotopic fractionation can occur along the trophic chain and through metabolic processes, its enrichment is generally minimal (Lahtinen et al., 2020) and can be effectively corrected using standard normalization protocols during mass spectrometry analysis (Fietzke and Eisenhaure, 2006; Frei and Frei, 2011; Spies et al., 2025). For this reason, the isotopic composition of low-mobility fauna and plants represents a valuable, though not exclusive, proxy for estimating the local bioavailable Sr isotope signatures (Spies et al., 2025). However, deciphering environmental heterogeneity and trophic complexity is not always straightforward; therefore, integrating multiple isotopic and environmental datasets is crucial for constructing robust isoscapes.

Over the last decade, advances in spatial modelling and machine learning techniques such as kriging and random forest algorithms (e.g. Bataille et al., 2018, 2020; Hoogerwerff et al., 2019; Gigante et al., 2023; Armaroli et al., 2024; Scaggion et al., 2025b; Spies et al., 2025), have significantly increased the predictive potential of isoscapes models, by linking spatial isotopic variability to environmental predictors (Bataille et al., 2020). One of the main challenges, however, remains the acquisition of spatially representative and geologically diverse samples, making well designed sampling strategies crucial for developing accurate predictive models.

Brazil represents a particularly interesting case study for Sr isoscape research, due to its geological diversity and rich archaeological and ecological history. In southern Brazil, the states of the Santa Catarina (SC) and Rio Grande do Sul (RGS) encompass various geological formations, from ancient Precambrian crystalline rocks to younger sedimentary basins, which provide a good opportunity to test the spatial resolution, robustness, and applicability of Sr isoscapes. Despite the relevance of Brazil for provenance studies in archaeology, ecology, and food authentication, to date a comprehensive isoscape of bioavailable Sr has not been produced.

In this study, we present a newly generated dataset of n= 233 plant samples, collected across SC and RGS, representing all major geological units of the region (ca. 370 000 km²). This dataset was integrated with a literature-based compilation of bioavailable Sr data from across South America to develop new isoscapes of bioavailable strontium (⁸⁷Sr $/$ ⁸⁶Sr) for South America and southern Brazil, with a focus on SC and RGS. Using a random forest machine learning model following Bataille et al. (2020), we incorporated global environmental variables described in Bataille et al. (2018, 2020) and Reich et al. (2024) to estimate spatial patterns of bioavailable ⁸⁷Sr $/$ ⁸⁶Sr ratios.

The sampled region represents a pivotal area for understanding the expansion of pre-Columbian populations, constituting pass-through areas that facilitated the spread of natives such as the Guaraní, who from approximately the 5th c. CE gradually occupied the continent from the Amazon to the Paraná Delta (Loponte et al., 2024) as well as the movements of Southern Brazilian Highland populations toward the coast. The geological heterogeneity of these states provides a unique natural context for investigating large-scale dynamics of mobility, settlement, and cultural interaction.

2 Geology

We focused our sampling on the southern portion of Brazil, targeting the main geological units in detail to increase the reliability and spatial resolution of the resulting Sr database and isoscape.

The study area, comprising the states of RGS and SC, represents a key segment of the South American continent. The region preserves a long and complex geological evolution, spanning from the Archean to the Cenozoic, encompassing processes of crustal growth, continental collision, sedimentary basin development, and large igneous province magmatism. A simplified overview of this geological complexity is provided in Fig. 1, modified from the 1:5 000 000 geological map of Brazil (Medeiros et al., 2025) downloaded from the Serviço Geológico do Brasil (https://geosgb.sgb.gov.br/geosgb/downloads.html, last access: 3 September 2025) for use in QGis software.

Figure 1 Simplified geological map of southern Brazil, modified from Medeiros et al. (2025) with sample locations collected in this study. The lithology codes are simplified from Medeiros et al. (2025) and grouped according to the two main geological domains of the area: the Mantiqueira Province and the Paraná Basin. The supergroups of the Paraná Basin are subdivided following Milani et al. (2007).

Two major geological provinces dominate this area: the southern Brazilian Shield, representing the Precambrian crystalline basement, and the Paraná Basin, a large Paleozoic-Mesozoic sedimentary and volcanic basin. In addition to these two major geological provinces, Pleistocene and Holocene sediments outcrop along the coast and the southern lagoons of RGS.

The southern Brazilian Shield, also known as the Mantiqueira Province (Heilbron et al., 2004) is composed of a mosaic of cratonic nuclei, accreted terranes, and mobile belts that were amalgamated during the Neoproterozoic, culminating with the Gondwana accretion (Brito Neves et al., 2014). In the state of RGS, the principal tectonic domains of the Mantiqueira Province are the Dom Feliciano Belt and the São Gabriel Terrane (Brito Neves et al., 2014; Philipp et al., 2018). The Dom Feliciano Belt extends in a NE-SW orientation and consists of high-grade metamorphic complexes, syn- to post-tectonic granitic intrusions (650–540 Ma), and metavolcanic-sedimentary sequences (2.35–0.78 Ga). It is interpreted as a Neoproterozoic collisional orogen, developed during the closure of oceanic domains between the Rio de la Plata Craton and the Kalahari-Congo Craton (Brito Neves et al., 2014; Philipp et al., 2018). The belt includes the voluminous granites of the Pelotas Batholith dated between 650 and 580 Ma. The São Gabriel Terrane, located west of the Dom Feliciano Belt, comprises volcanic-arc sequences, ophiolitic remnants, and accretionary complexes. Geochemical data indicate that it represents an intra-oceanic arc formed between 900 and 700 Ma, subsequently accreted to the continental margin during the Brasiliano Orogeny (Cerva-Alves et al., 2021). The Dom Feliciano Belt is also exposed along the coast in the eastern part of SC where it includes the Florianópolis granitic batholiths and the metavolcanic-sedimentary sequences (Brito Neves et al., 2014; Philipp et al., 2018).

The Paraná Basin is one of the world's largest intracratonic sedimentary basins, covering approximately 1.4 ×10⁶ km² across Brazil, Paraguay, Uruguay, and Argentina (Holz et al., 2010; Scherer et al., 2023). It developed over the consolidated Precambrian basement during a long-lived post-orogenic to extensional tectonic regime that lasted from the Ordovician to the Cretaceous.

Milani et al. (2007) divided the Paraná Basin succession into six different second-order supersequences (from oldest to youngest): Rio Ivaí, Paraná, Gondwana I, Gondwana II, Gondwana III (which includes the Serra Geral) and Bauru. In the studied area, the exposed units comprise: (1) Carboniferous-Permian volcanic-sedimentary units of the Gondwana I Supersequence (Itararé group and the Rio Bonito, Palermo and Irati formations); (2) Triassic sedimentary rocks of the Santa Maria Formation belonging to the Gondwana II Supersequence; and (3) Jurassic volcanic rocks of the Serra Geral Formation and aeolian sandstones of the Botucatu Formation, representing the Gondwana III Supersequence.

The volcanic-sedimentary units outcrop along the eastern margin of the basin, in contact with the rocks of Dom Feliciano Belt, while the volcanic rocks of Serra Geral Formation outcrop in the western part of the studied area (Fig. 1). The stratigraphic record documents alternating marine, glacial, and continental depositional environments. The major stratigraphic units include: the Itararé Group (Late Carboniferous-Early Permian) composed of glaciogenic diamictites, sandstones, and shales deposited during the Gondwanan glaciation (Holz et al., 2010); the Rio Bonito and Palermo formations (Permian) composed of coal-bearing strata deposited in deltaic and coastal-plain environments, marking post-glacial transgression (Holz et al., 2010); the Irati Formation (Late Permian) composed of organic-rich black shales and carbonates deposited under restricted marine conditions and considered one of the principal source rocks for hydrocarbons in the basin (Holz et al., 2010); the Santa Maria Formation (Triassic) composed of fluvial and aeolian sandstones (Scherer et al., 2023); and the Botucatu Formation (Jurassic) composed of cross-bedded aeolian sandstones indicative of an extensive desert environment prior to the onset of volcanism (Scherer at al., 2023).

One of the main geological domains in southern Brazil is the Serra Geral Formation which records the bimodal volcanism of the Paraná-Etendeka Large Igneous Province (LIP). This volcanic province covers an area of approximately 1.2×10⁶ km² in South America, extending across Brazil, Argentina, Uruguay and Paraguay, and constitutes the upper volcanic portion of the Paraná Basin (e.g. Piccirillo and Melfi, 1988). The magmatism has a bimodal character, with tholeiitic basalts and andesites as the dominant lithologies coupled with subordinate more silicic rocks (rhyodacite and rhyolite, also occurring as high-grade ignimbrites; Peate et al., 1992; Luchetti et al., 2018). The eruptive history is short, with ages constrained between 135.0 ± 0.6 Ma and 132.0 ± 0.2 Ma, with a time span of 1.6–3.0 Ma for the extrusive activity (Gomes and Vasconcelos, 2021). The Sr isotopic composition of this magmatism is extremely variable (e.g. ⁸⁷Sr $/$ ⁸⁶Sr₁₃₀ between 0.705–0.728; Peate et al., 1992), reflecting significant crustal contamination during magmas ascent and/or derivation from metasomatized mantle sources (e.g. Giovanardi et al., 2022). These volcanic rocks form a plateau in the western part of the studied area separated from other geological domains by a stiff cliff.

3 Material and Methods

3.1 Materials

A total of n= 233 samples of plant leaves were collected across the study area (Fig. 1), covering all of the geological units identified on state geological maps (downloaded from the Serviço Geológico do Brasil website https://geosgb.sgb.gov.br/geosgb/downloads_en.html, last access: 3 September 2025; geological map of Santa Catarina from Wildner et al., 2014; Rio Grande do Sul map from Wildner et al., 2008).

Sampling was designed to minimize potential contamination from agriculture, waste disposal or other anthropogenic sources. To this end, leaves were taken from trees as far as possible from anthropic activities (Spies et al., 2025). The sampled taxa include: Eucalyptus sp., Phytolacca dioica, Neltume nigra, Smilax campestris, Celtis australis, Fraxinus excelsior (Common ash) and Nerium oleander. Different plant taxa were intentionally collected at each location to account for different rooting depths at each location.

3.2 Isotopic analyses

Chromatographic separation of Sr was conducted at the MeGic lab (https://www.geochem.unimore.it/, last access: 10 October 2025) of the Department of Chemical and Geological Sciences (Unimore, Italy). Leaf samples were cleaned to remove major impurities using ultrapure MilliQ^® water and then dried for 24 h on a hotplate at 50 °C. Approximately 20 g of dried clean plant material was reduced to ash in a porcelain crucible in a muffle oven at 650 °C for 6 h. Five milligrams of ash were then dissolved in 0.5 ml of 3M HNO₃ and centrifuged prior to loading into the chromatographic column. Sr separation was performed using a standard procedure in ∼ 30 µL columns filled with the 100–150 µm bead size Eichrom Sr-spec Resin (Eichrom Technologies, LLC) as described in Argentino et al. (2021). The Sr was eluted with MilliQ^® water. Through the chromatographic separation runs, n= 5 blanks were produced to monitor possible contamination during the procedure.

Strontium isotope ratios (⁸⁷Sr $/$ ⁸⁶Sr) of 197 samples were measured with the Thermo Scientific NEPTUNE Plus multi-collector ICP-MS (MC-ICPMS) at the HUN-REN Institute for Nuclear Research (ATOMKI), Debrecen, Hungary. Isotopic ratios were corrected for instrumental mass discrimination using ⁸⁸Sr/⁸⁶Sr = 8.375209 and by applying an interference correction for ⁸⁷Rb⁺, ⁸⁴Kr⁺and ⁸⁶Kr⁺ with ⁸⁵Rb⁺ and ⁸³Kr⁺, respectively. The same method, shared between the two laboratories (Cavazzuti et al., 2021), was used to analyse n= 36 samples at the laboratory of the Centro Interdipartimentale Grandi Strumenti of Unimore using a Thermo Scientific NEPTUNE MC-ICPMS. Data reduction was performed with a modified version of the “SrDR” spreadsheet (Lugli et al., 2020). Data are available at https://doi.org/10.5281/zenodo.17988601 (Scaggion et al., 2025a).

3.3 Quality assessment of isotopic data

Poor chromatographic removal of Rb from the sample could lead to overestimation of the ⁸⁷Sr $/$ ⁸⁶Sr ratio due to the isobaric interference of ⁸⁷Rb on ⁸⁷Sr. Before isotopic analyses, the Rb and Sr concentration of the separated solutions were determined using a triple quadrupole ICPMS iCAP TQ (Thermo Fisher) at the laboratory of the Centro Interdipartimentale Grandi Strumenti of Unimore. A ¹¹⁵In solution was employed as an internal standard. Calibration curves were obtained using the multi-element standard IVICP- MS-71A (Inorganic Venture), in the range of 1–1000 ppb. Samples containing Rb higher than 0.5 ppb and with Rb/Sr ratio higher than 0.005 were reprocessed through a second chromatographic separation.

Potential contamination during the chromatographic separation was evaluated with laboratory blank samples, which were monitored during the isotopic analyses at the Centro Interdipartimentale Grandi Strumenti of Unimore. The resulting signals were indistinguishable from the analytical blank (i.e. 4 % suprapur HNO₃).

To ensure homogeneous reproducibility of the data, both laboratories followed the same analytical protocol, already used in previous publications (e.g. Cavazzuti et al., 2021). The international Sr carbonate isotope Standard Reference Material NBS987, produced by the National Institute of Standard and Technology (NIST), was used by both laboratories as an external standard for instrumental bias correction. Each session includes the analysis of three samples bracketed by the NBS-987 reference material. The instrumental bias was corrected independently for each analytical session by normalizing to the accepted NBS-987 value of 0.710248 (McArthur et al., 2001). The average ⁸⁷Sr $/$ ⁸⁶Sr value of the NBS-987 across all sessions was 0.710329 ± 0.000006 at the HUN-REN laboratory and 0.710250 ±  0.000007 at Unimore.

3.4 Dataset compilation

The dataset of the bioavailable Sr isotopic data analysed in this work and published on Zenodo (Scaggion et al., 2025a) includes: Sample ID; ⁸⁷Sr $/$ ⁸⁶Sr; 1SE; Latitude (in decimal degrees); Longitude (in decimal degrees); State; Country; Geographic Coordinate Systems; Coordinate Instrument; Coordinate Type; Bedrock Type; Bedrock Age; Geological Unit; Sample Type; Sample Part; Collected Sample Source; Preparation Step; Laboratory; Technique; Instrument; Reference Material; Normalizing Value; Reference of the Normalizing Value; Avg. ⁸⁷Sr $/$ ⁸⁶Sr Reference Material; 1SD.

3.5 Statistical analysis and geospatial modelling

We built the Sr isoscapes of South America by combining the new Sr isotope dataset of 233 plant samples with literature data. The latter includes the Isoarch South America dataset (Salesse et al., 2018, 2020), excluding rock samples, and additional published sources (see Table S1 in the Supplement) for a total of 1588 data points. For the continental-scale compilation, we included both environmental samples (plants, lichens, soils, and water) and archaeological human and faunal samples identified as local in the original publications. Non-local individuals were excluded. Metadata relating to coordinates and material categories were harmonized, and duplicates were consolidated by site. Uncertainty values not reported in the original publications or datasets were excluded. All primary datasets used in the compilation, including the open-access IsoArcH repository, Borges et al. (2021), Plomp (2021), Stantis et al. (2024) and Seferidou (2025), have been integrated with all additional sources and reported in Table S1. Overall data are from: Brass (1976), Palmer and Edmond (1989), Edmond et al. (1995, 1996), Hieronymus et al. (1993), Gaillardet et al. (1997), Négrel and Lachassagne (2000), Grove et al. (2003), Knudson et al. (2004, 2013), Brunet et al. (2005), Pasquini et al. (2005), Knudson and Price (2007), Knudson (2008), Martins et al. (2008), Andrushko et al. (2009), Bastos (2009, 2014), Fiege et al. (2009), Hermenegildo (2009), Poszwa et al. (2009), Bastos et al. (2011, 2014, 2015, 2016, 2019, 2021), Laffoon et al. (2012), Lee et al. (2013), Machado (2013), Malaspinas et al. (2014), Pouilly et al. (2014), Santos et al. (2014), Oppitz (2015), Strauss et al. (2016), Barberena et al. (2017, 2019, 2020, 2021), Chala-Aldana et al. (2018), Slovak et al. (2018), Plomp et al. (2019), Moquet et al. (2020), Serna et al. (2020), Azevedo et al. (2021), Washburn et al. (2021), Fernandez et al. (2022), Quaggio et al. (2022), Torres-Rouff et al. (2022), Avigliano et al. (2023), Seferidou et al. (2023), Silva et al. (2023), Kafino et al. (2024), Martinelli et al. (2025), Loponte et al. (2025), Scaggion et al. (2025a). The final dataset used for the model includes 1821 data from 883 sites (i.e., same coordinate samples) divided in 9 categories as follows: plant, fauna (archaeological remains, non-local excluded based on publication results), human (archaeological remains, non-local excluded), shell, snail, lichen, water and soil. The dataset is updated to September 2025. We constructed three different models, the first utilizing the whole dataset (hereafter “All” dataset), the second limited to plant + fauna + lichen + human (661 sites) and the third using only the plant + lichen samples (531 sites). After a preliminary run of the random forest model, we noticed remarkably large residuals for ⁸⁷Sr $/$ ⁸⁶Sr values higher than 0.900. We thus excluded from the models n= 5 “anomalous” values of river waters from the Parguaza Batholith (Venezuela; ⁸⁷Sr $/$ ⁸⁶Sr > 0.900) from Edmond et al. (1995). In the main text, we mainly focus on the plant + lichen isoscape, while the other models are reported in the Text S1 in the Supplement. Note that the quantile colour scale is the same for each model, but rescaled to the values of the highlighted areas.

The dataset was regressed using the random forest machine-learning model from Bataille et al. (2020) in the R free software (version 4.0.5), which reconstructs the Sr isotopic ratio variability with a resolution of 1 km². The method uses external predictors variables to generate multiple decision trees using each time a random subset of data and covariates. For the models generation we used the global variables presented in Bataille et al. (2018, 2020) and Reich et al. (2024) as external predictors (Table 1), filtering for each model the variables with high correlation ($-\mathrm{0.80}>r>\mathrm{0.80}$, Fig. 2) to avoid multicollinearity in the model and inflation of the variables importance (Strobl et al., 2008), as in Scaggion et al. (2025b). This resulted in the use of n= 14 external predictors, which were identical for the three models: r.elevation, r.bulk, r.fert, r.dust, r.pet, r.cec, r.meanage_geol, r.ssaw, r.ph, r.srsrq1, r.clay, r.age, r.volc and r.GUM.

We thus optimized the mtry parameter for the random trees construction at 5 random variables at the time. The spatial-uncertainty map was generated with the ranger package of Wright and Ziegler (2017), which calculates a quantile random forest regression, and then halving the random forest q0.84 − q0.16 difference (i.e., lower and upper limits of a ∼ 68 % interval: Funck et al., 2021, and Armaroli et al., 2024).

The models were evaluated with a 10-fold cross validation by calculating the Root Mean Squared Error (RMSE). This method requires partitioning of the data into multiple subsets (folds) which are used as unknown to train and test the model, ensuring that it is not biased by single subsets. Using this approach, the model's accuracy and precision are tested across multiple subsets of data, providing a comprehensive assessment of the model's predictive capabilities.

Table 1 List of global variables used in the model as external predictors before filtering.

Figure 2 Pearson's correlation coefficients (r) of the external predictors used in the random forest model.

4 Results

4.1 ⁸⁷Sr $/$ ⁸⁶Sr of novel plant samples

Collected samples from the studied area yielded a broad range of ⁸⁷Sr $/$ ⁸⁶Sr ratios, between 0.70521 and 0.76039 with a median of 0.71360 (Fig. 3). The dataset shows an asymmetric distribution (kurtosis = 13.03, skewness = 2.77) with a long tail towards more radiogenic values (Fig.  3). A Shapiro-Wilk test indicates significant deviation from normality (W=0.733, p<0.01). To facilitate interpretation, results of novel plant samples are presented here according to the underlying geological units or rock type as in Fig. 1 and Table 2. We stress however that plants were collected from soils developed on these units.

Figure 3 Histogram of the whole dataset utilized to construct the isoscape model and of the new plant data. Boxplots of the plant samples from SC and RGS clustered by geological units as in Fig. 1: Q4 – Holocene sediments; Q1 – Pleistocene sediments; EN – Paleogene sedimentary; K1α – Serra Geral Formation, felsic; K1β – Serra Geral Formation, basic; K1 – Cretaceous sedimentary; J – Jurassic sedimentary; T – Triassic sedimentary; P3 – Late Permian sedimentary; P1 – Early Permian sedimentary; Ppd – Permian andesite; C2 – Carboniferous sedimentary; Np3 – Ediacaran metasedimentary; Np3γ – Ediacaran metaigneous; Np2γ – Cryogenian metaigneous; Np1 – Tonian metamorphic; PP – Paleoproterozoic basement.

Table 2 Bioavailable Sr isotopic average composition of this work plant samples divided for the different geological super-units of the region as reported in Fig. 1. The order follows the stratigraphic order.

Plants from Holocene sediments (Q4) show values between 0.70843 and 0.74054, with an average of 0.71569 ± 0.00870 (n= 21; error calculated as standard deviation; Fig. 3; Table 2) while those from Pleistocene sediments (Q1, Q2, and Q3) have a narrower range from 0.70940 to 0.72650, with an average of 0.71306 ± 0.00481 (n=11), similar to only one Paleogene sample (EN) with a value of 0.71273 ± 0.00002 (Table 2). Plants from both these sedimentary groups include anomalous values, as outliers of the distribution (Fig. 3).

Samples from the Cretaceous felsic volcanics units of the Paraná LIP (K1α) show higher values (between 0.70766 and 0.72147; average of 0.71421 ± 0.00412, n=15) than those from the basic volcanics units (K1β; between 0.70666 and 0.71741; average of 0.71156 ± 0.00226, n=46; Fig. 3, Table 2). Plants from Cretaceous sediments of the Paraná Basin (K1) have similar K1α values, (0.71191 and 0.71578; average of 0.71351 ± 0.00202; n=3; Table 2), almost indistinguishable from those from Jurassic sediments (J; between 0.71045 and 0.71515; average of 0.71335 ± 0.00151; n= 8; Fig. 3, Table 2). Plants from other sedimentary units of the Paraná Basin (i.e., T, P1, and P3 units) show progressive radiogenic enrichment over time (averages at 0.71432 ±  0.00203; 0.71429 ± 0.00292; 0.71516 ± 0.00204, respectively). Samples from Permian volcanics (Ppd) have comparable values, between 0.71118 and 0.71554 with an average of 0.71336 ± 0.00309 (n=2; Table 2). In contrast, plants from Carboniferous sedimentary rocks (C2), have a distinct and more radiogenic signature between 0.71796 and 0.74832, average at 0.72604 ± 0.01263 (n=  5; Fig. 3, Table 2).

Samples collected from units of the Mantiqueira Province show broader ranges and generally higher values than those from Paraná Basin samples (except for C2 unit; Fig. 3). Plants from Ediacaran sediments of the Dom Feliciano Belt (Np3) yielded ⁸⁷Sr $/$ ⁸⁶Sr values between 0.70733 and 0.71871, with an average of 0.71364 ± 0.00428 (n= 8), while plants from metaigneous units (Np3γ) show higher and more variable values (0.7064 to 0.74790, with an average of 0.71892 ± 0.00830; n= 51; Fig. 3, Table 2). Samples from Cryogenian metaigneous units from the São Gabriel Terrane (Np2γ) show the widest range in the dataset (from 0.70987 to 0.76039; average of 0.72522 ± 0.01727; n= 10; Fig. 3, Table 2). Plants from the Tonian metasedimentary units (Np1) show values between 0.70769 and 0.72978 with an average of 0.71601 ± 0.00732 (n= 9; Fig. 3, Table 2). Finally, the sample from the Paleoproterozoic basement (PP) yielded the lowest value of the entire dataset (0.70521 ± 0.00002; Table 2).

4.2 Random forest models

The presence of isotopic outliers in our plant dataset, with ratios higher than those expected for unconsolidated coastal deposits, can be explained by the heterogeneity of the Holocene and Pleistocene deposits (fluvial, marine, aeolian) on which the plants grew, and the characteristics of the parental lithologies. Since these values reflect geological variability, they were retained in the model. In contrast, the extremely radiogenic values for river waters (> 0.900) were excluded.

The three constructed isoscape models clearly distinguish the most radiogenic areas of the SC and RGS, which often correspond to the oldest geological units in the region (Figs. 4 and 5). In all models, the associated errors are higher in these radiogenic areas (Figs. 4 and 5), indicating lower predictive power of the models in these regions.

Figure 4 Spatial distribution of strontium isotopes in South America based on plant + lichen samples (n= 531 sites). Isoscape random forest model (a), associated spatial uncertainty map (b) and sample distribution (c).

Figure 5 Spatial distribution of strontium isotopes in Santa Catarina and Rio Grande do Sul states based on plant + lichen samples (n= 531 sites). Isoscape random forest model (a) and associated spatial uncertainty map (b).

Evaluation metrics reported in Table 3 show significant differences between the datasets: the model based on the “All” dataset exhibits the lowest predictive power, with a RMSE of 0.0086 and an R² of 0.639, while the models built on the plant + fauna + lichen + human and the plant + lichen datasets show similar RMSE values of 0.0049 and 0.0054 respectively, and an R² of ca. 0.76 each. The lower variance observed in these latter models likely reflects the smaller number and lower variability of samples compared to the “All” dataset, which includes outliers with ⁸⁷Sr $/$ ⁸⁶Sr higher than 0.818 in water samples (Fig. 3).

Table 3 Evaluation metrics of the models.

The evaluation of the predictive power of the independent variables using the Increase in Node Purity (IncNodePurity; Fig. 6; Text S1 in the Supplement), shows that in the plant + fauna + lichen + human and plant + lichen models, the most relevant variables are r.srsrq1 (predicted first quartile of the global ⁸⁷Sr $/$ ⁸⁶Sr bedrock model; Bataille et al., 2018) and r.ph (soil pH in H₂O solution × 10; Bataille et al., 2020), both with IncNodePurity values just below 0.012 (Fig. 6), followed by r.meanage_geol (GLiM age attribute in Myrs; Bataille et al., 2020) and r.cec (cation exchange capacity; Bataille et al., 2020) with values of around 0.009, while all other variables have a contribution lower than 0.004 (Fig. 6).

Figure 6 Predictive variables and model accuracy for random forest model based on plant + lichen samples. Increase in Node Purity (IncNodePurity) of the n= 14 external predictors (a); cross-validation scatterplot of observed vs. predicted ⁸⁷Sr $/$ ⁸⁶Sr values for the test-validation sample splits (b); model residuals vs. predicted ⁸⁷Sr $/$ ⁸⁶Sr ratios (c).

Notably, a 2D-Partial Dependence Plot (PDP, Fig. 7) shows a nonlinear interaction between r.ph and r.srsrq1, and indicates that the RF model predicts higher ⁸⁷Sr $/$ ⁸⁶Sr ratios when pH is low and – as expected – when the r.srsrq1 bedrock model value is high (Fig. 7).

Figure 7 2D Partial Dependence Plot for predictive variables r.srsrq1 and r.ph for random forest model based on plant + lichen samples.

In contrast, within the “All” model the most predictive variable is r.pet (global potential evo-transpiration; Bataille et al., 2018) with an IncNodePurity value equal to 0.038, followed by r.ph at 0.020, r.age (Terrane age attribute in Myrs; Bataille et al., 2020) at 0.018, and r.cec and r.srsrq1 at 0.015 (Text S1 in the Supplement).

5 Discussion

5.1 Validation of the novel dataset against bedrock geology

Several inputs and processes govern the bioavailable Sr isotopic signature of local areas. However, considering a regional dataset, the bedrock geology is one of the main factors which govern the isotopic distribution (Vitousek et al., 1999; Voerkelius et al., 2010), and thus its relationship with the dataset could be a useful proxy to validate it.

The asymmetric statistical distribution of ⁸⁷Sr $/$ ⁸⁶Sr values observed in the novel dataset (Fig. 3), characterized by a long tail toward more radiogenic ratios, reflects the geological complexity of the area and the variety of contributing lithologies (Cryogenian igneous rocks, Ediacaran sedimentary sequences, Cretaceous volcanic flows, etc.). The presence of anomalous values, i.e., ratios higher than those expected for unconsolidated coastal deposits and identified as outliers of the distribution, among plants growing on Holocene and Pleistocene sediments can be explained by the heterogeneity of the deposits (fluvial, marine, aeolian) and the characteristics of the parental lithologies.

The more radiogenic values of plants samples growth on soils from felsic igneous and metamorphic units compared to igneous basic and sedimentary units (e.g. units K1α, magmatic felsic, and K1β, magmatic basic with the same age), are consistent with well-known geological differentiation processes, as Rb commonly behaves as a more incompatible element in magmatic and metamorphic processes than Sr, thus becoming enriched during fractionation, forming rocks with higher or lower Rb content, which over time produce higher or lower ⁸⁷Sr $/$ ⁸⁶Sr ratios, as in K1α and K1β, respectively.

Another aspect related to geology, is the increase of the ⁸⁷Sr $/$ ⁸⁶Sr ratio in rocks with time, which we expect to be reflected in the derived soil and then in the local vegetation. The plant dataset presented here follows this rule, with a general increase of the ⁸⁷Sr $/$ ⁸⁶Sr ratio as the bedrock is older in age (Fig. 3; Table 2) as observed for the Paranà Basin samples vs. Mantiqueira Province ones.

All this evidence suggests that the dataset is in overall agreement with bedrock geology which is the primary source of bioavailable Sr in nature.

5.2 Limits, advantages and predicting variables of the random forest models

Across the developed models, higher spatial uncertainties are associated with the most radiogenic areas of SC and RGS, generally corresponding to the oldest geological units in the area (Figs. 4, 5). This pattern, also reported in random forest models from the literature (e.g., Bataille et al., 2018; Gigante et al., 2023; Scaggion et al., 2025b), reflects the pronounced contrast in isotope composition between radiogenic bedrock and less radiogenic end-members such as aerosols, seawater, and precipitation. Such contrast complicates the mixing relationships processes that govern the bioavailable Sr isotope signature. Another parameter that influences model accuracy and associated errors is the uneven distribution of samples on a continental scale. Nonuniform sampling tends to introduce bias towards specific geological contexts, thereby limiting the ability to generalize predictions in unrepresented areas.

The improved performance observed in the plant + fauna + lichen + human and plant + lichen models, compared to the “All” dataset, can be attributed to the reduced variability of the data, which limits the influence of outliers and improves the overall model stability. We reiterate that these differences in RMSE and R² more accurately reflect the isotopic range and internal heterogeneity of the training datasets, rather than intrinsic differences in model quality. The absence of significant differences in the predictiveness between these two best performing models may reflect (i) nonuniform and incomplete geographic coverage associated with water and soil sampling strategies, and (ii) the limited representativeness of inorganic samples compared to biological and bioderived materials in representing the local terrestrial bioavailable Sr signature, as recently reported by Spies et al. (2025). In this regard, for random forest models, a smaller but more geologically representative dataset of biological samples is preferable over a larger but uneven dataset. Indeed, the latter reduces the representativeness of the isotopic variability and, consequently, compromises the accuracy and reliability of the models. Conversely, a smaller dataset anchored to the local geological context is more stable, robust, and useful for modelling (Marsh et al., 2025).

The predicted distribution map of the bioavailable ⁸⁷Sr $/$ ⁸⁶Sr isotope ratio shows a consistent correlation with the data obtained as a function of the variables r.srsrq1 and r.ph. The combined marginal effect (Fig. 7) between these two independent variables reveals that isotopic predictions are influenced by non-linear interactions between strontium isotope ratios (r.srsrq1) and soil pH (r.ph). Higher ⁸⁷Sr $/$ ⁸⁶Sr ratios correlate with more acidic soils, a pattern consistent with pedogenesis on ancient silicate igneous-metamorphic rocks. In the Amazon region, low pH values result from intense decomposition of organic matter, which leads to strong leaching of the bedrock and, consequently, increases Sr availability in the soil. The PDP (Fig. 7) explains how the geological context, pedogenic processes, and environmental conditions shape the spatio-temporal distribution of Sr isotope values.

Compared to existing Sr isoscapes, this study substantially improves both the spatial coverage and resolution of bioavailable predictions for South America, particularly southern Brazil. Although global or continental-scale isoscapes are essential for broader applications, developing regional models based on dense, geologically informed sampling offer much higher resolution and are essential for fine-scale provenance studies, such as intra-site archaeological mobility or ecological habitat use. Despite being focused on SC and the RGS regions, our dataset allowed us to obtain a reasonably robust continental-scale model with an R² of 0.764 (Fig. 4). Our model compares positively with other continental-scale isoscapes (e.g. the bioavailable Sr isoscape of Australia, with an R² of 0.590; Dosseto et al., 2025). With this comparison, we do not intend to downplay the importance of large-scale modeling approaches, which are essential for global applications, but rather to highlight how carefully planned sampling strategies at the local level can provide robust, useful, and informative models at larger scales.

This highlights how carefully designed local sampling strategies not only strengthen regional provenance studies but also provide a reliable base for modeling at broader geographical scales.

6 Data availability

The bioavailable Sr plant isotope dataset of Santa Catarina and Rio Grande do Sul states analysed in this work, the R script and the raster files of the plant + lichen model are available at https://doi.org/10.5281/zenodo.17988601 (Scaggion et al., 2025a). Check the latest version of the dataset. The literature dataset is reported in Table S1 in the Supplement.

7 Conclusions

This work presents a comprehensive dataset of 233 bioavailable ⁸⁷Sr $/$ ⁸⁶Sr values from plant samples collected across the main geological units of Santa Catarina and Rio Grande do Sul states (southern Brazil). The measured ⁸⁷Sr $/$ ⁸⁶Sr ratios in the terrestrial biosphere are consistent with lithological patterns, reflecting the pronounced geological heterogeneity of the region.

By integrating these new data with a multi-material dataset of literature data, we developed three isoscapes of bioavailable Sr distribution for South America and the studied region (ca. 370 000 km²) using advanced random forest modelling. Among the tested approaches, the plant + fauna + lichen + human and plant + lichen models achieved the best predictive performance, while the addition of water and soil reduced model performance due to their lower representativeness of locally bioavailable Sr.

Despite the limited dataset, the resulting models show robust predictive performance even at the continental scale. This result demonstrates how well-designed local sampling strategies can ensure reliable regional models that maintain their relevance and applicability at larger spatial scales. Beyond providing a new baseline for environmental, archaeological, and palaeoecological research in southern Brazil, our study highlights an important methodological insight: strategically constructed regional datasets can play a key role in refining continental-scale isoscapes. Future efforts should focus on expanding high resolution datasets across underrepresented geological domains.