Climatic suitability assessment of pistachio (Pistacia vera L.) under current and future climate scenarios in the Mediterranean Region

1 Introduction

Climate change is driving profound and accelerating transformations in agricultural systems worldwide, with the Mediterranean Basin emerging as a highly vulnerable region (Ali et al., 2022; Buchanan et al., 2025). Mediterranean surface temperature is now approximately 1.5 °C above pre-industrial levels, and warming is projected to continue faster than the global average (MedECC, 2020). Annual mean warming over Mediterranean land is projected to reach 0.9–5.6 °C by the end of the century, depending on the emission scenario, while precipitation is projected to decrease by 4.0–22.0% across most of the region (MedECC, 2020; Ali et al., 2022). These compound trends—rising heat, reduced and more variable rainfall, and more frequent drought—are intensifying climatic constraints in areas characterised by high seasonality and water deficits (Noto et al., 2023; Granata et al., 2025).

Perennial cropping systems are especially exposed because their productivity depends on stable climatic conditions across multiple years and on the fulfilment of phenological requirements through distinct growth stages (Meza et al., 2023; Caspersen et al., 2025). The Mediterranean Basin is characterised by strong latitudinal and longitudinal environmental gradients, complex interactions among temperature, precipitation, and topography, and diverse agroecological conditions (Lionello et al., 2014; Materia et al., 2022). This heterogeneity, combined with the basin’s status as a climate-change hotspot (MedECC, 2020; Ali et al., 2022), makes it a valuable region for examining how climate-driven shifts affect agricultural systems and reshape the spatial boundaries of viable cultivation zones (Costanza et al., 2024; Osorio-Marín et al., 2024). Pistachio (Pistacia vera L.) is a perennial nut crop with considerable tolerance to drought and high temperatures in semi-arid and transitional environments (Osku et al., 2025). It is widely cultivated because of its economic value and nutritional importance (Khadivi et al., 2025). Pistachio phenology is closely associated with temperature regimes: budbreak and flowering require sufficient winter chilling, whereas fruit development requires adequate heat accumulation during the growing season (Benmoussa et al., 2017; Alipour et al., 2025). Precipitation patterns, particularly during dry periods, also affect tree vigour and long-term persistence (Marino et al., 2018). Consequently, pistachio suitability reflects a complex trade-off among mean temperature, seasonal variability, and water availability, making the species sensitive to long-term climatic trends and interannual variability.

Historically, pistachio production has been concentrated in the Mediterranean Basin, the Middle East, and Central Asia, where climatic conditions closely match the ecological requirements of the species (Mir-Makhamad et al., 2022; Khadivi et al., 2025). Climate change is disrupting these established geographical patterns of cultivation (Doğan and Tekeş Düdükçü, 2025). To support strategic agricultural decisions and reduce climate-related risks, it is important to identify emerging suitable zones. During the past decade, land-suitability assessment and species distribution modelling have become important tools for analysing crop-climate interactions and projecting spatial changes in suitability (Li et al., 2023; Costanza et al., 2024; Çelik et al., 2025). These methods enable environmental data to be analysed at large spatial scales and allow the relative influence of numerous climatic and landscape factors to be quantified.

Parametric and index-based frameworks, such as the FAO land-evaluation system, rely on predefined thresholds and weighted agronomic criteria, which can limit their capacity to represent nonlinear environmental responses (Navidi et al., 2022). Multi-criteria evaluation methods combining the analytic hierarchy process with weighted overlay techniques allow flexible spatial integration but remain sensitive to expert-assigned priorities and often omit uncertainty treatment (Talukdar et al., 2022). MaxEnt is suitable for this application because it combines presence records with background environmental samples to estimate relative environmental suitability, accommodates nonlinear crop–environment relationships, and does not require confirmed absence records (Merow et al., 2013; Campos et al., 2023).

Several pistachio suitability studies are geographically limited to local or national case studies and therefore do not represent the full range of environmental gradients in the Mediterranean Basin (Mercan and Acibuca, 2023; Çelik et al., 2025). Methodological inconsistency across studies also complicates comparisons, particularly when climate variables are included without assessing redundancy or collinearity (Dormann et al., 2013). Existing research provides limited basin-scale evidence on emerging future cultivation zones under alternative climate scenarios, while projections based on individual global climate models (GCMs) do not represent inter-model climate uncertainty (Belvisi and TranChina, 2026). Consequently, a gap remains in identifying robust patterns of future pistachio climatic suitability across the Mediterranean Basin using multi-GCM ensembles. This study assesses current and future climatic suitability for pistachio cultivation across the Mediterranean Basin under intermediate (SSP2-4.5) and high (SSP3-7.0) climate change scenarios. The specific objectives include: (i) mapping and quantifying the current geographical distribution and extent of climatic suitability for pistachio cultivation; (ii) projecting and quantifying changes in climatic suitability using a five-GCM ensemble under SSP2-4.5 and SSP3-7.0 for 2021–2040, 2041–2060, and 2081–2100; and (iii) identifying potential newly suitable areas and comparing their spatial and elevational characteristics with currently suitable areas. By addressing these objectives, the study contributes to understanding the adaptation of perennial crops to climate change and supports strategic land-use planning in Mediterranean agroecosystems.

