1 Introduction

Testudo graeca ibera, commonly referred to as the “spur-thighed tortoise” or “İber Kaplumbağası” in Turkish, is a terrestrial chelonian belonging to the Testudinidae family. [1] This subspecies is naturally distributed across the eastern Mediterranean region, including Anatolia, the Balkans, and parts of the Caucasus, where it occupies a wide range of Mediterranean and semi-arid habitats. [2] Its presence in these regions reflects long-term ecological adaptation rather than recent human-mediated introduction. Notably widespread in temperate and Mediterranean climate zones, Testudo graeca ibera exhibits considerable ecological plasticity, allowing it to persist under varying temperature and precipitation regimes. [3] The species commonly inhabits scrublands, open woodlands, grasslands, and agricultural landscapes, demonstrating tolerance to habitat fragmentation and seasonal climatic fluctuations. Due to its sensitivity to environmental change and its broad yet structured distribution, this subspecies represents an informative model organism for ecological niche modeling studies aimed at understanding climate–habitat relationships. # Methods

1.1 Data

1.1.1 Acquiring Species Occurrence Data

To acquire the occurrence data for Testudo graeca ibera, we used the Global Biodiversity Information Facility (GBIF) database, which provides globally curated species occurrence records. Occurrence data were obtained by querying the GBIF database for records corresponding to Testudo graeca ibera and filtering the results to include only observations with available geographic coordinates. The downloaded dataset was imported into R as a comma-separated values (CSV) file. To prepare the data for ecological niche modeling, only the geographic coordinate information—specifically the decimalLongitude and decimalLatitude columns—was retained. Records containing missing or incomplete coordinate values were removed using the na.omit function to ensure data quality and reliability. The resulting cleaned dataset consisted of georeferenced occurrence points representing the observed distribution of Testudo graeca ibera. These occurrence records formed the basis for subsequent ecological niche modeling analyses conducted in this study.

1.1.2 Formatting Data for biomod2

To format the occurrence data for ecological niche modeling using the biomod2 package, the dataset was prepared by combining species presence records of Testudo graeca ibera with environmental variables. The following steps outline the data preparation process:

  1. Preparing the Species Occurrence Data: The geographic coordinates (decimalLongitude and decimalLatitude) of Testudo graeca ibera presence records, obtained from GBIF, were formatted as a DataFrame. This DataFrame was then converted into a spatial object using the vect function from the terra package, with the Coordinate Reference System (CRS) defined as EPSG:4326.

  2. Preparing the Environmental Variables: Bioclimatic variables were downloaded from the WorldClim database at a spatial resolution of 10 minutes. The variables were stored as .tif files and loaded into R as a SpatRaster object using the rast function from the terra package. These variables were used as environmental predictors for modeling the distribution of Testudo graeca ibera.

    Selected bioclimatic variables include:

    • BIO1 Annual Mean Temperature
    • BIO5 Max Temperature of Warmest Month
    • BIO6 Min Temperature of Coldest Month
    • BIO12 Annual Precipitation
    • BIO15 Precipitation Seasonality
    • BIO17 Precipitation of Driest Quarter
  3. Integrating Species Data with Environmental Variables: Environmental variable values associated with the species occurrence locations were obtained using the extract function. This procedure produced a dataset in which geographic coordinates were linked with the corresponding bioclimatic values for each occurrence point.

  4. Formatting Data for Modeling: Occurrence records were prepared for analysis with the biomod2 package through the BIOMOD_FormatingData function. To obtain a balanced dataset, pseudo-absence locations were randomly sampled within the defined study area. Species presence records were coded as 1, whereas pseudo-absence records were assigned a value of 0. The main parameters applied during the data formatting process are outlined below:

    • resp.var: Species presence data.

    • expl.var: Environmental predictor variables (bioclimatic variables).

    • resp.xy: Coordinates of species presence points.

    • PA.nb.absences: Number of pseudo-absence points generated (set to 10,000).

