Connecting protected areas in the North Mesopotamian steppes: can this ensure the survival of the Arabian Sand Gazelle (Gazella marica)?

The Arabian Sand Gazelle (Gazella marica Thomas, 1897) is endemic to the Middle East and its populations in South-eastern Anatolia in Turkey are in danger of extinction. Here, we focus on its protected habitats in the Mesopotamian steppes and examine how to create connectivity between habitats under changing climate conditions. We use ecological niche modelling for four time periods, as simulated in five General Circulation Models (GCMs), to assess the distributional shifts of G. marica over the decades and predict its future ecological connection effectivity based on a least-cost model. We suggest establishing three potential corridors between core protected areas. The findings indicate that there will be a southward and westward shift in the distribution of the species due to future climate change. Overall, the results of the GCMs suggest that climate change will trigger habitat loss for G. marica. Unless conservation measures take the forecasted changes in habitat into consideration, the extinction of this species will be unavoidable.


Introduction
The Arabian Sand Gazelle (Gazella marica Thomas, 1897) is the largest extant terrestrial mammal species native to the Middle East, with isolated populations in Saudi Arabia, the United Arab Emirates, Oman, Jordan, and South-eastern Turkey. The species is currently categorised as Vulnerable (VU) by the IUCN Red List (IUCN, 2017). The total number of Sand Gazelles has been estimated at 1,750-2,150 mature individuals with the largest numbers being found in the United Arab Emirates (approx. 800-1,000 individuals) and Saudi Arabia (1,000-1,100) (IUCN, 2017). The species, which has been subject to heavy poaching throughout its distribution area in recent years, has in fact been hunted ever since the Neolithic era (Bar-Oz et al., 2011;Lang et al., 2013): until the adoption of livestock husbandry during the Neolithic (8500 BCE), human groups in the Middle East hunted gazelles for meat (Horwitz et al., 1999;Lang et al., 2013).
In Turkey, the distribution area of the Arabian Sand Gazelle extended all the way from Cizre to Gaziantep in the early 1900s, but it is now limited to Şanlıurfa, where the species is protected at the Kızılkuyu Wildlife Reserve (Durmuş, 2010). Even with are located on the Turkish side of the border with Syria and extend to the Taurus Mountains. The province of Şanlıurfa includes three core areas, of which two -the Kızılkuyu Wildlife Reserve (WR) and the Tektek Mountains National Park (NP) -are protected areas and one -the Karacadağ Reserve Area (RA) -is an "important biodiversity reserve area". The average altitude in the region, which is also known as the Şanlıurfa plateau, is 518 m. The terrain is rugged and generally composed of limestone and volcanic basalt. A semi-arid, cool Mediterranean climate prevails, with an average annual precipitation of 463.3 mm (Emberger, 1955) and an average temperature of 18.5°C. The highest average temperature is in July (32.0 °C) and the lowest in January (5.6 °C). The province of Şanlıurfa is home to 1,668 plant taxa belonging to 88 families with an endemism rate of 6% (Kaya & Karataş, 2019).
The steppes in the area, which form an important part of the Mesopotamian steppe ecosystem, are under threat due to a high degree of habitat loss, improper land practices, and climate change (Ceylan, 2019). Due to the diminishing ecological connectivity, we focused on the three core areas among which connectivity might be established. Of these, the Karacadağ Reserve Area (RA) (60,000 ha) is located to the north-east, the Kızılkuyu Wildlife Reserve (WR) (15,337 ha) to the south-west, and the Tektek Mountains National Park (NP) (19,335 ha) to the east of the city centre of Şanlıurfa (FAO-MAF, 2022). The Kızılkuyu WR is the largest habitat of G. marica. A large portion of the Tektek Mountains NP is dominated by turpentine trees (Pistacia terebinthus) (FAO-TOB, 2019). Occurrence data: A total of 25 occurrence records were obtained for G. marica from the Global Biodiversity Information Facility (GBIF: 12 records, www.gbif.org/), from the available literature (10 records), and during the course of this study (3 records). All records were georeferenced using a WGS84 coordinate system and checked for accuracy with ArcGIS (v10.7, ESRI, California, USA). To minimize sampling bias which might otherwise result in the overestimation of the predicted distribution (Merow et al., 2013), and to reduce spatial autocorrelation (Boria et al., 2014;Fourcade et al., 2014), we drew a 10-km buffer area around each occurrence record using the spThin ver. 0.2.0 (Aiello-Lammens et al., 2015) package, reducing the total number of records from 25 localities to 21.
We hypothesized the study area accessible to the Sand Gazelle during present and future time periods (Soberón & Peterson, 2005). In order to predict the future projections of the species in the Middle East, we determined the study area for the species using a five-degree buffer around known occurrence data points and cropped it as a polygon. We masked the climatic variables for the south of South-eastern Anatolia and the whole of the Arabian Peninsula as a polygon between 9.90°N and 39.11°N, and 34.33°E to 60.74 o E representing the study area (Soberón & Peterson 2005).
To reduce the negative effect that might result from multicollinearity and high correlation (r>0.75 or < -0.75) among the bioclimatic variables (Dormann et al., 2013;Heikkinen et al., 2006), we removed the variance inflation factor lower than 10 and employed a correlation threshold of 0.75, as recommended by Guisan et al. (2017) using the 'usdm' package (Naimi et al., 2014). Based on the ecological requirements of the species, six bioclimatic variables were selected as predictors -namely: BIO2 = Mean Diurnal Range (Mean of monthly (max temp -min temp)); BIO3 = Isothermality; BIO5 = Maximum Temperature of Warmest Month; BIO12 = Annual Precipitation; BIO14 = Precipitation of Driest Month, and BIO15 = Precipitation Seasonality (Coefficient of Variation).
The future habitat suitability of the G. marica was estimated by projecting the fitted models to five different global circulation models - ( (Shiogama et al., 2019) -so as to account for a reasonable degree of uncertainty in the climate model projections (Sanderson et al., 2015). Future datasets were obtained for the time periods 2021-2040, 2041-2060, 2061-2080, and 2081-2100 from the 6th Climate Model Intercomparison Project (CMIP6, www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6) with three shared socioeconomic pathways (SSPs) (optimistic -SSP245, middle of the road -SSP370 and worst -SSP585).

Ecological niche modelling:
We built ecological niche models to forecast the present-time and future habitat suitability of the G. marica using an ensemble method in the sdm package (Naimi & Araújo, 2016) in the R v3.6.3 environment. We implemented five algorithms using different approaches: BIOCLIM (Busby, 1991), the generalized linear model (GLM;McCullagh & Nelder, 1989), boosted regression trees (BRT) (Friedman, 2001), random forests (RF) (Breiman, 2001), and maximum entropy (MaxEnt) (Phillips et al., 2006) with 10,000 randomly selected pseudoabsences. We ran all algorithms with standard parameterization from the sdm package (Naimi & Araújo, 2016). Due to the small sample size, model evaluation was performed with the subsample and bootstrapping resampling methods (Hastie et al., 2009) and split into 70-30% subsets for model calibration and testing. We ran each model 100 times with calculated true skill statistics (TSS) and the area under the receiver operating characteristic (AUC). In total, we calculated 1,000 individual models (5 algorithms x 2 resampling methods x 100 replications) to determine the present-time habitat suitability of the species. We selected the best models with TSS>0.7 and AUC>0.9 to build ensemble models for each scenario. We made use of the mean of predicted presence-absence values method (wherein the predicted probability of occurrences is first converted to presence-absence using a threshold, then averaged) for ensembling the best models. Subsequently, the selected models were projected into the current and future scenarios. We therefore implemented 60 projections (5 GCMs × 3 SSPs × 4 time-periods) and produced the ensemble rasters for the GCMs scenarios and time periods. The outputs represent habitat suitability on a scale from 0 (unsuitable) to 1 (suitable). We transformed the ensemble suitability models into binary maps of suitable environmental conditions and applied them by maximizing the sum of sensitivity and specificity (maxSSS) as recommended by Liu et al. (2005). The results were visualized with the RasterVis package (Lamigueiro et al., 2022).

Corridor model:
In the design of ecological corridors, landscape resistance maps were prepared based on assumptions of the extent to which the landscape matrix affects or limits the mobility of Figure 2. Average prediction of climate habitat suitability maps for the Arabian Sand Gazelle (Gazella marica) projected to the present day. Red dots show occurrence records. The probability of occurrence ranges from 0 (dark purple, low probability) to 1 (yellow, highest probability). the target species. Forms of land use, roads and railways, and hydrological structure and irrigation channels, were taken as variables in the creation of resistance surfaces (Appendix 1, Zeller et al., 2012). Previous studies on different species using similar methodology were used to determine the resistance values of the variables (e.g. Loro et al., 2016;Özcan & Erzin, 2020). To increase resistance sensitivity, the resistance values were set at between 0 and 1,000. The hydrological structure was used as an indicator in determining the resistance values. The topographic structures formed on the terrain by water can generally help wild animals move more easily. Places with hydrological networks running through them were accorded lower resistance values than places with various forms of land use. Irrigation canals create significant obstacles to the gazelle in the region. The intensive agriculture practiced in the Harran Plain, which is located between two of the core areas (Kızılkuyu WR and Tektek Mountains NP), is based entirely on irrigation canals that constitute a major obstacle to wildlife mobility. Forests have the lowest landscape resistance values, while settlements, mines and irrigation channels have the highest values. Forest plans and LCC (CORINE) maps were used to determine the land use (Copernicus Land Monitoring Service, 2022). The dataset of the General Directorate of Highways was used tom determine the average daily traffic density of the transport network and the State Hydraulic Works (DSI) dataset was used to obtain information about the irrigation canals. The scale used was 1:100,000 and the raster resolution 100 m.
In determining ecological connections and the locations of potential corridors, the circuit model was used following the graph theory approach (West, 2001). Here, graphs are networks consisting of clusters of nodes (connections that represent habitat patches, populations, or cells in a raster landscape) connected by edges. The simplest connectivity measure derived from circuit theory is the resistance distance (Klein & Randić, 1993). A convenient feature of the resistance distance is that it includes multiple paths connecting nodes, since the resistance distances meas-ured between pairs of nodes decreases over time as more connections are added. In other words, this indicator incorporates both the minimum travel distance/cost and the availability of alternative routes . The GIS-based Least-Cost Path (LCP) tool was used to estimate ecological connections through the graph method . The LCP analysis based on graph theory uses a raster-based algorithm that measures the minimum travel cost between a resource and a target cell (Bunn et al., 2000). We calculated the corridors between each resource point pair using the LCP in the ArcGIS (v10.7, ESRI, California, USA) Spatial Analyst Tool with the least cumulative cost and direction raster. We used two different layers: the spatial layer of core areas (Protected Areas) and the resistance layer within the landscape against the mobility of the target species. The width of the corridor was only one pixel. The least-cost path should be interpreted as a potential pathway that minimizes the cost of mobility, rather than the functional expression of the dispersal process (Gurrutxaga et al., 2010;Theobald, 2006).

Results
Overall, our ecological niche models (ENMs) have an average AUC of 0.860 (SD = 0.133) and an average TSS of 0.729 (SD = 0.258). The best models demonstrated good performance, with an average AUC of 0.995 (SD = 0.018) and an average TSS of 0.988 (SD = 0.040). Our ENM for present-day conditions indicated that the habitats suitable for the Arabian Sand Gazelle are patchily spread across the south of South-eastern Anatolia and the Arabian Peninsula ( Figure 2). Our results show that the habitat suitability of the species is explained by the maximum temperature of the warmest month (BIO5, 23%), mean diurnal range (BIO2, 20%), precipitation seasonality (BIO15, 18%), precipitation of the driest month (BIO14, 16%), annual precipitation (BIO12, 15%), and isothermality (BIO3, 8%). Except for isothermality, the other environmental variables shape the distribution of these areas almost equally.
The future projections estimate that the potential distribution of the species will diminish significantly in the 2021-2100 period (Figure 3). Under all the scenarios (optimistic -SSP245, middle of the road -SSP370, and worst -SSP585), it is predicted that the potential habitats of the species, especially those in the centre and south of the Arabian Peninsula, will diminish. In the moderate and pessimistic scenarios, it is estimated that these areas will almost disappear and that there will also be serious losses in South-eastern Anatolia, Syria, and Jordan -i.e., in the north of the species range. The potential habitats of the species in Yemen and Oman are predicted to shrink in 2061-2080 and 2081-2100. In conclusion, the future projections point to varying but serious losses in the potential habitats of the species and indicate that more stringent protection activities will be necessary in the future.
Two different corridors were identified for the Sand Gazelle between the Kızılkuyu WR and the Tektek Mountains NP (Figure 4). The most suitable corridor between the Kızılkuyu WR and the Tektek Mountains NP passes through the border area between Turkey and Syria. Since the border zone and its surroundings are closed to human use, their landscape resistance is low. This corridor is therefore treated as an alternative connection. This route lies towards the north of the Harran Plain. Likewise, this corridor passes through the west of the province of Şanlıurfa, which is located in the north-west of the Harran plain. According to CORINE data, approximately 80% of the corridor passes through areas of pasture and sparse vegetation, which may be described as natural. The corridor identified to link the Kızılkuyu WR and Karacadağ RA traverses the north of Şanlıurfa and the south of Siverek. A large portion of this corridor also passes through natural areas according to CORINE data. The corridor between the Tektek Mountains NP and the Karacadağ RA connects these two locations through the north of Ceylanpınar and the west of Viranşehir.
As a result of the corridor model, three corridors (not including the alternative corridor between the Kızılkuyu WR and the Tektek Mountain NP) were identified to link the core areas. The major barriers were also identified. The shortest corridor in the study area, 73 km in length, was between the Kızılkuyu WR and the Tektek Mountains NP (Corridor 1 in Figure 4), while the longest corridor, 146 km in length, was between the Kızılkuyu WR and the Karacadağ RA (Corridor 2 in Figure 4). The average length of the corridor between the Tektek Mountains NP and the Karacadağ RA was 117 km.  (2021-2040, 2041-2060, 2061-2080, and 2081-2100) and three shared socioeconomic pathways (optimistic -SSP245, middle of the road -SSP370 and worst -SSP585). The probability of occurrence ranges from 0 (dark purple, low probability) to 1 (yellow, highest probability). The lengths of the identified potential corridors by CORINE land cover classes (LCC) are presented in detail in Appendix 2. Most of the route of Corridor 1 passes through agricultural areas (122 km), followed by forest and semi-natural areas (22 km). Most of these agricultural areas, accounting for 75 km of the route, are pastures, which are natural habitats. Corridor 2 passes mainly through forest and semi-natural areas, which account for 75 km of the route. Agricultural areas follow with 42 km. The CORINE LCC that accounts for the greatest stretch of Corridor 3 is forest and semi-natural areas, with a similar total length to Corridor 2 (48 km). What is common to all three corridors is that the lengths of corridor passing through areas classed as artificial surfaces are rather short compared to the lengths of corridor passing through other classes of land. Meanwhile, none of the three corridors pass through areas classed as artificial, nonagricultural vegetated areas or as wetlands and waterbodies.

Discussion
The wild populations of gazelle species have decreased dramatically in the last 50 years, and almost all gazelle species are considered endangered (IUCN Red List, 2021). Kemal (1955) stated in his review that hundreds of Arabian Sand Gazelles used to migrate between Ceylanpınar and Mount Abdülaziz, Syria. However, it seems that this migration route disappeared completely after 1960. The processes that accelerated the decline of the species in South-eastern Anatolia were the rapid transformation of natural habitats into agricultural areas, the conflict with sheep breeders in the remaining pastures, and increasing hunting activity. The dramatic decline in Arabia, where the largest populations of the species are located, occurred in the 1970s, at almost the same time as the northernmost distribution (Nader, 1989). Al Hikmani et al. (2015) stated that decades of heavy poaching had reduced the presence of the species not just in Oman but throughout its entire range in the Arabian Peninsula. Subsequently, the species was taken under protection and efforts were made to reintroduce it through breeding and breeding-and-reintroduction programmes. Today, most populations of the species are to be found in protected areas. Between 2,650 and 3,050 individuals live in four populations in the protected areas of Al-Khunfah, Harrat al Harrah, Mazahat as-Sayd and Uruq Bani Ma'arid in Saudi Arabia (Cunningham & Wacher, 2009). Gazelles are also under protection in Jordan, and protected areas such as the Yarmouk Forest Reserve have been used for conservation measures (Eid & Mallon, 2021). In Turkey, which has the second largest population, 350 individuals were living in the Kızılkuyu WR as of 2019 and three individuals were identified in the Tektek NP in the same year. Some members of the species are also living under protection in reproduction farms and new reintroduction sites (FAO-TOB, 2019).
Overall, the numbers of Sand Gazelles living in the wild has been reduced drastically due to heavy poaching and the range of the species has declined noticeably due to habitat degradation/loss. The total population in the Middle East somewhat exceeds 3,000 -or approximately 2,100 mature individuals. When captive or "managed" populations are taken into account, the total number of individuals is probably more than 100,000 (IUCN, 2017).
With climate change, new problems await the gazelles in Turkey, particularly because the country forms the northernmost distribution of the species. Species can survive this rapid process of change either by changing their distribution areas or by evolving to adapt to new local climatic conditions (Berg et al., 2010). For this reason, predicting and managing the responses of species to climate change requires a full consideration of the distribution characteristics of the species in question and of how these characteristics may alter with climate change .
To some extent, efficient management of the protected areas will buffer the consequences of global warming. Currently, almost all the populations of G. marica live in controlled protected areas. The management approaches adopted in the protected areas focus on maintaining and improving the connectivity between these areas to reduce and stabilise the pressure of habitat loss and fragmentation on biodiversity and to increase the resistance of reserve networks to potential threats associated with climate change (Bhola et al., 2021). In this study, we used the least-cost method to determine three corridors and one alternative corridor for connecting the three core areas within an ecological network. The estimated lengths of the corridors identified in our analysis are within a suitable distribution distance range for gazelles. Legge & Rowley-Conwy (1987) used desert kites to show that the migration routes of gazelles on the North-South axis could reach up to 500 km.
The most important problems affecting the connections between the protected areas in the study were linear barriers -namely, roads and irrigation channels. Measures to be taken for conservation planning are inevitably costly (Hetherington et al., 2008). Where highways intersect with potential ecological corridors, a sufficient number of suitable crossing points need to be established to permit the dispersal of medium-large mammals (Özcan & Erzin, 2020). At least two roads with a very high traffic density must be crossed along the two corridors that connect the Karacadağ RA with the other areas. Wildlife overpasses and underpasses may provide connections, improve connectivity and reduce road accidents (Pierik et al., 2016). Although the linear distance between the Tektek Mountains NP and the Kızılkuyu WR is only 40 km, the Harran Plain, which has the largest irrigation system in Turkey, lies directly between them and extends the route considerably. Furthermore, there are two major main irrigation canals on these routes. These irrigation canals must also be bridged using wildlife overpasses similar to those to be provided for the highways. There are already many overpasses with widths of 150-750 m on the irrigation canals. For example, the Kızılkuyu WR is located across the overpasses in the south of the village of Çatallı, where the corridor between the WR and the Tektek Mountains NP intersects with the irrigation channels.
The route of the alternative corridor between the Kızılkuyu WR and the Tektek Mountains NP identified in this study passes through a buffer zone which runs along the border between Turkey and Syria and which is as much as one kilometre wide in places. Since these areas are closed to all human use, the land use has not changed and they have the properties of natural pastureland. However, the border zone is surrounded by agricultural fields, which makes it impossible for gazelles to use this corridor. Nevertheless, we find it important in that it offers an alternative connection between the two core areas: while we have focussed on the gazelles, corridors are also important for other wild animals.
Reserve systems expanded by ecological connections will play a critical role in enhanced protection efforts and efforts to reduce the impacts of climate change on biodiversity and ecosystem services. Even now, gazelles are on the verge of extinction unless effective conservation measures are implemented. Prioritising protection efforts for this species will also provide secondary protection to other animal and plant species of the Mesopotamian steppes. The management approaches in protected areas focus on maintaining and improving connectivity between protected areas to reduce the pressure of habitat loss and fragmentation on biodiversity and to increase the resistance of reserve networks to potential threats associated with climate change (Bhola et al., 2021). This study predicts that reduced connectivity might have adverse impacts on habitat functions, as there will be significant losses in the habitats of G. marica. Furthermore, previous studies have shown that climate change has adverse impacts on many large mammals including G. marica (Hetem et al., 2012;Isaac, 2009). This study indicates that there will be a southward and westward shift in the distribution of the species due to climate change in the future. The Sand Gazelles in our study area are found in the northernmost distribution of the species. In its current potential distribution, the species is most widely dispersed in the Syrian desert. According to future climate scenarios, even in the best emission scenario, the potential habitats in the Syrian desert are moving south and west. The species will increasingly be squeezed up against the Al Nusayriya and Lebanon mountain ranges that border the desert to the west. Future climate projections for the Syrian deserts (Mourad & Alshihabi, 2016;Homsi et al., 2020) are consistent with our model. Overall, the results of the GCMs suggest that climate change will trigger habitat loss for G. marica. If conservation measures do not take both the potential corridors and the habitat suitability under future climate scenarios into account, the species will inevitably become extinct.

Supplementary Material
Supplementary Material is provided in a Supplementary Annex, which is available via the "Supplementary" tab on the online page of the article.