2 Materials and methods

2.1 Study area

The study encompassed the Mediterranean region, a major agroecological area for perennial crops and a climate change hotspot in which pistachio may provide an adaptation option. The Mediterranean Basin includes southern Europe, North Africa, and parts of the Middle East, and contains heterogeneous environments shaped by strong latitudinal and longitudinal gradients (Lionello et al., 2014). This region is characterised by hot, dry summers, mild to cool and relatively humid winters, and pronounced seasonal changes in temperature and rainfall (Todaro et al., 2022). These climatic features support diverse perennial crops, including pistachio, olive, almond, and grapevine, whose productivity depends on crop-specific thermal requirements, such as winter chilling and growing-season heat accumulation, as well as seasonal water availability (Materia et al., 2022). Pistachio is already cultivated widely in Mediterranean semi-arid environments because it tolerates prolonged dry periods, although yield remains dependent on adequate winter chilling, summer heat, and pollination conditions (Alipour et al., 2025; Burgos-Miranda et al., 2025). Climate change threatens Mediterranean agriculture by altering temperatures, rainfall, and seasonal variability across the basin (Ali et al., 2022; Noto et al., 2023). All spatial datasets were harmonised and clipped to this domain (Figure 1).

2.2 MaxEnt modelling framework

MaxEnt was selected because the GBIF dataset contained presence-only records without reliable absences; locations without records may reflect incomplete or uneven sampling rather than unsuitable conditions. MaxEnt estimates relative environmental suitability by contrasting conditions at occurrence locations with those available across a defined background and identifies the maximum-entropy distribution subject to constraints derived from occurrence-site predictors (Elith et al., 2011; Merow et al., 2013). Regularisation helps limit overfitting, while linear and nonlinear feature classes accommodate complex crop–environment relationships across heterogeneous climatic gradients (Merow et al., 2013; Phillips et al., 2017; Kass et al., 2021). The background represented the accessible environmental space within the Mediterranean study region, allowing relative suitability to be estimated without treating unrecorded locations as true absences. The MaxEnt model was implemented using maxnet through ENMeval, tuned and evaluated through four-fold spatial cross-validation before projection across 30 combinations of five GCMs, two SSPs, and three future periods. Five-GCM ensemble predictions and extrapolation diagnostics were subsequently used to quantify future suitability and spatial changes (Figure 2).

2.3 Species occurrence data

Occurrence records of Pistacia vera L. were downloaded from the Global Biodiversity Information Facility (GBIF) on 2 January 2026 using a species-level taxonomic query (GBIF.org, 2026). The download contained 4,331 records. Records outside the Mediterranean modelling region and duplicate coordinate records were removed, leaving 1,252 records. Removal of 607 records with empty or invalid coordinates left 645 records. Fourteen records with reported coordinate uncertainty exceeding 5 km were excluded; 202 records without reported uncertainty were retained after coordinate validation, producing 631 cleaned occurrences. Extraction against the common environmental mask excluded one boundary record lacking complete predictor data, leaving 630 occurrences for model calibration and evaluation (Figure 1). Spatial non-independence was addressed through four-fold spatial cross-validation rather than distance-based thinning.

2.4 Selection and screening of environmental indicators

Nineteen bioclimatic variables representing the 1970–2000 baseline were obtained from WorldClim v2.1 at 2.5-arc-minute resolution (Fick and Hijmans, 2017). Slope was derived in degrees from the WorldClim elevation layer at 30-arc-second resolution and subsequently aggregated by the mean to the 2.5-arc-minute grid. The 19 bioclimatic variables and slope formed an initial set of 20 candidate predictors (Table 1). All layers were clipped to the Mediterranean modelling region and harmonised to a common WGS 84 coordinate system, spatial extent, resolution, origin, grid alignment, and valid-cell mask. Predictor screening was conducted using 10,000 randomly sampled cells from the Mediterranean environmental background. Pairwise Pearson correlation identified highly correlated predictor groups using |r| ≥ 0.70 (Figure 3), while principal component analysis of standardised predictor values was used diagnostically to characterise the dominant environmental gradients. PC1 and PC2 explained 61.72% and 13.91% of the variance, respectively, accounting for 75.63% cumulatively (Supplementary Figure S1). Predictor retention combined correlation structure, ecological relevance, future-data availability, and variance inflation factors, with VIF < 5 considered acceptable. This process retained BIO1, BIO12, BIO15, and slope as candidate predictors. Although BIO1 and BIO12 were correlated (r = −0.844), their VIF values in the four-predictor candidate set were 3.784 and 3.618, respectively; therefore, both were retained for subsequent model comparison. Alternative predictor sets were evaluated during model calibration to test the effects of temperature, precipitation, and slope independently. Detailed PCA contributions, VIF values, and predictor-selection outcomes are provided in Supplementary Table S1.

Table 1. Candidate bioclimatic and topographic predictors considered before variable screening for pistachio suitability assessment.