This procedure prepared the dataset for species distribution modeling with the biomod2 package. By linking species presence records with environmental predictors and incorporating pseudo-absence points, a comprehensive dataset suitable for ecological niche modeling was obtained.

1.1.3 Acquiring Historical & Future Bioclimatic Variables

Historical and future bioclimatic datasets were obtained to characterize the ecological niche of Testudo graeca ibera and to assess its potential distribution under present and projected climate conditions. The data acquisition process followed the steps outlined below:

  1. Historical Bioclimatic Variables:

    • Past bioclimatic variables were obtained from the WorldClim database at a spatial resolution of 10 minutes.

    • The variables, provided in .tif format, were imported into R as a SpatRaster object using the rast function from the terra package.

    • The following bioclimatic variables were selected for this study:

    • BIO1 Annual Mean Temperature

    • BIO5 Max Temperature of Warmest Month

    • BIO6 Min Temperature of Coldest Month

    • BIO12 Annual Precipitation

    • BIO15 Precipitation Seasonality

    • BIO17 Precipitation of Driest Quarter

    • These variables were selected based on their established influence on habitat suitability and environmental tolerance in terrestrial species.

  1. Future Bioclimatic Variables:

    • Future bioclimatic variables were obtained from the WorldClim database for the 2081–2100 time period. The data were derived from the UKESM1-0-LL climate model under the SSP585 scenario, which represents a high greenhouse gas emissions pathway.

    • Future bioclimatic variables were acquired at the same spatial resolution as the historical dataset (10 minutes) to ensure consistency across analyses.

    • The future climate .tif files were imported into R and aligned with the historical variable names to ensure compatibility during modeling and projection steps.

  2. Variable Preparation:

    • Both historical and future bioclimatic datasets were organized as SpatRaster objects using the rast function from the terra package, enabling their efficient use as environmental predictors in the modeling process.

    • The Coordinate Reference System (CRS) of the future climate variables was harmonized with that of the historical datasets to maintain spatial consistency across analyses.

  3. Selection of Bioclimatic Variables:

    • To optimize model performance and limit dimensionality, only four key bioclimatic variables were selected based on their ecological relevance for Testudo graeca ibera. These variables represent critical environmental factors influencing habitat suitability and persistence of terrestrial tortoise species.

By incorporating both historical and future bioclimatic variables, the model was able to characterize the current ecological niche of Testudo graeca ibera and assess potential shifts in habitat suitability under projected climate conditions. This approach is essential for evaluating the potential effects of climate change on the species’ spatial distribution.

## Modeling

Species distribution models were fitted using the BIOMOD_Modeling function from the biomod2 package, applying Generalized Linear Models (GLM). Model calibration employed a random cross-validation approach with three replicates, allocating 80% of the data to model training. Model accuracy was assessed using the True Skill Statistic (TSS) and the Area Under the Receiver Operating Characteristic Curve (ROC). A total of 10,000 pseudo-absence points were generated using a random sampling strategy.

1.2 Projections

1.2.1 Current Projection

Current habitat suitability was projected using historical bioclimatic variables (BIO1, BIO5, BIO6, BIO12, BIO15, and BIO17) through the BIOMOD_Projection function.

1.2.2 Future Projection

Future climate projections for the 2081–2100 period under the SSP585 scenario were prepared for analysis. Following the harmonization of layer names with the historical dataset, habitat suitability projections were generated using the same GLM model.

1.3 Model Evaluation

Model performance was evaluated using the get_evaluations() function, with mean TSS and ROC values calculated across all cross-validation replicates.

2 Results

2.1 Present Projection

