Lippia alba is used for its sedative, antidepressant, analgesic, antiviral, antimicrobial, anti-inflammatory, anthelmintic, antioxidant, antimalarial, and cytostatic effects (Carvalho et al. 2003; Aguiar et al. 2008; Hennebelle et al. 2008; Gomes et al. 2019). Adapted to a wide climatic range—from tropical to temperate zones—it thrives in diverse environments such as forests, swamps, and open fields, often growing along the margins of rivers, ponds, and lakes (Moldenke 1965; Table 1). Lippia turbinata is used as a diuretic, emmenagogue, tonic, sedative, and for stomach infusions, with antimicrobial activity reported (Hernández et al. 2000; Coll Aráoz and Ponessa 2007; Pérez-Zamora et al. 2016). It occurs in the Monte, Espinal, and Chaco biogeographic provinces, with smaller occurrences in the Yungas and Pampean provinces (Múlgura et al. 2012; Table 1). Lippia turbinata inhabits alluvial terraces of mountain streams in certain areas (Andersen et al. 2006) and performs well in shallow soils (Table 1). Infusions of Salimenaea integrifolia leaves and flowers are used for diuretic, emmenagogue, antibiotic, and sedative effects, and to treat bronchopulmonary conditions (Catalán et al. 2021); effects against gastritis inflammation were reported (Marcial et al. 2014). This plant is used in food preparation, including appetizers, non-alcoholic beverages, and teas (Código Alimentario Argentino 1969). This species inhabits stony hills or sites with gentle to very steep slopes, arid fields, or xeric shrublands, thriving in incipient soils with poorly developed profiles and low levels of organic matter and nutrients (Table 1).

The area of this study (4–32°S, 38–70°W) covers approximately 8,788,573 km², from eastern Bolivia and Brazil, Paraguay to central and northern Argentina (Fig. 1; Table 1). This spans the complete distribution of Lippia turbinata and Salimenaea integrifolia; for L. alba it covers the distribution in southern South America. Following Cabrera and Willink (1980), it comprises different biogeographic provinces from two domains.

The Chacoan Domain, encompassing the Caatinga, Chaco, Espinal, Monte, Pampean, and Prepuna provinces, is defined by a predominantly continental climate with moderate to scarce rainfall, mild winters, and warm summers. Vegetation varies from deciduous xerophilous forests and shrublands to grass steppes and xerophytic flora. The Caatinga in north-eastern Brazil features clear forests and open shrublands, with 400–750 mm of annual rainfall and temperatures of 26–27°C. The Chaco, spanning northern Argentina, central Paraguay, south-eastern Bolivia, and parts of Brazil, receives 500–1200 mm of rainfall annually, with temperatures of 20–23°C, and is dominated by deciduous xerophilous vegetation. The Espinal, encircling the Pampean province in central-eastern Argentina, has thorny and microphyllous vegetation, 340–1170 mm of rainfall, and temperatures of 15–20°C. The Monte, a dry steppe in western Argentina, consists of xerophytic shrubs, receiving 80–350 mm of rainfall with temperatures of 13–15.5°C. The Prepuna, in north-western Argentina’s Andean foothills, lies at 1000–3400 m a.s.l., with summer rains and shrublands mixed with tree-like cacti. Finally, the Pampas, covering eastern Argentina, Uruguay, and southern Brazil, features a temperate climate, 600–1200 mm of rainfall (decreasing southward), and temperatures of 13–17°C.

The Amazonian Domain, comprising the Amazonian, Atlantic, Cerrado, Paranaense, and Yungas provinces (Cabrera and Willink 1980), spans much of South America and is characterized by a warm, humid climate with dense vegetation rich in biodiversity. The Amazonian province in northern Brazil and nearby South American countries has a warm, humid climate with steady 26°C temperatures and 2000–2600 mm of annual rainfall, supporting lush tropical rainforests. The Yungas province spans the eastern Andes from Venezuela to northern Argentina (500–3000 m a.s.l.), featuring cloud forests, montane forests, and grasslands in a cool, humid climate with abundant rainfall and fog that decrease with altitude. The Cerrado, dominating central Brazil and extending into northern Paraguay, is an open woodland region at elevations of 500–1000 m, with annual rainfall of 1200–2000 mm, average temperatures of 21–25°C, and low forests. The Paranaense province spans southern Brazil, eastern Paraguay, and north-eastern Argentina, with subtropical rainforest, annual rainfall of 1500–2000 mm, temperatures of 16–22°C, and altitudes from 1000 to 3000 m a.s.l. The Atlantic province forms a narrow strip along Brazil’s eastern slopes, characterized by tropical rainforests, a hot and very humid climate, and annual rainfall exceeding 2000 mm, reaching up to 4000 mm in some areas, with temperatures of 19–25°C.