BIO, bioclimatic variable. BIO1-BIO19 were obtained from WorldClim v2.1 for 1970–2000 at 2.5-arc-minute resolution (Fick and Hijmans, 2017). Slope was calculated in degrees fromnthe 30-arc-second WorldClim elevation layer and aggregated by the mean to 2.5-arc-minute resolution.

2.5 Calibration and configuration of the MaxEnt model

The four candidate predictor sets and 630 presence records were used to fit candidate MaxEnt models using the maxnet implementation through ENMeval in R v4.4.2 (Phillips et al., 2017; Kass et al., 2021). Four predictor sets were evaluated: Set A (BIO1, BIO12, BIO15 and slope), Set B (BIO1, BIO15 and slope), Set C (BIO12, BIO15 and slope), and Set D (BIO1, BIO12 and BIO15). Each predictor set was evaluated using linear (L), linear–quadratic (LQ), and linear–quadratic–hinge (LQH) feature-class combinations and regularisation multipliers (RM) from 0.5 to 4.0 at 0.5 intervals, producing 96 candidate models. Models used 10,000 randomly sampled background cells from the Mediterranean modelling region and complementary log–log (cloglog) output for suitability prediction. Four spatial folds were generated by k-means clustering of the occurrence coordinates. Each background point was assigned to the nearest occurrence-cluster centroid so that occurrence and background records followed the same spatial partitions. Each fold was withheld once for validation. Candidate models were compared using the Akaike information criterion corrected for small sample size (AICc), validation area under the receiver operating characteristic curve (AUC), 10th-percentile training-presence omission rate, validation continuous Boyce index (CBI), and number of non-zero parameters. Although the Set D model with LQH features and RM = 0.5 produced the minimum AICc, RM = 1.0 was retained as a compromise between model fit and parsimony because it reduced the number of non-zero parameters from 23 to 17, slightly reduced omission from 0.200 to 0.192, increased validation CBI from 0.735 to 0.762, and maintained equivalent validation AUC (0.892). Models with RM values of 1.5 or greater increased CBI and reduced omission rates but incurred AICc penalties of at least 26.3 relative to RM = 1.0 (37.0 at RM = 2.0), supporting retention of RM = 1.0 as a balance between fit and complexity. At RM = 1.0, the slope-containing Set A had an AICc only 0.45 lower than Set D but a slightly lower validation AUC (0.889 versus 0.892); therefore, slope was excluded by parsimony. The final model comprised BIO1, BIO12, and BIO15 with LQH features and RM = 1.0 and was refitted using all 630 occurrences and 10,000 background cells.

2.6 Evaluation of model performance and predictor importance

Model performance was evaluated using four-fold spatial cross-validation. For each fold, the model was trained with occurrences and background records from the other three folds and evaluated with the withheld records. Performance was assessed using the area under the receiver operating characteristic curve (AUC), the continuous Boyce index (CBI), and the 10th-percentile training-presence omission rate. The partitioned-background AUC was used as the primary measure of spatial predictive performance. Evaluation metrics were summarised as the mean and standard deviation across the four folds. Predictor importance was evaluated separately within each withheld fold. Values of each predictor in the withheld evaluation data were randomly permuted 30 times, and importance was calculated from the resulting reduction in AUC relative to the unpermuted model and normalised to sum to 100%. Jackknife analysis was performed by refitting models with each predictor used alone and with each predictor omitted; model performance was compared using the training deviance ratio (dev.ratio). Marginal response curves were generated for each fold by varying one predictor across a common sequence of 100 values while holding the other predictors at their training-data medians. The mean response and standard deviation across the four folds were then calculated to represent the fitted relationship and its uncertainty.

2.7 Future climate data and ensemble projections

Future suitability was projected using WorldClim v2.1 bias-corrected and downscaled CMIP6 bioclimatic data at 2.5-arc-minute resolution (Eyring et al., 2016; WorldClim, 2025). Five global climate models (GCMs)—ACCESS-CM2, INM-CM5-0, IPSL-CM6A-LR, MPI-ESM1-2-HR, and MRI-ESM2-0—were evaluated under SSP2-4.5 and SSP3-7.0 (O’Neill et al., 2016) for 2021–2040, 2041–2060, and 2081–2100, representing near-, mid-, and late-century conditions, respectively. SSP2-4.5 and SSP3-7.0 were selected to contrast intermediate and high forcing without relying on the more extreme SSP5-8.5 pathway. This design produced 30 GCM–scenario–period combinations. The three predictors retained in the final model (BIO1, BIO12, and BIO15) were extracted from each precomputed bioclimatic dataset, clipped to the Mediterranean modelling region, and aligned with the current predictor grid in coordinate reference system, resolution, extent, origin, and valid-cell mask. The calibrated maxnet model was projected to each combination using complementary log–log output with clamping enabled, while differences between clamped and unclamped predictions were retained as diagnostics. For each scenario and period, the five GCM predictions were summarised cellwise using the median as the ensemble estimate; the interquartile range, standard deviation, minimum, and maximum described inter-model variation. Multivariate environmental similarity surface (MESS) analysis compared each future predictor stack with current environmental conditions, with negative values identifying cells containing at least one predictor outside its current range. Using the 10th-percentile threshold defined below, GCM agreement was calculated as the number of models (0–5) predicting suitability in each cell, and agreement among at least three models defined majority suitability. Novelty frequency was calculated as the number of GCMs (0–5) with negative MESS values in each cell. The climatic spread across GCMs was characterised using geodesic-area-weighted changes in BIO1, BIO12, and BIO15 relative to the 1970–2000 baseline (Supplementary Figure S2).