The projected current habitat suitability for Testudo graeca ibera reveals a spatial pattern largely concentrated in the eastern Mediterranean basin and adjacent regions. The GLM-based model identifies high suitability across Anatolia, the southern Balkans, and parts of the eastern Mediterranean coastline, with additional suitable areas extending into the Caucasus region and surrounding lowland habitats. Moderate suitability is also observed in parts of the Middle East and western Asia, reflecting the species’ tolerance to semi-arid and Mediterranean climatic conditions. Overall, the model highlights a relatively constrained but ecologically coherent distribution, closely aligned with the known native range of Testudo graeca ibera. Areas characterized by extreme climatic conditions, such as northern Europe and highly arid desert regions, are consistently predicted as unsuitable. These results suggest that temperature extremes and precipitation patterns play a key role in limiting the current distribution of the species, with suitable habitats primarily restricted to regions exhibiting mild winters, warm summers, and moderate seasonal precipitation variability.

2.2 Future Projection

Under the SSP585 scenario for the 2081–2100 period, future habitat suitability projections for Testudo graeca ibera indicate a pronounced shift in suitable areas toward higher latitudes and more northern regions. Compared to current conditions, suitability across the eastern Mediterranean basin and parts of Anatolia is markedly reduced, with many coastal and southern areas becoming increasingly unsuitable. In contrast, regions in the northern Balkans, parts of Eastern Europe, and higher-altitude areas show relatively higher suitability under future climate conditions. Overall suitability across North Africa and the Middle East declines substantially, reflecting the species’ sensitivity to increased temperatures and altered precipitation regimes. The projected distribution suggests a contraction of suitable habitats from traditionally occupied southern and lowland regions, accompanied by a relative expansion toward cooler and more temperate zones. These patterns indicate that future climate change may substantially constrain the geographic range of Testudo graeca ibera, potentially forcing the species toward climatically marginal but cooler habitats.

The following maps illustrate the spatial distribution of the bioclimatic variables used in the species distribution model for Testudo graeca ibera. Each map represents the geographic variation of a specific climatic factor, such as temperature and precipitation, across the study area. Color gradients indicate differences in variable intensity, providing an environmental context for understanding the species’ predicted habitat suitability.

The map below represents the predicted current distribution of Testudo graeca ibera*, derived from species distribution modelling using present-day bioclimatic variables.

2.3

3 Discussion

3.1 Future Distribution of Testudo graeca ibera

Based on future climate projections for the 2081–2100 period under the SSP585 scenario, Testudo graeca ibera is expected to undergo a pronounced contraction and spatial shift in habitat suitability. Compared to current conditions, suitability markedly decreases across much of the eastern Mediterranean region, including large parts of Anatolia and surrounding lowland areas. While some relatively suitable habitats persist in higher-latitude regions of Eastern Europe and the northern Balkans, overall suitability within the species’ traditional Mediterranean range declines substantially. Regions characterized by increasingly warm and dry conditions, particularly North Africa, the Middle East, and southern coastal zones, show a strong reduction in habitat suitability. In contrast to expectations of broad range expansion, the projections indicate only limited potential persistence in cooler and more continental environments rather than widespread northward expansion. Overall, the results suggest that future climate change may restrict Testudo graeca ibera to a narrower range of climatically suitable areas, leading to a loss of habitat in regions that currently support stable populations, especially within Mediterranean and semi-arid ecosystems.

3.2 Comparison to Literature

In light of the updated future projections, the bioclimatic preferences of Testudo graeca ibera—which favor temperate and Mediterranean conditions—are only partially maintained, while notable departures from earlier expectations become evident. Although previous assumptions often implied that warming climates might allow persistence or gradual shifts toward cooler regions, the present model suggests that altered temperature extremes and changing precipitation patterns, including increased drought intensity and seasonal variability, may substantially limit the species’ future distribution. Rather than retaining suitability across southern and coastal Mediterranean regions, the projections indicate a pronounced decline in habitat suitability throughout much of Anatolia, the Middle East, and North Africa. Limited areas of relative suitability persist in more northern and continental regions, such as parts of Eastern Europe and the northern Balkans, where climatic conditions remain comparatively moderate. This pattern contrasts with earlier expectations of stable suitability within traditional Mediterranean habitats and highlights a potential contraction toward climatically marginal but cooler environments. These discrepancies may be attributed to differences among climate models, the choice of emission scenarios, and the exclusion of additional ecological factors such as habitat fragmentation, land use change, and biotic interactions, which can strongly influence real-world species persistence. Consequently, while Testudo graeca ibera continues to exhibit an affinity for temperate climatic conditions, its projected future distribution appears more restricted than previously anticipated, reflecting the complex and region-specific impacts of climate change on Mediterranean and semi-arid ecosystems.