We performed studies of three MAPs (Fig. 1; Table 1). The georeferenced points were compiled from herbarium records BAB, BAF, CORD, CTES, SI, LPB (abbreviations follow Index Herbariorum: Thiers 2025) and from the Documenta Florae Australis base (www.darwin.edu.ar/iris), a database of Argentinean vascular flora. The accurate identification of each specimen was confirmed, and any specimens with unverified identities or missing collection site information were omitted. To ensure the accuracy of the geographic coordinates, the entries were reviewed, resulting in the removal of duplicate points. To avoid the spatial autocorrelation effect in the SDM, records less than 1 km apart were removed using the R package Wallace v.1.0.6.1 (Kass et al. 2018) implemented in R v.3.6.1. After filtering 693 occurrences, a total of 376 points were obtained: 127 for L. alba, 178 for L. turbinata, and 71 for S. integrifolia (Fig. 1; Suppl. material 1). The resulting occurrence points are adequate for conducting maximum entropy (MaxEnt) analyses (Hernandez et al. 2006; Phillips et al. 2006).

For the development of current SDMs, environmental layers including 19 bioclimatic variables were sourced from the WorldClim2 database (http://www.worldclim.org) at a spatial resolution of 30 arc-seconds per pixel (approximately 1 km²) (Fick and Hijmans 2017). To avoid overestimation of climatic influences due to multicollinearity, highly absolute correlated variables (r = 0.8) were excluded via Pearson’s correlation analysis (Suppl. material 2) (Herrando-Moraira et al. 2020; Scrivanti and Anton 2020; Zhang et al. 2023). Correlation tests were performed in Infostat v.2020 (Di Rienzo et al. 2020). In addition, we used edaphic variables extracted from the SoilGrids database (https://www.isric.org). Out of seven possible edaphic variables, two were selected: SND (sand content) and CLYPPT (clay content). This choice was made because sandy and clayey soils are the most suitable for the development of these plant species (Múlgura et al. 2012). After a secondary analysis using Maxent, only the SND variable was retained for the study, as CLYPPT was excluded due to its low permutation importance. For all species we selected the following 5 to 6 variables: Bio1 (annual mean temperature), Bio3 (isothermality), Bio7 (temperature annual range), Bio15 (precipitation seasonality), and SND (sand content); for L. turbinata and S. integrifolia, we added Bio12 (annual precipitation). Therefore, climatic and edaphic variables that were considered biologically important and directly relevant to these species were selected (Ali et al. 2020; Alem et al. 2022; Hu et al. 2024). All selected layers were clipped according to the calibration area, which was defined based on the Biotic-Abiotic-Mobility (BAM) theoretical approach (Soberón and Peterson 2005; Peterson et al. 2008). For this study, the calibration area was defined as a polygon with a buffer zone of 200 km around the occurrence points of each species.

For future projections, three Atmospheric-Ocean Global Circulation models (AOGCMs) were used: CCSM4 (Community Climate System Model, version 4), CCM3 (Community Climate Model, version 3), and HadGEM3 (Hadley Center Global Environmental Model, version 2). These models have been used to evaluate the impact of climatic change on the medicinal plant Valeriana carnosa Sm. in southern South America (Nagahama and Bonino 2020). For each model, three climate change scenarios (Representative Concentration Pathways Scenario, RCP) 2.6, 4.5, and 8.5 W.m^-2 were evaluated for the period 2070 (2061–2080). These scenarios are distinguished by their radiative forcing projections for the year 2100, with values of 2.6 and 8.5 representing the most optimistic and most pessimistic emission pathways, respectively (Meinshausen et al. 2011; Van Vuuren et al. 2011). Climate variables for the future were obtained from the CHELSA database (Climatologies at High Resolution for the Earth’s Land Surface Areas) at a spatial resolution of 30 arc-seconds. These three different AOGCMs and RCPs were chosen to account for the uncertainty they introduce when modelling future conditions (Diniz-Filho et al. 2009; Nori et al. 2011). Since climate change generally has little impact on soil variables, we used SND as stabilizing factors in the development of future model projections (Gilani et al. 2020; Liu et al. 2022).

Finally, current and future altitude predictions were calculated.

To construct the current and future projections, the maximum entropy algorithm was implemented using MaxEnt v.3.3.3k (Phillips et al. 2017), following the recommendations of Phillips et al. (2006). This method is widely used and highly regarded for its robustness, as well as its broad application in modern scientific research, particularly in assessing the potential and future distribution of plant species (Ngarega et al. 2024; Xia et al. 2023; Recopuerto-Medina et al. 2024). The data were analysed using the following parameters: a maximum of 500 iterations, a maximum of 10,000 background points, a convergence threshold of 10^-5, and a regularization multiplier of 1. For model training, 75% of the random localities were used, while 25% were allocated for testing the model through bootstrap with 10 replicates. To evaluate the models, the AUC test (area under the receiving operating characteristic curve) was used according to Peterson et al. (2007) and Lobo et al. (2008). An AUC of 0.5 indicates that the model predictions are no better than random predictions, below 0.5 worse than random, 0.5–0.7 low performance, 0.7–0.9 moderate performance, and above 0.9 high performance (Peterson et al. 2011). Variability across AOGCMs and RCP scenarios was assessed through the standard deviation of suitable distribution areas.

The variables’ contribution to the SDM was assessed using percent contribution, permutation importance, the jackknife test, and the response curves generated by MaxEnt (Phillips et al. 2017). The 10 models of each taxon were stored in ASCII raster format and imported into QGIS v.3.24.1 (QGIS Development Team 2022) to produce a strict consensus map, which yielded four levels of habitat suitability: excellent (0.6–1), good (0.4–0.6), fair (0.2–0.4), and poor (< 0.2) (Abolmaali et al. 2018; Xie et al. 2021). Finally, the MaxEnt projections were reclassified to convert the continuous output into a binary presence/absence (0–1) map. To obtain the total potential area in km² for present and future predictions, we employed the 10^th percentile training presence threshold. This means that MaxEnt ranked the habitat suitability values of the grid cells species with records that were used to train the model, and the threshold was set at the 90^th percentile of these values. This approach assumes a 10% omission error in the presence records and is commonly used in conservation studies (Abba et al. 2012; Scrivanti and Anton 2021). Subsequently, to estimate the expansion, retraction, and stable areas in the most suitable regions, we applied a threshold of 0.6 to construct binary maps, averaging the results from the three models for each scenario (Peralta et al. 2024).

We obtained current and future potential distribution for Lippia alba, L. turbinata, and S. integrifolia based on climatic and soil factors. The SDM results indicated that the models performed satisfactorily, with AUC training and testing values ranging from 0.802 to 0.894 of current, RCP 2.6, 4.5, and 8.5 scenarios (Suppl. material 3); these values were higher than the random prediction AUC value to assess the robustness of projected distribution changes, we evaluated variability across AOGCMs and RCPs by calculating the standard deviation of suitable areas, which ranged from 0.016 to 0.044 (Suppl. material 3). Mean values were 0.0289 for Lippia alba, 0.0305 for L. turbinata, and 0.0335 for S. integrifolia, indicating limited dispersion around mean projections and supporting model reliability (Suppl. material 3).

The projected potential current distribution of L. alba covers an area of approximately 2,804,320 km² (Fig. 2; Suppl. material 4). The most suitable areas (habitat suitability value > 0.6) cover 2,550,640 km² (Table 2) and are found in central-eastern Argentina, western Uruguay, Paraguay, central Bolivia, and south-eastern and north-eastern Brazil predominantly in Chacoan, Espinal, Pampean, Paranaense, Caatinga, Atlantic, and Amazonian biogeographic provinces. The distribution in Cerrado and Yungas provinces had lower probability.

For L. turbinata, the potential current area is approximately 768,314 km² (Fig. 2; Suppl. material 4) and the most suitable areas (habitat suitability value > 0.6) cover 671,851 km² (Table 2) and are predominantly distributed in central and western Argentina in Chacoan, Espinal, Monte, Pampean, and Yungas provinces.

For S. integrifolia, the area is about 240,187 km² (Fig. 2; Suppl. material 4); the most suitable areas (habitat suitability value > 0.6) cover approximately 196,412 km² (Table 2) and are primarily located in the mountainous regions of central and north-western Argentina, in Sierras Pampeanas and Sierras Subandinas, covering the Monte, Chacoan, and Prepuna/Puna regions, with lower probabilities in the Espinal.

Predictor variable importance was assessed using percent contribution, permutation importance, and jackknife tests (Table 3; Suppl. material 5). For L. alba, Bio15 (precipitation seasonality) and Bio3 (isothermality) were the variables that contributed the most to the models, with percentage ranges between 17–53% and 35–57%, respectively (Suppl. material 6A). Concerning L. turbinata, the variables Bio1 (mean annual temperature) and Bio12 (annual precipitation) were the most important for the models. Bio1 ranges from 13–22°C and Bio12 covers a range of 180–1200 mm (Suppl. material 6B). For S. integrifolia, Bio1 and Bio12 contributed the most to the models. Bio1 has a temperature range of 15–19°C and Bio12 covers a range of 97–630 mm (Suppl. material 6C). For the three species, Bio15 and SND have a significant influence on the modelling outcome.

We performed future projections of species distribution in the 2070s (2061–2080) under the RCPs scenarios 2.6, 4.5, and 8.5, focusing on areas with the highest suitability (habitat suitability value > 0.6) (Fig. 3; Table 2); for modelling with a 10^th percentile training presence logistic threshold, see Suppl. materials 4 and 7.

For L. alba, a decrease in suitable distribution areas is predicted under all three RCP scenarios, with reductions ranging from 4.60 to 23.23%. These reductions are expected primarily along the northern edge of its distribution, as well as in the eastern and southern regions of Brazil, particularly within the Atlantic and Caatinga provinces. Both the most optimistic low-emission scenario (RCP 2.6) and the extreme high-emission scenario (RCP 8.5) showed the greatest decrease in area. However, the species could expand into the southern Paranaense province and western Argentina, extending into the Yungas.

The potential distribution of L. turbinata under the three RCP scenarios predicts a net increase in environmental suitable areas (15.12 to 19.86%), with expansion toward all boundaries of its distribution. An exception is the loss of suitable areas in central-western Argentina, including the Precordillera and some Andean mountains in Mendoza and San Juan, as well as around the eastern Bermejo River and the confluence of the Paraná and Uruguay rivers.

For S. integrifolia, a slight net increase in suitable areas is predicted under the three RCP climate scenarios (1.48 to 3.46%), primarily at the eastern edge of its distribution, extending into the valleys and lowlands of the Sierras Pampeanas and Subandinas mountain ranges. On the other hand, there is a retraction of suitable areas along the western edge of its distribution.

In the projections, the current maximum elevation is overestimated relative to the actual altitudinal range of each of the three species (Table 1; Suppl. material 8). Looking ahead, L. turbinata tends to increase its maximum elevation under all future scenarios, whereas L. alba and S. integrifolia show a tendency to decrease their upper elevation limits as climatic conditions worsen across scenarios (Suppl. material 8).

SDM predictions under current climatic conditions for all species align with the observed and published plant distributions (Mulgura et al. 2012; O’Leary et al. 2023) providing additional support for the accuracy of the models.

Potential distribution of L. alba partially matches the Seasonally Dry Tropical Forests from the Brazilian Caatinga to western Bolivia (Prado 2000), but also encompasses other provinces, such as Chaco and Espinal, all of which share precipitation seasonality, characterized by well-defined dry and wet seasons (Cabrera and Willink 1980; Prado 2000). In this sense, precipitation seasonality (Bio15) together with isothermality (Bio3) are critical factors for its establishment and survival; blooming and fruiting are associated with warm and humid periods, while seed germination requires temperatures above 25°C (Múlgura et al. 2012; Bonilla et al. 2013).

The spatial distribution of L. turbinata is primarily shaped by mean annual temperature (Bio1) and annual precipitation (Bio12). This species thrives in mesophytic and xerophytic environments with precipitation ranging from 180 to 1200 mm and an optimal temperature range of 13–22°C. Its predicted range aligns with expected distribution in central and northern Argentina, from the northern Monte to Yungas (Múlgura et al. 2012). It enters the pre-flowering stage from mid to late spring and blooms in early summer (Múlgura et al. 2012) and germination requires temperatures above 5°C and high light intensity (Galíndez et al. 2017).

The current distribution of S. integrifolia is also influenced by mean annual temperature and annual precipitation, with ranges lower than L. turbinata, that allow it to grow in mesophytic and xerophytic environments in mountain systems from central Argentina and Bolivia up to 3000 m a.s.l. Since it requires temperatures above 5°C to germinate (Galíndez et al. 2017), its intolerance to extremely low temperatures would prevent it from expanding to higher latitudes. It has a prolonged flowering period, entering the pre-flowering stage from mid to late spring and blooming from late spring to early summer (Múlgura et al. 2012; Brunetti 2017; Leiva and Brunetti 2022). Since L. turbinata and S. integrifolia share several ecological requirements, interspecific competition should be considered when interpreting the modelling results (Hu et al. 2022). Climatic variables such as isothermality and precipitation seasonality delimit suitable areas for different MAPs that are threatened or endangered worldwide (Karami et al. 2022; Zou et al. 2023; Hosseini et al. 2024; Sarikaya et al. 2024).

Climate change simulations in the study area indicate a general warming trend in the coming decades (Guerrero and Agnolin 2016; Zeballos et al. 2020; Salariato et al. 2022). Both RCP 4.5 and 8.5 scenarios predict a temperature increase of 0.5–1°C by mid-century and 1–4°C by the end of the century, with more pronounced warming from the northern and western regions of Argentina to Brazil (Barros et al. 2015; Brêda et al. 2020; IPCC 2023). Additionally, an intensification of summer temperatures and extreme temperature events is expected (Castellanos et al. 2022). In terms of precipitation, the northern and central-western regions of Argentina are projected to experience an increase of approximately 100 mm per year (Barros et al. 2015; Brêda et al. 2020; IPCC 2023). In contrast, central and north-eastern Brazil are expected to face increased drought, dryness, and aridity, with a reduction in mean precipitation towards the eastern regions (Castellanos et al. 2022; IPCC 2023).

Despite we expected a contraction of their distribution area, the species studied presented both gains and losses of territory according to different climate change scenarios. Lippia alba is the species predicted to contract its range and shift to lower elevations. The predicted reduction of L. alba at the northern edge of its distribution, particularly within the Caatinga and Atlantic provinces, aligns with regions projected to undergo substantial temperature rises and heightened aridity in the coming decades (Brasil 2016; Castellanos et al. 2022). The connection zone between the Caatinga and other regions along Brazil’s east coast was not included in our models due to the lack of reliable georeferenced points for calibrating the buffer zone; however, a loss of suitable habitats is anticipated. The sensitivity of L. alba to isothermality and precipitation seasonality suggests it would be vulnerable to extreme temperature fluctuations, as well as changes in drought and precipitation regimes. Various MAPs would be exposed to local or regional extinctions linked to climate change (Rodríguez-Cravero et al. 2017; Karami et al. 2022; Peralta et al. 2024). The predicted loss of suitable areas in the Caatinga has also been observed for other species (Cavalcante et al. 2020; Simões et al. 2020). Additionally, ecosystems such as the Caatinga, Chaco, and Espinal face severe environmental pressures, including deforestation and habitat loss driven by human activities (Beuchle et al. 2015; Rubio et al. 2022). Therefore, conservation strategies should prioritize the preservation of key areas and the promotion of biological corridors linking the Caatinga to the west of Bolivia. In the Caatinga, protected areas like Catimbau and Chapada Diamantina National Parks would provide some level of safety for the species. Furthermore, we recommend ex situ conservation of germplasm in vulnerable areas, along with in situ cultivation in favourable regions such as the Chaco and Espinal.

For both L. turbinata and S. integrifolia, an expansion into novel geographic areas is predicted, highlighting their adaptability to changing environmental conditions. However, both species are expected to contract in the western mountainous regions, where higher temperatures are anticipated compared to the central areas (Barros et al. 2015). Geographical features, such as bodies of water or mountain chains, may limit or even prevent their expansion (Croteau 2010; Nakazawa 2013). Specifically, L. turbinata is predicted to retract in areas near the Bermejo River and the confluence of the Paraná and Uruguay rivers; the latter are expected to experience an increase in the frequency and duration of fluvial floods (Barros et al. 2015). Regarding altitude, projections differ: despite the overestimation of net values, L. turbinata is expected to expand its altitudinal range, whereas S. integrifolia is projected to shift to lower elevations, coinciding with its movement towards valleys and lowlands. Although it is generally expected that mountain species migrate to higher elevations (Elsen and Tingley 2015), various shifting patterns have been predicted in montane plants from Africa, Asia, and America (DeChaine et al. 2013; Forester et al. 2013; Rodríguez-Cravero et al. 2017; You et al. 2018; Asase and Peterson 2019). Finally, these species could potentially compete, as their distributions overlap partially and they share some similar ecological requirements.

The apparent absence of specialized dispersal structures in the three species analysed is a key factor when interpreting climate change projections. The fruit, enclosed within a persistent calyx and dehiscing into two readily separable mericarps (Múlgura 2003), suggests a schizocarpic dispersal unit with limited adaptation for anemochory or zoochory. This morphology is typically associated with passive, unspecialized secondary dispersal, implying short dispersal distances. Such a trait could significantly limit these species’ capacity for rapid range shifts in response to climatic changes.

Furthermore, several reproductive traits—such as pollination dependence, flowering variability, and seed dormancy—may also constrain their distribution and are not accounted for in climate-only models. For example, L. alba and possibly S. integrifolia exhibit obligate outcrossing systems with self-incompatibility, relying mainly on insect pollinators such as Hymenoptera, followed by Lepidoptera, Diptera, and others. The composition of floral visitors is influenced by key climatic variables, including temperature, humidity, wind speed, and light availability (Venâncio et al. 2016). In addition to its pollination constraints, L. alba shows low seed germination rates under experimental conditions, likely due to dormancy or reduced viability (Pimenta et al. 2007; Brunetti 2017). Its germination is light-dependent (Galíndez et al. 2017), a likely adaptation that ensures small seeds germinate only when near the soil surface (Bond et al. 1999). For L. turbinata, Morales and Galetto (2003) reported a low rate of spontaneous self-pollination and a high proportion of fruits resulting from natural pollination, although the exact nature of its compatibility system remains unresolved.

Research on SDMs for MAPs in South America is still limited. Such approaches are essential for the development and design of environmental conservation and climate adaptation policies (Wani et al. 2024b). Available models suggest that climate change may differentially influence the distribution ranges of various species. Patterns of expansion, stability, and retraction have been documented in plants such as Stevia Cav. and Hedeoma multiflora Benth. in mountain regions of north-western and central Argentina (Rodríguez-Cravero et al. 2017; Peralta et al. 2024) and Valeriana carnosa in Patagonia (Nagahama and Bonino 2020). Conversely, expansion patterns were observed in other medicinal Lippia plants, such as the shrub L. graveolens Kunth in the arid regions of Mexico (Martínez-Sifuentes et al. 2022). These patterns reveal the complexity in the response of native vegetation to climate change and underscore the importance of considering species-specific responses to climate changes. Furthermore, other factors such as interspecific competition, seed dispersal ability, topography, geology, and vegetation and land use could influence species distribution in the coming decades (Tsiftsis et al. 2024). Consequently, future studies must assess the importance of these factors in MAPs distribution to propose accurate conservation strategies.

In conclusion, these findings can drive local and regional efforts to address climate change and protect medicinal and aromatic flora, also supporting the sustenance of different regional economies. They also provide a strong basis for researching Verbenaceae conservation and guiding resource management and biodiversity protection.