2.8 Suitability classification and change analysis

Continuous cloglog predictions for current conditions and the six future ensemble medians were classified into five suitability classes: very low (0.0–0.2), low (>0.2–0.4), moderate (>0.4–0.6), high (>0.6–0.8), and very high (>0.8–1.0). For binary analyses, current, individual GCM, and ensemble-median predictions were classified using the 10th-percentile training-presence threshold from the final full model (cloglog = 0.3371397935). Cross-tabulation of current and future ensemble-median binary maps produced four transition classes: stable unsuitable, suitability loss, persistence, and suitability gain. Geodesic cell areas were calculated on the WGS 84 grid using terra::cellSize(unit = “km”) and summed by class; the total modelling domain was 8,469,574 km². Elevation was extracted from the WorldClim v2.1 2.5-arc-minute elevation raster after alignment to the model grid. Geodesic-area-weighted mean, standard deviation, and median elevation were calculated for currently suitable cells and newly suitable (gain) cells to characterise altitudinal redistribution.

3 Results

3.1 Model discrimination and spatial reliability

The final model (Set D; LQH features; RM = 1.0) had 17 non-zero parameters and an AICc of 14,576.9 (ΔAICc = 16.6 from the AICc-minimum RM = 0.5 configuration; Supplementary Table S2). The 10th-percentile training-presence omission rate averaged 0.192 ± 0.168 and ranged from 0.034 to 0.420 among folds. Fold 1 comprised 100 records, mainly from France (61) and Italy (26), with smaller groups from Spain and Tunisia (six each) and Monaco (one); its omission rate was 0.420 compared with the nominal 0.10 reference (Figure 4). The final model had a training AUC of 0.924. Under the primary partitioned-background evaluation, mean AUC was 0.886 ± 0.034, with fold-specific values ranging from 0.840 to 0.917 (Figure 5); the corresponding ENMeval full-background validation AUC was 0.892 ± 0.021. Validation CBI averaged 0.762 ± 0.134 and ranged from 0.580 to 0.890. Notably, Fold 2 combined the highest partitioned-background AUC (0.917) with the lowest CBI (0.580), showing that the two metrics captured complementary aspects of model evaluation. Thus, discrimination remained above random expectation in every fold, but the pronounced heterogeneity in omission indicates spatial variation in model transferability across the Mediterranean modelling region.

3.2 Environmental variable importance and response patterns

Spatial-fold permutation tests identified annual precipitation (BIO12) as the strongest predictor, accounting for 69.8% of total permutation importance (mean held-out AUC reduction = 0.213 ± 0.031 across folds). Annual mean temperature (BIO1) accounted for 28.7% (0.087 ± 0.064), whereas precipitation seasonality (BIO15) accounted for 1.5% (0.005 ± 0.010). The BIO15 AUC reduction was small, highly variable, and not consistently positive among folds. Normalised absolute coefficient sums were highest for BIO1 (43.5%), followed by BIO12 (29.0%) and BIO15 (27.5%). The contrast between these diagnostics shows that coefficient magnitude did not translate directly into independent predictive information, particularly for BIO15 (Figure 6).

Jackknife refits showed that BIO1 produced the greatest training deviance ratio when used alone (0.151), followed by BIO12 (0.136), whereas BIO15 alone produced little explanatory gain (0.010). Relative to the all-variable model (0.155), omitting BIO1 reduced the deviance ratio to 0.127 and omitting BIO12 reduced it to 0.133; omission of BIO15 produced only a small reduction to 0.151. The contrasting rankings arise because the diagnostics quantify different properties: permutation importance measures the loss of held-out discrimination after a predictor is disrupted, whereas jackknife deviance ratios measure training fit when a predictor is used alone or omitted. Accordingly, BIO12 contributed more uniquely to held-out discrimination, while BIO1 provided stronger standalone training fit and caused a larger reduction when omitted (Figure 7).

Marginal response curves, with the other predictors held at their median values, showed a broad maximum in suitability at annual mean temperatures of approximately 15.1–16.6 °C, with suitability declining towards both cooler and warmer conditions. Suitability increased sharply with annual precipitation to a maximum near 405 mm and then declined progressively at higher precipitation. The response to BIO15 (coefficient of variation, expressed as %) was non-monotonic, with broad maxima near 38% and 92%. Across-fold variability was lowest for BIO1 (mean SD = 0.058; maximum SD = 0.158) and markedly greater for BIO12 (mean SD = 0.154; maximum SD = 0.421) and BIO15 (mean SD = 0.162; maximum SD = 0.395). Because BIO1 and BIO12 were strongly correlated (r = −0.844 in the predictor-screening sample), their one-dimensional marginal sweeps may include climatic combinations that are sparse or absent in the current Mediterranean domain. Their joint fitted response, together with current environmental support and occurrence environments, is shown in Supplementary Figure S3. The marginal curves therefore describe conditional model responses and should not be interpreted as causal physiological thresholds or as the complete fitted surface where predictors are strongly correlated (Figure 8).