3.3 Model Performances

##                                  full.name      PA    run algo metric.eval
## 1        Testudo.graeca.ibera_PA1_RUN1_GLM     PA1   RUN1  GLM         TSS
## 2        Testudo.graeca.ibera_PA1_RUN1_GLM     PA1   RUN1  GLM      AUCroc
## 3        Testudo.graeca.ibera_PA1_RUN2_GLM     PA1   RUN2  GLM         TSS
## 4        Testudo.graeca.ibera_PA1_RUN2_GLM     PA1   RUN2  GLM      AUCroc
## 5        Testudo.graeca.ibera_PA1_RUN3_GLM     PA1   RUN3  GLM         TSS
## 6        Testudo.graeca.ibera_PA1_RUN3_GLM     PA1   RUN3  GLM      AUCroc
## 7      Testudo.graeca.ibera_PA1_allRun_GLM     PA1 allRun  GLM         TSS
## 8      Testudo.graeca.ibera_PA1_allRun_GLM     PA1 allRun  GLM      AUCroc
## 9  Testudo.graeca.ibera_allData_allRun_GLM allData allRun  GLM         TSS
## 10 Testudo.graeca.ibera_allData_allRun_GLM allData allRun  GLM      AUCroc
##    cutoff sensitivity specificity calibration validation evaluation
## 1   654.0      98.993      95.588       0.946      0.938         NA
## 2   781.0      98.322      96.025       0.985      0.988         NA
## 3   693.0      98.322      95.975       0.943      0.923         NA
## 4   780.5      97.651      96.400       0.986      0.982         NA
## 5   782.0      97.651      96.150       0.941      0.921         NA
## 6   780.5      97.987      96.138       0.986      0.984         NA
## 7   625.0      98.660      95.670       0.943         NA         NA
## 8   709.0      98.391      95.900       0.985         NA         NA
## 9    14.0      97.319      95.330       0.925         NA         NA
## 10   18.5      97.051      95.840       0.984         NA         NA

4 References

-   [1] Fritz, U., & Havaš, P. (2007). Checklist of chelonians of the world. Vertebrate Zoology, 57(2), 149–368.
-   [2] Vamberger, M., Stuckas, H., Sacco, F., D’Angelo, S., Arculeo, M., Cheylan, M., & Fritz, U. (2017). Differences in gene flow and genetic structure in the spur-thighed tortoise (Testudo graeca) across the Mediterranean region. Zoology in the Middle East, 63(2), 113–128.
-   [3] IUCN. (2023). Testudo graeca. The IUCN Red List of Threatened Species. International Union for Conservation of Nature.
-   [4] Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., & Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology, 25(15), 1965–1978.
-   [5] Fick, S.E., & Hijmans, R.J. (2017). WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12), 4302–4315.
-   [6] Thuiller, W., Georges, D., Engler, R., & Breiner, F. (2023). biomod2: Ensemble platform for species distribution modeling. R package version 4.2-5.
-   [7] Araújo, M.B., & Peterson, A.T. (2012). Uses and misuses of bioclimatic envelope modeling. Ecology, 93(7), 1527–1539.
-   [8] Urban, M.C. (2015). Accelerating extinction risk from climate change. Science, 348(6234), 571–573.