3.3 Current pattern of climatic suitability for Pistacia vera in the Mediterranean Basin

Under current climatic conditions, suitability was strongly skewed towards the lowest class across the 8,469,574 km² modelling domain (Figures 9A, B). Very low suitability (cloglog 0.0–0.2) covered 6,595,540 km² (77.9%), while low (>0.2–0.4) and moderate (>0.4–0.6) suitability occupied 555,779 km² (6.6%) and 484,980 km² (5.7%), respectively. High (>0.6–0.8) and very high (>0.8–1.0) suitability covered 637,877 km² (7.5%) and 195,397 km² (2.3%), together representing 833,274 km² (9.8%). Using the 10th-percentile training-presence threshold (cloglog = 0.337), 1,464,982 km² (17.3%) was classified as suitable. This larger estimate includes cells above the occurrence-derived threshold but below the 0.4 class boundary.

High and very high suitability formed discontinuous concentrations in interior Iberia, the Tell Atlas and northern Morocco, interior Sicily and southern Italy, central Greece and the Aegean, and western, central, and south-eastern Anatolia, with smaller patches along western Syria and Lebanon (Figure 9A). Moderate suitability generally formed transitional fringes around these areas. Very low suitability prevailed in both the hyper-arid Sahara and Syro-Arabian interior and the cooler, humid regions of northern France, the Alpine arc, and Atlantic northern Iberia. The climatic distribution of the domain was consistent with these mapped contrasts. Area-weighted median annual precipitation was 55 mm, compared with a modelled suitability maximum near 405 mm, while median annual mean temperature was 20.6 °C, above the response maximum of approximately 15.1–16.6 °C. At the 630 calibration locations, median predicted suitability was 0.715, and 567 records (90.0%) exceeded the 0.337 threshold.

3.4 Projected climatic suitability under future scenarios

The ensemble median projected an increase in threshold-defined suitable area under both pathways (Figure 10, Table 2). Relative to the current 1,464,982 km² (17.3% of the domain), suitable area under SSP2-4.5 reached 1,809,591 km² (21.4%) in 2021–2040, 1,865,053 km² (22.0%) in 2041–2060, and 1,991,116 km² (23.5%) in 2081–2100. Corresponding values under SSP3-7.0 were 1,777,821 km² (21.0%), 1,910,075 km² (22.6%), and 1,979,047 km² (23.4%). Differences between pathways did not exceed 0.6 percentage points at any horizon and were smaller than the 1.9–3.1 percentage-point ranges among individual GCMs. By 2081–2100, very low suitability had declined from 77.9% under current conditions to 70.9% under SSP2-4.5 and 70.5% under SSP3-7.0.

Table 2. Ensemble-median areas of projected climatic suitability classes and binary suitability transitions under SSP2-4.5 and SSP3-7.0.

Most of the increase occurred in the low and moderate classes. By 2081–2100, moderate suitability increased from 484,980 km² to 853,615 km² under SSP2-4.5 and 846,271 km² under SSP3-7.0, while low suitability increased from 555,779 km² to 761,349 and 775,870 km², respectively. Across all future combinations, high suitability varied between 613,668 and 652,985 km², close to its current extent of 637,877 km², and very high suitability ranged from 226,274 to 249,177 km². The expansion therefore consisted mainly of transitional suitability rather than a marked enlargement of the two highest classes. New moderate and high suitability was concentrated in southern and central France, northern Iberia, the Po basin and northern Adriatic, the Balkans, and central and northern Anatolia. High suitability contracted along parts of the Tell Atlas and Levantine corridor, whereas interior Iberia, Sicily and southern Italy, central Greece, and south-eastern Türkiye remained favourable across all periods (Figure 10).

Inter-GCM spread generally increased with projection horizon. Mean ensemble IQR was 0.014–0.016 cloglog units in 2021–2040 and reached 0.035 under SSP2-4.5 and 0.044 under SSP3-7.0 in 2081–2100; corresponding local maxima were 0.348 and 0.469. Areas classified as suitable by all five GCMs represented 88.2–90.2% of ensemble suitable area in 2021–2040, declining to 71.6% under SSP2-4.5 and 67.3% under SSP3-7.0 by 2081–2100. Negative MESS values affected 0.8–4.2% of the domain in 2021–2040 and 12.0–26.3% under SSP3-7.0 in 2081–2100, depending on the GCM. Across all 30 projections, negative-MESS cells did not overlap areas classified as suitable at the 0.337 threshold, indicating limited extrapolation within projected suitable areas.

3.5 Spatial dynamics of suitability change

Cross-tabulation of current and future binary surfaces separated change into suitability loss, persistence, gain, and stable unsuitable area (Figure 11, Table 2). Because the current binary surface was fixed, persistence plus loss equalled the current suitable extent of 1,464,982 km² in every comparison. Under SSP2-4.5, persistent suitable area declined from 1,337,368 km² in 2021–2040 to 1,271,037 km² in 2041–2060 and 1,133,136 km² in 2081–2100, retaining 91.3%, 86.8%, and 77.3% of currently suitable area. Loss increased from 127,614 to 193,945 and 331,846 km², while gain increased from 472,223 to 594,015 and 857,980 km². Under SSP3-7.0, persistence declined from 1,344,045 to 1,236,627 and 894,382 km²; by 2081–2100, loss reached 570,600 km² (38.9% of currently suitable area) and gain reached 1,084,665 km² (74.0%). Gains exceeded losses in every comparison. Net increases ranged from 344,609 to 526,134 km² under SSP2-4.5 (23.5–35.9% of current suitable area) and from 312,839 to 514,064 km² under SSP3-7.0 (21.4–35.1%). Gains occurred mainly across France, northern Iberia, the Po basin, the Balkans, and central and northern Anatolia, whereas losses were concentrated along the Atlas belt and the Levantine corridor (Figure 11). Persistent suitability remained most extensive in interior Iberia, Sicily and southern Italy, central Greece, and central and south-eastern Türkiye. By 2081–2100 under SSP3-7.0, losses extended across much of the North African suitable belt, demonstrating that the net basin-wide increase masked substantial regional turnover.

Gain areas also occurred at higher elevations than the current suitable area (Table 3). Currently suitable cells had an area-weighted mean elevation of 510 ± 382 m (median 453 m). Mean elevations of gain areas were 693 ± 501, 728 ± 522, and 667 ± 533 m under SSP2-4.5 and 712 ± 496, 691 ± 523, and 747 ± 580 m under SSP3-7.0 for 2021–2040, 2041–2060, and 2081–2100, respectively. These means were 157–237 m above the current value, although the large standard deviations indicate substantial overlap among elevation distributions. The projected increase in total suitable area therefore coincided with redistribution towards higher latitudes and elevations and with losses along the southern basin margin.

Table 3. Elevational characteristics of currently suitable and projected suitability-gain areas under SSP2-4.5 and SSP3-7.0.

Means, standard deviations (SD), and medians were weighted by geodesic cell area. Future rows describe cells classified as suitability gain using the 10th-percentile training-presence threshold (cloglog = 0.337). Differences are relative to the current mean elevation (510 m).

4 Discussion

4.1 Model performance and reliability of the projections

The final model retained useful discriminatory and ranking ability under spatial cross-validation, with a partitioned-background AUC of 0.886 ± 0.034 and a validation continuous Boyce index (CBI) of 0.762 ± 0.134. In presence–background models, AUC measures separation of occurrences from the selected background and is sensitive to the extent and composition of that background. CBI provides a complementary test of whether withheld occurrences are concentrated in cells assigned higher suitability (Heumann et al., 2011; Lissovsky and Dudov, 2021). The retained regularisation multiplier (RM = 1.0) was not the AICc minimum, but it preserved the validation AUC of the RM = 0.5 model while reducing the number of non-zero parameters from 23 to 17, slightly lowering omission and increasing CBI. Its selection therefore represents a stated compromise between fit and parsimony, evaluated with spatial rather than random partitions (Phillips et al., 2017; Kass et al., 2021).

Performance varied substantially among folds. Fold 1, which contained 100 records mainly from France and Italy, had an omission rate of 0.420, compared with the nominal 0.10 reference. Fold 2 combined the highest partitioned-background AUC (0.917) with the lowest CBI (0.580). These results show that discrimination and suitability ranking were not spatially uniform and suggest weaker transfer to parts of the northern calibration margin. This point is relevant because several future gains are projected in the same broad region; local predictions there require greater caution than the basin-wide mean metrics imply.

The extrapolation analysis provides a separate qualification. Negative MESS values affected 0.8–4.2% of the domain in 2021–2040 and 12.0–26.3% by 2081–2100 under SSP3-7.0, depending on the GCM, but did not overlap threshold-defined suitable cells in any of the 30 projections. Future suitable areas therefore did not rely on predictor values outside the univariate ranges used for calibration. MESS does not, however, exclude every novel combination among correlated predictors, so this result should not be read as proof of complete climatic analogy (Merow et al., 2013). Agreement among all five GCMs declined from 88.2–90.2% of ensemble suitable area in 2021–2040 to 67.3% under SSP3-7.0 by 2081–2100. Late-century spatial detail and magnitude are consequently less certain than the broad direction of redistribution.

4.2 Climatic controls on pistachio distribution

Annual precipitation (BIO12) and annual mean temperature (BIO1) accounted for 69.8% and 28.7% of permutation importance, respectively, whereas precipitation seasonality (BIO15) accounted for 1.5% and produced an AUC change close to zero across folds. The jackknife ranked BIO1 slightly above BIO12 when each was fitted alone, while permutation ranked BIO12 first in held-out discrimination. These diagnostics measure different properties and together indicate that thermal and moisture regimes carry most of the model information. Their individual percentages should not be interpreted as causal shares because BIO1 and BIO12 were strongly correlated (r = −0.844).

The precipitation response peaked near 405 mm, consistent with the occurrence of pistachio in dry Mediterranean environments and with the importance of water supply for commercial production (Marino et al., 2018). The decline in fitted suitability above this value describes the climatic distribution of recorded cultivation; it does not imply that reduced rainfall improves orchard performance. Irrigation remains important where rainfall is insufficient, and a process-based analysis in a Mediterranean-climate production region projected an increase in pistachio crop evapotranspiration under warming (Kandhway et al., 2025).

BIO1 showed a broad suitability maximum at approximately 15.1–16.6 °C, but annual mean temperature cannot resolve the sequential thermal requirements of pistachio. Winter chill regulates dormancy release and flowering, whereas growing-season heat controls fruit and kernel development. Warm winters can delay or disrupt budbreak and reduce yield, and mismatched male and female flowering can limit pollination (Benmoussa et al., 2017; Elloumi et al., 2024; Burgos-Miranda et al., 2025). At the cooler margin, adequate heat accumulation is required to complete nut development, with marked differences among cultivars (Núñez et al., 2024; Álvarez et al., 2026). Because chill portions and growing-degree measures were not included, these mechanisms are represented only indirectly by BIO1.

4.3 Current suitability and correspondence with production geography

Under current conditions, high and very high suitability occupied 9.8% of the modelling domain, and 17.3% exceeded the 10th-percentile training-presence threshold. The principal suitable zones in interior Iberia, Sicily and southern Italy, central Greece and the Aegean, and western to south-eastern Anatolia broadly coincide with established pistachio-growing areas. This pattern is consistent with regional land-suitability studies in Türkiye (Mercan and Acibuca, 2023; Çelik et al., 2025). Predicted suitability was also high at the calibration locations: the median was 0.715 and 90.0% of records exceeded the threshold. This agreement is a useful check on the fitted surface but is not independent validation because the occurrence data contributed to model calibration.

The interpretation is further constrained by the status of Pistacia vera as a cultivated crop across most of the study region. GBIF records represent a cultivation footprint shaped by human dispersal, irrigation, land use and market history as well as climate. The species was domesticated in Central Asia and dispersed westwards through long-term cultivation, whereas wild Pistacia vera communities are centred in Middle and south-western Asia (Mir-Makhamad et al., 2022; Nowak et al., 2024). The map therefore identifies climates resembling those of recorded cultivation sites; it should not be interpreted as a direct estimate of physiological tolerance or expected yield.

4.4 Projected redistribution: expansion at the margin, erosion at the core

Threshold-defined suitable area increased from 17.3% at present to 23.5% under SSP2-4.5 and 23.4% under SSP3-7.0 by 2081–2100. This net increase was concentrated in marginal and transitional classes. Moderate suitability expanded by 76.0% and 74.5% under the two pathways, whereas high suitability remained between 613,668 and 652,985 km², no more than 3.8% below or 2.4% above its current extent. Very high suitability increased from 195,397 km² to 226,274–249,177 km² but remained a small share of the domain. The result is therefore better described as an expansion of the climatic margin than enlargement of the climatic optimum. Expansion was accompanied by substantial turnover. By 2081–2100, 22.7% of currently suitable land was lost under SSP2-4.5 and 38.9% under SSP3-7.0. Gains were concentrated across France, northern Iberia, the Po basin, the Balkans, and central and northern Anatolia, while losses occurred mainly along the Atlas belt and the Levantine corridor. Gain areas were also located at higher elevations than the current suitable distribution, with mean elevation differences of 157–237 m. Comparable reorganisation towards northern marginal areas has been projected for Mediterranean olive, and temperature-related shifts are widely reported in temperate fruit and nut systems (Meza et al., 2023; Osorio-Marín et al., 2024; Belvisi and TranChina, 2026). Climate projections for Pistacia atlantica in the Zagros likewise show substantial redistribution and contraction of the highest-suitability class (Zadeh et al., 2024).

The concentration of losses along the southern and eastern margins agrees with the broader Mediterranean climate signal, where warming exceeds the global mean and drying is strongest across many southern sectors (MedECC, 2020; Ali et al., 2022). Increasing drought frequency and hydrological stress are also projected unevenly across the basin (Noto et al., 2023; Granata et al., 2025). The North African and Levantine losses are therefore consistent with combined thermal and moisture pressure. Nevertheless, the close late-century totals under SSP2-4.5 and SSP3-7.0 should not obscure their different spatial turnover and uncertainty.

4.5 Implications for adaptation and land-use planning

The northern gains identify locations for further agronomic assessment rather than land immediately available for orchard development. Field observations from Castilla y León show that chilling and heat requirements can be met in a cold-winter Spanish region, although cultivar choice affects yield, blank nuts and shell splitting (Núñez et al., 2024). Measurements from central Spain also demonstrate cultivar-specific heat requirements for kernel development (Álvarez et al., 2026). Pistachio has been considered as an alternative crop in Mediterranean olive areas affected by Xylella fastidiosa (Greco et al., 2024). Before investment, however, climatic screening must be followed by assessment of soils, water availability, spring-frost exposure, polliniser synchrony and orchard establishment constraints.

In established southern production areas, the priority is to retain productivity where climate suitability declines. Low-chill cultivars may reduce exposure to warmer winters, while rootstock selection offers scope to improve tolerance to drought and salinity (Akbari et al., 2024; Alipour et al., 2025; Osku et al., 2025). Plant-based monitoring of stem water potential can support more precise irrigation decisions where water becomes limiting (Mortazavi et al., 2025). These measures complement, rather than replace, suitability mapping. Within climate-smart planning, the maps can be used to prioritise field evaluation, guide cultivar and rootstock trials, and identify areas where long-lived orchard investments face increasing climatic risk (Buchanan et al., 2025).

4.6 Limitations and future research

The principal limitation is the absence of explicit winter-chill and growing-season heat predictors. BIO1 cannot distinguish warming that alleviates a heat deficit near the northern margin from warming that reduces chill accumulation in established orchards. Future calibration should include chill portions and cultivar-relevant growing-degree measures and test whether they alter the projected balance between northern gains and southern losses (Benmoussa et al., 2017; Elloumi et al., 2024; Núñez et al., 2024). Process-based phenology and water-balance modelling would also provide a mechanistic test of the correlative temperature and precipitation responses (Ahmad et al., 2025; Kandhway et al., 2025). The occurrence data introduce a second limitation. They are presence-only records of a cultivated species and may retain spatial variation in recording effort and orchard accessibility. Spatial cross-validation tests geographic transfer more strictly than random partitioning, but it does not remove bias in the underlying observation process. Independent orchard inventories, yield records and planned field surveys are needed to test the predictions. Correlation between BIO1 and BIO12 also limits causal interpretation of their marginal responses, even though VIF values within the final set were below five.

The model represents climate alone. Soil properties, irrigation access, salinity, cultivar, rootstock and socio-economic conditions determine whether climatic suitability becomes viable production (Mercan and Acibuca, 2023; Seyedmohammadi et al., 2023). Elevation was used as a post-hoc descriptor of gain areas, not as a retained model predictor; the observed upslope shift therefore reflects the geography of projected climate rather than an independently estimated topographic effect. Genetic improvement may also change the realised tolerance envelope, reinforcing the value of conserving Pistacia genetic resources (Maxted et al., 2010).

Finally, the ensemble samples five GCMs and two forcing pathways but retains one modelling algorithm and assumes that the fitted cultivation–climate relationship remains stable. Declining GCM agreement and increasing novelty show that uncertainty grows towards the end of the century. MESS identifies values outside the calibration range but does not fully diagnose novel combinations of correlated predictors. Further work should compare alternative modelling algorithms, apply complementary extrapolation diagnostics, and validate emerging northern and upland zones agronomically before the maps are used for site-level investment decisions.

5 Conclusion

This study provides a basin-scale, five-GCM assessment of climatic suitability for Pistacia vera across the Mediterranean Basin using a spatially cross-validated MaxEnt model. Under current conditions, high and very high suitability occupy 9.8% of the domain and 17.3% is suitable at the occurrence-derived threshold. Suitable areas are concentrated in interior Iberia, the Atlas belt, Sicily and southern Italy, central Greece, and Anatolia. Their broad correspondence with established cultivation zones supports the geographical plausibility of the model, although it is not independent validation because occurrence records informed calibration. Under SSP2-4.5 and SSP3-7.0, ensemble-median suitable area increases through the century, reaching 23.5% and 23.4% of the domain, respectively, by 2081–2100. Most of this increase occurs in low and moderate suitability, while high suitability remains close to its current extent. At the same time, 22.7–38.9% of currently suitable land is projected to be lost by late century. Gains occur mainly across France, northern Iberia, the Po basin, the Balkans, and central and northern Anatolia, at mean elevations 157–237 m above the current suitable distribution, whereas losses concentrate along the North African and Levantine margins. Net expansion therefore coincides with substantial geographical relocation and erosion of established southern suitable areas. These patterns require different responses across the basin. Adaptation in established southern regions should prioritise low-chill cultivars, drought- and salinity-tolerant rootstocks, and efficient water management. Northern and upland gains identify areas for further assessment rather than immediate orchard development because soil, water, frost and pollination constraints were not represented. Explicit winter-chill and growing-season heat metrics, independent orchard validation, and integration of edaphic and water-balance information are the main priorities for refining these projections and supporting climate-resilient pistachio planning.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fagro.2026.1965432/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Andrew Juan Challinor, University of Leeds, United Kingdom

Reviewed by: Milena Jovica Gocić, University of Nis, Serbia

Iftikhar Ali, University of Agriculture, Faisalabad, Pakistan