The Elusive Turkestan Lynx at the Northwestern Edge of Geographic Range: Current Suitable Habitats and Distribution Forecast in the Climate Change

Abstract

The Turkestan lynx (Lynx lynx isabellina Blyth, 1847) is a rare and understudied subspecies of the Eurasian lynx occupying the mountains of Central and South Asia. This elusive felid’s northwestern range includes the Tien Shan and Zhetisu Alatau mountains in the border region of Kazakhstan, China, Kyrgyzstan, and Uzbekistan. As the first step to conserve this vulnerable carnivore, we have conducted the first full-scale research from 2013 until 2022 on its distribution in this region. Using 132 environmental predictors of 359 lynx sightings, we have created species habitat distribution models across the lynx’s northwestern range using machine learning approaches (Maximum Entropy—MaxEnt). Additionally, we created species distribution forecasts based on seven bio-climatic environmental predictors with each three different future global climate model scenarios. To validate these forecasts, we have calculated the changes in the lynx distribution range for the year 2100, making the first species distribution forecast for the Turkestan lynx in the area. Additionally, it provides insight into the possible effects of global climate change on this lynx population. Based on these distribution models, the lynx population in the Northern and Western Tien Shan and Zhetisu Alatau plays a significant role in maintaining the stability of the whole subspecies in its northwestern and global range, while the distribution forecast shows that most lynx distribution ranges will reduce in all future climate scenarios, and we might face the Turkestan lynx’s significant distribution range decline under the ongoing and advancing climate change conditions. For a future (year 2100) warming scenario of 3 deg. C (GCM IPSL), we observe a decrease of 35% in Kazakhstan, 40% in Kyrgyzstan, and 30% in China as the three countries with the highest current predicted distribution range. For a milder temperature increase of 1.5–2 deg. C. (GCM MRI), we observe an increase of 17% Kazakhstan, decrease of 10% in Kyrgyzstan, and 57% in China. For a cooling scenario of approx. 1–1.5 deg. C (GCM MIROC), we observe a decrease of 14% Kazakhstan, increase of 11% in Kyrgyzstan, and a decrease of 13% in China. These modeled declines indicate the necessity to create new and expand the existing protected areas and establish ecological corridors between the countries in Central and South Asia.

1. Introduction

The development of industry and agriculture, as well as a significant increase in the world’s population, have led to an increase in anthropogenic pressure on natural ecosystems, fauna and flora [1,2,3]. In particular, human activities have precipitated habitat loss and fragmentation, followed by isolation, loss of genetic diversity and biodiversity, decline in population sizes, and extinction of wild animals [4,5,6,7,8]. Among vertebrates, a number of large carnivorous mammals are threatened with extinction due to their lower population density and reproductive rates. Moreover, some of the main threats to carnivores are related to people’s prevailing negative attitude towards these animals, which may result in their direct or indirect extermination [9].

The Eurasian lynx (Lynx lynx L., 1758) is a felid species especially affected by human activities [10]. For instance, as a species that was once widespread throughout Europe, it became extinct in the region in the 19th century due to habitat loss and persecution [11,12]. The lynx is still considered a rare or endangered species, even after its re-introduction in a number of countries [13,14]. The case of lynx disappearance in Europe is a good example that indicates the need for measures to study and preserve this species and its subspecies in other habitats, taking into account the evolving industry and intensive developing of natural areas.

The Turkestan lynx, also known as Himalayan lynx, Central Asian lynx (Lynx lynx isabellina, also referred to as Lynx lynx isabellinus Blyth, 1847) is an elusive subspecies of the Eurasian lynx that occupies the mountains of Central and South Asia. Being rare and understudied throughout all 11 countries of its habitat, the lynx is listed in Appendix II of the International CITES Convention. In Central Asia, this carnivore is listed in the Red Data Book of Kazakhstan as “a rare subspecies, the range and number of which are decreasing” (III category) [15]; Kyrgyzstan as a “subspecies close to a vulnerable position” (NT, VI category), [16]; Uzbekistan as a “vulnerable declining subspecies” (II category—VU: D) [17]; Turkmenistan as “an extinct or endangered species or subspecies” (Category I) [18]; Tajikistan as an “endangered subspecies” (EN) [19]; China as an “endangered subspecies” (EN A1cd) [20]; and in Afghanistan, the lynx is considered “vulnerable” [21,22]. In South Asia, the lynx is common in the northern part of Pakistan, and it is considered to be of “Least Concern” (LC) in the region [23]. However, the current status and distribution in Pakistan is unknown [24], which confirms the need to update the Red List in Pakistan [23,25]. In India, the lynx is on the verge of extinction [26,27], and is listed in Plan I—protected by the National Indian Wildlife Protection Act since 1972. The lynx is listed in the Red Data Book of Nepal as a “vulnerable subspecies” (Category Vulnerable B1a; D2; protected under the National NPWC Act 2029) [28]. In Bhutan, there is no information available to determine lynx status (Data Deficiency).

The study is complicated by the inaccessibility of lynx habitats, as well as its secretive lifestyle. As it is exceedingly elusive, this form of the lynx is poorly studied, especially in the northwest of its global range, where it inhabits the mountains of Northern and Western Tien Shan and Dzhungarian (Zhetisu) Alatau [29]. The mentioned mountains are located in the border region of Kazakhstan, China, Kyrgyzstan, and Uzbekistan.

Most of the early studies conducted on a local scale in Kazakhstan in the 20th century were mainly related to the general distribution of the species, indicating the lynx presence in a particular mountain range [30,31,32], or contain fragmentary data on lynx sightings obtained while studying other mammalian species [33,34,35,36,37,38,39,40,41]. The situation is similar in neighboring Kyrgyzstan [16,40,42,43,44,45,46,47,48,49,50,51,52,53] and Uzbekistan [54,55,56,57,58,59,60,61,62,63,64,65]. More extensive studies in China by Ablimit et al. [66] indicate that, in the northwestern region, which includes the Xinjiang Autonomous Region, suitable lynx habitats in this region have decreased by 50% since the 1950s. At the same time, data on lynx in China are still insufficient [67], and detailed studies on lynx distribution are needed in the country [68].

In the northwest of its distribution, one of the key lynx habitats in Central Asia, the Turkestan lynx populations are understudied and heavily affected by habitat degradation and fragmentation, prey base depletion, poaching, and conflict with livestock breeders (fur trade and retaliatory killings due to livestock depredation) [69]. With the intensive regional development, the habitats may reduce in size and quality, thus increasing the extinction rate of lynx populations and subpopulations in these mountains [70]. As the first step in preserving this felid, we will predict the subspecies’ distribution, as it is an important component of populations and species conservation planning, and a variety of modeling techniques have been developed for this purpose [71].

Our study goals are twofold: (i) relying on 132 relevant environmental characteristics and our recently obtained more accurate lynx occurrence records, assess the habitat suitability for the lynx populations and create the first subspecies distribution model (SDM) in its northwest range; (ii) based on seven bio-climatic environmental characteristics, modelling a forecast on the lynx habitat distribution—subspecies distribution forecast (SDF) for the year 2100, thus evaluating the effect of climate change on the lynx’s distribution range. For this purpose, we will apply the above-mentioned models, also used in different modifications by other researchers to study the Eurasian lynx distribution in the world [72,73,74,75].

2. Materials and Methods

2.1. Study Area

We collected the data from the Northern Tien Shan, Western Tien Shan, and the Zhetisu Alatau mountains (Figure 1). We conducted our research within Kazakhstan, located in the most northwestern edge of the Turkestan lynx’s geographical range, and analyzed the survey and literature data from China, Kyrgyzstan, and Uzbekistan.

• In the Northern Tien Shan, our study area covered the mountain ranges—Ile Alatau (N 42°58′–43°47′, E 74°53′–78°54′), Kungei Alatau (N 42°47′–43°01′, E 77°18′–78°59′), and Terskei Alatau (N 42°28′–42°46′, E 79°13′–80°07′), located in the border area of Kazakhstan and Kyrgyzstan, and Uzynkara (Ketmen) (N 43°02′–43°25′, E 79°17′–80°40′), located in the border area of Kazakhstan and China. In the Xinjiang part of Northern Tien Shan, we considered the felid occurrences in the Tomur (N 41°54′, E 80°18′), Kalajun–Kuerdening (N 42°59′–43°05′, E 82°22′–82°59′), Bayinbuluke (N 42°47′, E 84°09′) and Bogda (N 43°51′, E 88°09′) areas. We also studied Kyrgyz Alatau (N 42°38′–42°40′, E 73°14′–74°35′), which is located between Kazakhstan and Kyrgyzstan, and is an interjacent ridge between Northern and Western Tien Shan.
• In the Western Tien Shan, we checked the lynx occurrences in the ranges of Karatau (N 42°53′, E 69°58′), located in Kazakhstan; Karzhantau (N 41°50′–42°03′, E 69°50′–69°58′), Talas Alatau (N 42°19′–42°25′, E 70°22′–71°51′) and Ugam (N 41°49′–42°04′, E 70°16′–70°25′), located in the border region of Kazakhstan and Uzbekistan; Pskem (N 41°45′–42°02′, E 70°22′–70°41′) and Chatkal (N 41°10′–41°23′, E 69°49′–70°29′), located in the border between Uzbekstan and Kyrgyzstan.
• In the Zhetisu Alatau, located in the border region of Kazakhstan and China, we surveyed the ridges of Sholak (N 43°56′, E 77°56′), Degeres (N 44°00′, E 78°10′), Matai (N 44°10′, E 78°26′), Altynemel (N 44°13′, E 78°34′), Koyandytau (N 44°21′, E 78°49′), Toksanbay (N 44°43′, E 79°04′), Baskantau (N 44°04′, E 80°26′), Karatau (N 44°57′, E 79°38′), Bedzhintau (N 45°02′, E 79°58′), and Tyshkantau (N 44°37′, E 80°08′).

The Tien Shan ridges have a latitudinal or nearly latitudinal strike [76]. The climate in the mountains is continental; nevertheless, the complexity and vegetation of the relief cause contrasts in temperature and the degree of moisture [77,78]. Thus, the average annual air temperatures in Tien Shan decrease to 6.5 °C below zero [79]; at the same time, the southern slopes of the ridges can be 5–10 °C warmer than the northern ones [80].

Most of the plant species are found in the mid-altitude mountain forest belt. Above 2000–2200 m, deciduous forests give way to spruce on mountain–forest dark-colored soils with a high (up to 10%) humus content [81,82]. The important component of the spruce forests on the northern slope of the Northern Tien Shan ridges, particularly for the lynx distribution, is the endemic species, the Schrenk’s spruce (Picea shrenkiana) [83,84]; in the Western Tien Shan—the juniper shrubs (Juniperus pseudosabina, J. sibirica); and in the Zhetisu Alatau—the Schrenk’s spruce mixed with the Siberian fir (Abies sibirica) [30].

In addition to the Turkestan lynx, many other mammalian species inhabit Tien Shan and Zhetisu Alatau. Among them, there are carnivores (potential lynx competitors), in particular, representatives of the family Canidae—the wolf (Canis lupus) and the fox (Vulpes vulpes); the family Ursidae—the Tien Shan brown bear (Ursus arctos isabellinus); the Mustelidae—the short-tailed weasel (Mustela erminea), the least weasel (Mustela nivalis), stone marten (Martes foina), and Asian badger (Meles leucurus); and Felidae—snow leopard (Panthera uncia), manul (Otocolobus manul), and African wildcat (Felis lybica). The lynx preys on the tolai hare (Lepus tolai) and mountain hare (Lepus timidus), the Siberian roe deer (Capreolus pygargus), the grey marmot (Marmota baibacina), the Siberian red squirrel (Sciurus vulgaris exalbidus), and the Turkestan red pika (Ochotona rutila), as well as the Siberian ibex (Capra sibirica), wild boar (Sus scrofa), red deer (Cervus elaphus), and birds [2,30,31,32,41,44,45,67,85,86].

We collected the majority of occurrence data in specially protected natural areas (PAs): Ile-Alatau State National Natural Park (SNNP), Almaty State Nature Reserve (SNR), Kolsai Kolderi SNNP, Sharyn SNNP, Altyn Emel SNNP, Zhongar Alatau SNNP, Aksu-Zhabagly SNR, Sairam-Ugam SNNP, and Karatau SNR, and in the adjacent territories. The lynx populations in the PAs were comparably more stable, with a lower level of human disturbance and activities noted. General environmental conditions given in average for lynx locations within our study area are presented in Appendix A (

Table A1. Environmental predictors.

Predictor Name	Source	Explanation	Citation
BIO1_2_5min–BIO19_2_5min; tmin1–tmin12; tmax1–tmax12; tavg1–tavg12; srad1–srad12; prec1–prec12; Wcaltitude	Worldclim (https://www.worldclim.org/data/worldclim21.html (accessed on 1 June 2022))	These datasets represent most of the climate data utilized for the SDM.	[105]
FAOCC	FAO Geonetwork (http://www.fao.org/geonetwork/) (accessed on 1 June 2022)	This predictor represents the global climate classes.
LC12asc2; VE4	Geospatial Information Authority of Japan (https://www.gsi.go.jp/kankyochiri/gm_global_e.html) (accessed on 1 June 2022)	These predictors represent the global land cover (LC12asc2), and the global vegetation cover (VE4).
GlobalRiversProxy2	A Simple Global River Bank full Width & Depth Database (http://gaia.geosci.unc.edu/rivers/) (accessed on 1 June 2022)	This predictor represents all global mid-size and large rivers.	[106]
GlobalBigRivers11	Global major rivers (https://www.arcgis.com/home/item.html?id=44e8358cf83a4b43bc863646cd695945) (accessed on 1 June 2022)	This predictor represents all global large rivers.
GlobalCities2	Global Cities (https://hub.arcgis.com/datasets/6996f03a1b364dbab4008d99380370ed_0?geometry=-65.394%2C25.931%2C73.737%2C49.818) (accessed on 1 June 2022)	This predictor represents all global cities.
GlobalLakes2	Global Lakes and Wetlands Database (GLWD) (http://www.fao.org/land-water/land/land-governance/land-resources-planning-toolbox/category/details/es/c/1043160/) (accessed on 1 June 2022)	This predictor represents all global lakes and wetlands.
GlobalSnowCoverMonthJan2021_7; FFJan2020_3; FFFeb2020_3; FFMar2020_3; FFMay2020_3, FFJun2020_3; FFJul2020_3; FFAug2020_3; FFSep2020_3; FFOct2020_3; FFNov2020_3; FFJan2021_3	Global Snow Cover and Forest fires (https://neo.sci.gsfc.nasa.gov/view.php?datasetId=MOD10C1_M_SNOW and https://neo.sci.gsfc.nasa.gov/view.php?datasetId=MOD14A1_M_FIRE&year=2020) (accessed on 1 June 2022)	These predictors mainly represent the Global snow cover in the month of January, and the forest fire information for nearly all months of 2020 with exception of April and December as these months were not available.
WorldSoil2	Global Soil characteristics map (https://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/HWSD_Data.html?sb=4) (accessed on 1 June 2022)	This predictor represents all global soil types and its characteristics.
WorldProtectedAreasMerged4	Global Protected areas. (https://www.protectedplanet.net/en/search-areas?geo_type=region&filters%5Bis_type%5D%5B%5D=terrestrial) (accessed on 1 June 2022)	This predictor represents all global protected areas merged into one shapefile.	[107]
WorldMammaldensity4; WorldRodentDensity3; WorldThreatenedMammalDensity3; GlobalBirdDensity2	Global Mammal density. Proximity maps for the world mammal density, world rodent density, world threatened mammal density (https://biodiversitymapping.org/index.php/mammals/) (accessed on 1 June 2022)	These predictors mainly represent global biodiversity densities. In detail, they contain the world mammal density, world rodent density, world bird density, and the world’s threatened mammal density.	[108,109]
GlobalRoadsProxy2	Global Roads—Socioeconomic data and applications center (SEDAC)—Data center in NASA’s Earth Observatory System Data and Information System (EOSDIS) (https://sedac.ciesin.columbia.edu/data/set/groads-global-roads-open-access-v1/data-download) (accessed on 1 June 2022)	This predictor represents the global proximity to all world’s roads. Minor roads may not be included.
HII1	Human Influence Index (HII). (https://sedac.ciesin.columbia.edu/data/set/wildareas-v2-human-influence-index-geographic/data-download) (accessed on 1 June 2022)	This predictor represents the global Human Influence index.
WorldSlope1	Slope. (https://scholarworks.alaska.edu/handle/11122/7151) (accessed on 1 June 2022)	This predictor represents the global terrestrial and aquatic slope.	[94]
World_MAX_RH_JAN–World_MAX_RH_DEC; World_MIN_RH_JAN–World_MIN_RH_DEC	Global Monthly Relative Humidity. (http://palebludata.com/?q=data) (accessed on 1 June 2022)	This predictor set represents the global maximum and minimum relative humidity for the months January to December of the year 2020.	[110]

).

2.2. Occurrence Data

2.2.1. Records Search and Analysis

We searched extensively for observational records within our study area. The literature data used in this study were collected mainly from zoological and ecological reports of the Institute of Zoology of Kazakhstan and Specially Protected Natural Areas (PAs), zoological articles, books, conference materials, methodological manuals, etc. In order to determine the modern lynx distribution, in the Northern Tien Shan, we frequently used the records in the diaries of the inspectors (rangers) of protected areas. In the Xinjiang part of Tien Shan, we mainly relied on the data given in the papers and reports of Ablimit et al. [42], the Leading Group for Application of WNH of Xinjiang Uygur Autonomous Region [87], and Tancheng Dong et al. [88], and we also examined the locations mentioned in craniological notes of lynx skulls by A. Regel (collected in 1880; Skulls No. 1164, 1284, 1287, 1325; Zoological Institute of Russian Academy of Science, St. Petersburg, Russia). One location presented in the gbif.org database was examined as well (accessed on 16 April 2022). In the Western Tien Shan, we often relied on the occurrence data by Bykova et al. [67] and Shakula et al. [89] and the oral communication of S.V. Baskakova. In the Zhetisu Alatau, one location given by V.A. Zhiryakov and B.R. Baidavletov [32] was considered for the lynx occurrence data. The reliability of all regions of lynx encounters, provided in scientific papers and reports, were reviewed and verified. Oral communications were taken from experienced researchers in the area.

2.2.2. Field Data Collection

During research from 2013 until 2022, we relied on traditional methods of field mammalogical research [90,91], which included visual observations and identification of various lynx and their prey traces (paw and hoof prints, prey remains, feces, scratches on trees and rocks, etc.) with their registration on the GPS. We typically collected the data on field trips during winter, as it was more attainable to observe tracks on snow.

Within the Kazakh part of our study area, we completed horse-walking routes with a total length of 1426 km, with 586 km in the Ile Alatau, 185 km in the Kungei Alatau, 74 km in the Terskei Alatau, 32 km in the Uzynkara, and 54 km in the Kyrgyz Alatau ridges of Northern Tien Shan; 108 km in the Talas Alatau and 47 km in the Ugam ridge of Western Tien Shan; and 340 km in the Zhetisu Alatau.

2.2.3. Camera Trapping

Research with the use of camera traps was mainly carried out in areas where the level of anthropogenic pressure was reduced. The camera traps (models Bushnell, Reconyx, Seelock, ScoutGuard, BolyGuard, Browning) were installed in 143 locations from 2013 to 2022 in the most potential lynx habitats (Figure 2), which we selected taking into account our own and survey data, as well as the presence of its prey. For the calculation of lynx and its prey abundance in the study area, we used the formula below (see the results in the Section 4.5).

The average abundance was calculated using this formula:

$$\mathrm{SIC}=\frac{\mathrm{Total} \mathrm{IC}}{\mathrm{Total} \mathrm{camera} \mathrm{trap} \mathrm{days}}\times 100,$$

where IC means independent captures, SIC means IC per season (IC per 100 trap days), and camera trap days represents the number of applied camera traps multiplied by the number of days they were installed for.

2.2.4. Data Classification Methodology

We considered all field data obtained both by using traditional field methods and camera trapping, as well as data from a survey of PA staff about encounters of tracks and the lynx individuals themselves. Considering the fact that the study area occupied a large scale, we followed the data interpretation criteria according to the SCALP categories for the classification according to their level of reliability (C1; C2; C3) [92,93] (see

Table 1. The lynx observations in the study area (according to the SCALP categories) from 2013 to 2022.

Ridge	Occurrence			TOTAL
	C1 (n)	C2 (n)	C3 (n)
Northern Tien Shan	122	127	36	285
Western Tien Shan	31	18	13	62
Zhetisu Alatau	9	2	1	12
TOTAL	162	147	50	359

and Figure 3).

C1 (Category 1): “Confirmed occurrence”—our observations (traditional field methods, camera trapping) and survey (questionnaire) data: verified and reliable data such as (1) dead lynx, (2) captured lynx, (3) videos and photographs of lynx, or (4) samples (e.g., fur).

C2 (Category 2): “Probable occurrence”—data checked and confirmed by a specialist (inspector (ranger), hunting biologist, protected area employee, our data), such as (1) remains of livestock or (2) a wild animal killed by a lynx, (3) documented traces of a lynx or other signs of vital activity, (4) feces, or (5) documented (recorded) and confirmed lynx calls.

C3 (Category 3): “Unconfirmed data” of category 2, such as remains of domestic or wild animals, footprints, feces, or lynx calls, and all unverifiable information, such as encounters of lynx by local residents without attached evidence.

In order to present a more complete modeled distribution, here we merged all occurrence data, and created the SDM using a total of 359 occurrence points (see Figure 4). For completion, and possible future research, we also created the same SDMs, using the same predictor set, for all the three different C-categories separately, presented in Appendix B.

2.3. Environmental Predictors

2.3.1. Current Environmental Predictors

In this study, we created species distribution models (SDMs) based on 132 environmental predictors. These 132 environmental predictors aim to best describe the surrounding environmental and habitat conditions of the studied subspecies (Turkestan lynx (Lynx lynx isabellina)) (see [94]). These are the building blocks of the SDMs, apart from the collected occurrence data [95]. These 132 predictors have been retrieved by MS from Steiner and Huettmann (in-press) [8]. They describe, e.g., the temperature, precipitation, relative humidity, soil conditions, snowfall, proximity to streets, cities and rivers, wildfires, vegetation cover, climate classes, prey densities, etc., on a global scale. All details and sources of the 132 predictors can be found in Appendix A, Table A1 (Reproduced from Steiner and Huettmann (in-press)) [8]. Each environmental predictor has an accuracy of 2.5 km² with a pixel size of 0.04166666666666666435x − 0.04166666666666666435 decimal degrees (CRS: EPSG:4326—WGS 84). This is believed to be the most complete set of environmental predictors publicly available to this date. For details on how to standardize and prepare environmental predictors for SDMs please see Steiner and Huettmann, in-press (Chapter 3) [8].

2.3.2. Projected (Future) Climatic Predictors

We created species distribution forecasts (SDFs) based on seven environmental predictors. These SDFs aim to predict the distribution of Turkestan lynx for the year 2100 using bio-climatic environmental predictors only. These predictors follow the same technical details (pixel size, accuracy) as the ones used for the SDMs above. The reason for only using 7 predictors for these SDFs compared to the 132 predictors for the SDMs is that significantly less data are available for the year 2100 than 2000–2020. Since 2100 is in the relatively far future, for the models here we used three different scenarios. These three scenarios aim to cover three different climate directions that are likely to be approached by 2100. The three scenarios utilized in this study are Global Climate Models (GCMs) that have been downloaded from www.WorldClim.org (accessed on 1 June 2022) (MIROC6, MRI-ESM2-0, and IPSL-CM6A-LR). Similarly, as described by Steiner and Huettmann (in-press—Chapter 7) [8], the MIROC6 may be considered as the low temperature increase scenario, or even as a certain cooling scenario [96]. MRI-ESM2-0 is considered as a low–medium temperature increase as an approximate global increase of 2 degrees Celsius [97]. One might refer to it as a “business as usual” model. Lastly, IPSL-CM6A-LR is considered as a medium–high increase in temperature by approximately three degrees Celsius [98]. Additionally, we created an SMD for the year 2000 with the same predictor set as the SDFs for 2100 in order to best compare current and predicted future distribution ranges.

2.4. MaxEnt

For this study, we used the Maximum Entropy software MaxEnt (https://biodiversityinformatics.amnh.org/open_source/maxent/) (accessed on 1 June 2022) version 3.4.4 to create the SDMs and SDFs (following Machine Learning approaches) (see [99] and references within). A thorough description of the workflow for MaxEnt SDMs can be found in Chapter 3 (Steiner and Huettmann, in-press) [8]. In coarse, the data must first be compiled, cleaned up, and standardized before they can be used for SDMs in MaxEnt. Once the models have been created with MaxEnt, their visualization has been improved and enriched by MS using ArcGIS Pro (version 2.7) (https://pro.arcgis.com/en/pro-app/2.8/get-started/download-arcgis-pro.htm) (accessed on 1 June 2022), and Open-Source GIS—QGIS (versions 3.10 and 3.22.7) (https://qgis.org/en/site/forusers/download.html) (accessed on 1 June 2022). For these models, we have used the standard settings for Maxent; this includes the number of iterations set to 500, and maximum number of background points set to 10,000. Additionally, for all other settings we have used the default option, which includes the feature selection, where we used “Auto Features”, including “Hinge Features”, “Product Features”, “Quadratic Features”, and “Linear Features”.

2.5. Habitat Distribution Percentage

In order to assess the spatial distribution of the models and to provide an overview of in which countries the predicted distribution occurrences can be found, we summarized the spatial information and presented it in percentage per country. This has been carried out by utilizing the QGIS tool “Zonal Statistics”, which calculates the sum of all raster cells (pixels) within a given polygon (in our case countries). Subsequently, the outcomes of this analysis have been converted into percentage per country for better interpretation.

3. Results

3.1. GPS Mapping Results

3.1.1. Species Distribution Models (SDMs)

SDMs have been created for all occurrence datasets available. They have been created for all C-datasets individually (C1, C2, and C3), and for all available data merged. The individual SDMs for all C-datasets can be found in Appendix B (Figure A6, Figure A7 and Figure A8). For a more holistic overview, here we present the SDMs of all available data merged. Figure 4 illustrates the predicted distribution range for the Tien Shan–Alai region.

The SDM indicates that the highest predicted occurrence of Turkestan lynx can be found near the Parkent district of Tashkent Region (Uzbekistan), in the Uzgen district (Kyrgyzstan), near the Naryn River in the Naryn region (Kyrgyzstan), around the Issyk Kul (Kyrgyzstan), in Kemin (Kyrgyzstan), and throughout the Almaty region (Kazakhstan) except the northwest areas near Balkhash Lake. For details, see the Discussion Section 4.3. The top predictors contributing most to the creation of Figure 4 are World Soil Characteristics (with 41.9%), Slope (15.8%), Minimum Relative Humidity of March 2020 (7.1%), Maximum Relative Humidity of October 2020 (5.4%), Maximum Relative Humidity of September 2020 (4.1%), Altitude (3.6%), Solar Radiation of June 2021 (3.1%), and Proximity to World Protected Areas (3%). The reported Area Under the Curve (AUC) for this model was 0.991.

3.1.2. Species Distribution Forecasts (SDFs)

Similar to the SDM, SDFs have been created for all C-datasets, which can be found in Appendix B (Figure A9, Figure A10, Figure A11, Figure A12, Figure A13, Figure A14, Figure A15, Figure A16, Figure A17, Figure A18, Figure A19 and Figure A20). In Figure 5, Figure 6, Figure 7 and Figure 8, we present the SDFs for the year 2100 using seven bio-climate environmental predictors for all three included climate scenarios (IPSL, MRI, and MIROC). Similar to the above, here we present an overview of the Tien Shan–Alai region. For details about Figures below, see the Discussion Section 4.4. The top predictors and the corresponding contribution percentages for Figure 5, Figure 6, Figure 7 and Figure 8 can be found in

Table A3. Percent contribution of 7 bio-climatic environmental predictors to the SDF, year 2000, all data combined.

Variable	Percent Contribution	Permutation Importance
Bio11_2000	41	47.6
Bio17_2000	28.8	14.6
Bio7_2000	27.1	29.7
Bio16_2000	1.6	2.4
Bio10_2000	1.3	3.6
Bio12_2000	0.3	2.1
Bio1_2000	0	0

,

Table A4. Percent contribution of 7 bio-climatic environmental predictors to the SDF, IPSL scenario, year 2100, all data combined.

Variable	Percent contribution	Permutation importance
Bio17_IPSL2100	35.8	14.1
Bio11_IPSL2100	28.6	31.1
Bio7_IPSL2100	20.7	24.6
Bio16_IPSL2100	7.3	2.7
Bio10_IPSL2100	6.9	2.8
Bio12_IPSL2100	0.4	1.9
Bio1_IPSL2100	0.3	22.9

,

Table A5. Percent contribution of 7 bio-climatic environmental predictors to the SDF, MRI scenario, year 2100, all data combined.

Variable	Percent contribution	Permutation importance
Bio7_MRI2100	37	27.7
Bio11_MRI2100	29.1	41.3
Bio17_MRI2100	27.5	5.4
Bio10_MRI2100	2.9	9.5
Bio12_MRI2100	2.3	9
Bio16_MRI2100	1.1	7
Bio1_MRI2100	0	0

and

Table A6. Percent contribution of 7 bio-climatic environmental predictors to the SDF, MIROC scenario, year 2100, all data combined.

Variable	Percent contribution	Permutation importance
Bio17_MIROC2100	27.9	7.6
Bio7_MIROC2100	24.2	13.7
Bio11_MIROC2100	23	50.6
Bio1_MIROC2100	19.2	15.9
Bio12_MIROC2100	2.2	4.2
Bio10_MIROC2100	2.2	1.4
Bio16_MIROC2100	1.2	6.7

.

The reported Area Under the Curve (AUC) values for these SDFs were 0.982, 0.973, 0.982, and 0.984 for Figure 5, Figure 6, Figure 7 and Figure 8, respectively.

3.2. Habitat Suitability

Based on the contribution percentage of each country in the SDM and SDF models seen in

Table 2. Lynx habitat distribution percentage by country.

Turkestan Lynx	Region	Country	All Data Merged 132 Predictors Contribution Percentage	SDF Current Contribution Percentage	SDF IPSL Sum Contribution Percentage	SDF MRI Sum Contribution Percentage	SDF MIROC Sum Contribution Percentage
	Central Asia	Kazakhstan	28.5	5.63	3.66	6.8	4.85
		Kyrgyzstan	21.71	1.21	0.73	1.1	1.36
		China	21.28	9.49	6.58	4.1	8.23
		Uzbekistan	2.16	0.37	0.21	0.25	0.28
		Tajikistan	3.33	0.84	0.42	0.54	0.64
		Turkmenistan	0.16	0.07	0.02	0.05	0.04
		Afghanistan	0.03	0.77	0.49	0.35	0.53
	South Asia	India	0.09	0.4	0.08	0.17	0.35
		Pakistan	0.08	0.74	0.21	0.25	0.67
		Nepal	0.04	0.03	0	0	0.01
		Bhutan	0.01	0	0	0	0
Countries with other lynx species and subspecies			19.59	80.45	87.6	86.39	83.04

below, we evaluated which countries provide the most favorable habitats for the lynx’s present and predicted future distribution. Table 2 was created using the “Zonal statistics” tool in QGIS.

It can be seen that Kazakhstan, Kyrgyzstan and China have the highest contribution in lynx distribution in both SDM and SDF models. In the case of using 132 predictors, models showed that Kazakhstan has the greatest percentage (28.5%) of suitable lynx habitat, followed by Kyrgyzstan (21.71%) and China (21.28%). The countries mentioned have the Northern Tien Shan and Zhetisu Alatau mountains with suitable vegetation cover, water access, and an abundance of smaller prey species. The Western Tien Shan mountains, shared between Kazakhstan, Kyrgyzstan, Uzbekistan and the northmost part of Tajikistan, also provide a favorable habitat for the lynx populations, which also indicates the contribution of Uzbekistan and Tajikistan (2.16% and 3.33%, respectively). According to the models’ estimation, the contribution percentage of other countries in Central and South Asia differ insignificantly and provide less habitat area for the felid. The mentioned row of countries with other species and subspecies was based on models automatically running the calculations for suitable habitats throughout the world in all of the other countries, including other continents, and is not relevant for our evaluation. For the evaluation of predicted lynx movements (SDF) based on seven bio-climatic predictors, see the Discussion Section 4.4 and Section 4.5.

4. Discussion

4.1. Data Classification

For the standardization process and usage of data obtained on lynx sightings over a large area, it is important to collect and classify the data on lynx at the local scale according to their reliability. First proposed by Molinari-Jobin et al. [92], the method to classify the data in three categories was approved and applied by a number of researchers studying lynx in mountainous conditions [91,100,101]. For this research, we also used three sets of data according to their reliability—C1 (Confirmed data), C2 (Probable occurrence) and C3 (Unconfirmed data). After evaluating the models created on these three sets of data, we noted the accuracy of all three categories of data collected laid over models. Moreover, the models of C3 data in the area of Kazakhstan and Kyrgyzstan, shown in Figure A8 and Figure A17, seem to be highly precise and even more accurate compared to models of C2 data, illustrated in Figure A7 and Figure A13, in the aforementioned countries. Notedly, models calculated on the C2 and C3 data do not provide the distribution information of lynx in China. As we did not obtain the C2 and C3 data from Xinjiang and the country has a different set of environmental predictors classification, the extrapolation did not cover the Xinjiang Tien Shan area for our distribution models. The accurate models created based on combined data provide the necessary insight into the predicted current and future lynx habitat distribution and dispersal.

4.2. Notes on the Turkestan Lynx’s Taxonomic Status

For the study area of this manuscript, we excluded the Altai region (N 48°45′–51°55′, E 83°05′–90°35′), located in the border region of Kazakhstan, Russia, Mongolia and China. This is due to the unclear taxonomical position of the Turkestan and Altai lynx (Lynx lynx wardi Lydekker, 1904) and where the border between the two subspecies is located in the region. According to our own observations, the differences in physical characteristics of the lynx inhabiting the Altai mountains and Zhetisu Alatau–Tien Shan mountains are distinct, with the Altai lynx being larger in size than the Turkestan lynx. Still, we included the Saur (N 47°12′, E 85°05′) and Tarbagatai (N 47°14′, E 82°13′) mountain system in the SDM and SDF maps, as the Saur–Tarbagatai lynx population has less distinguishable traits from the Turkestan lynx. Yet, some researchers believe the Saur–Tarbagatai population can still be identified as the Altai lynx, as the felid in this region frequently migrates to the Altai mountains (oral commination by Starikov S.). Thus, it is legal in Kazakhstan to hunt the lynx in the hunting grounds of the Saur mountains.

It should be noted that the subspecies differentiation between the Turkestan and Altai lynx has never been specially studied. At the same time, there is an assumption about the possibility of the Altai subspecies being identical or very close to the Turkestan lynx [29], while some researchers do not classify the Altai lynx as a separate subspecies and also consider this felid to be the Turkestan lynx [102]. For conservation purposes, it is essential that we conduct special studies in this direction. In order to determine the position of the Turkestan lynx as a subspecies and possibly strengthen its conservation status, we have conducted comparative morphometric research on lynx skulls and started molecular–genetic research, the results of which might help us to gain insight into the taxonomic position of the lynx in the region.

4.3. Analysis of Present Distribution

Considering the increasing ecological challenges for the Turkestan lynx populations such as habitat loss and fragmentation, mostly caused by overgrazing, deforestation, prey base depletion, infrastructure development, unregulated tourism, poaching, and retaliatory killings by herders [69,103], the evaluation of the elusive felid’s distribution and suitable habitats holds essential value for the creation of conservation strategies. Furthermore, assessment of the potential current lynx distribution and predicted dispersal according to a changing climate with the subsequent protection measures will also lead to the conservation of important habitats for other plant and animal species. In the northwest part of Turkestan lynx’s global range, where there has been no prior specialized and large-scale research conducted, we used our data and survey occurrence records of various lynx populations and created the first highly predictive habitat distribution models for the lynx in the region. For the current and predicted distribution forecast models, we relied on 132 environmental and 7 bio-climatic predictors, respectively. Thus, the models based on the predictors mentioned, illustrated in Figure 4 and Figure 5, provide an ostensive characterization of present favorable lynx habitat that will accelerate the application of conservation management to areas crucial to the Turkestan lynx.

Figure 4 shows the model based on all lynx observation records obtained, combined with 132 environmental predictors, where we can visually see that lynx sightings fall within regions that our model confirms to be highly potential habitats. In particular, these regions include the Northern Tien Shan and Zhetisu Alatau areas in Kazakhstan, and the Western Tien Shan area in the border of Kazakhstan and Uzbekistan.

Notwithstanding the lack of spatial data from other Tien Shan ridges in Kyrgyzstan and China and their absence from the Western Tien Shan and Alai ridges in Tajikistan due to the fact that we did not collect data from this country, our model suggests that the distribution of lynx populations could include suitable areas in those regions as well. In particular, in Kyrgyzstan, the model presents the most propitious habitats to be located in the Ile Alatau and Kungei Alatau north to the Issyk Kul in the Northern Tien Shan, and especially in the Terskei Alatau east and south of Issyk Kul between the Northern and Inner Tien Shan, as well as in the Talas Alatau, Chatkal and Koksu ranges in the Western Tien Shan and the Fergana range between Western and Inner Tien Shan. In China, the model shows that the most highly favorable habitat is located in the Kalajun area of Northern Tien Shan in Xinjiang, near the Tekes River’s mouth. In Tajikistan, the average level of predicted occurrence is estimated to be in the Alai range at the upper stream of Zeravshan River, near the Komarou Nature Refuge.

4.4. Analysis of Predicted Movements (SDF)

The figures of the SDF models, based on seven bio-climatic predictors, indicate that the highest predicted distribution index of the Turkestan lynx, according to the current distribution, can be found in the following areas: the Ile Alatau, Kungei Alatau, Terskei Alatau and Uzynkara ranges in the Almaty region (Kazakhstan); in the Ile Alatau and Terskei Alatau ranges in the Issyk-Kul region (Kyrgyzstan); in the Ugam, Pskem and Chatkal ranges in the Tashkent region (Uzbekistan) and Jalal-Abad region (Kyrgyzstan); and Talas Alatau range in the Turkistan region (Kazakhstan). Additional hotspots can be found in the central part of Tajikistan, the northeastern part of Afghanistan, and the far north of Pakistan. Inspecting Figure 6, it can be identified that the described hotspots seem to increase towards the south-western distribution (near Tajikistan and the far western part of Tien Shan), and also near Almaty. However, the distribution seems to decrease in the northern parts (near Zhetisu Alatau (Kazakhstan)) and in the very southern parts (near Pakistan). Inspecting Figure 7, and comparing it to Figure 5, it can be identified that the described hotspots seem to change following similar patterns observed for the IPSL scenario (Figure 6), with the difference being that there is a significant increase observed near the area of Taldykorgan (Kazakhstan), which is not visible for the IPSL scenario. Lastly, for the MIROC scenario, it seems that the distribution hotspot near Almaty increases and densifies in the far western part of Tien Shan. Additionally, the southern distribution seems to decrease.

Overall, it can be observed that following the IPSL scenario, the Turkestan lynx distribution is expected to densify, move northwards when following the MRI scenario, and again densify for the MIROC scenario. Details of the observed shifts in the predicted distribution ranges (likelihood of occurrence) can be seen in

Table 3. The movement trends of lynx according to the three global climate model scenarios.

Scenario	Move Northwards	Move Southwards	Move towards Higher Altitudes	Move towards Lower Altitudes	Predicted Distribution Range Decreased	Predicted Distribution Range Increased
MIROC	NO (it centralizes)	NO (it centralizes)	YES	YES (especially in the southern range)	YES	YES (only around the Northern Tien Shan Mountain range)
MRI	YES	NO	YES	YES (but not in mountainous regions towards central Kazakhstan)	YES (in southern areas)	YES (in northern areas)
IPSL	NO	YES	YES (in southern areas)	YES (in northern areas)	YES (in northern areas)	YES (in southern areas)

below.

According to our SDF models based on seven bio-climatic predictors, shown in Table 2, the greatest contributions to the lynx populations are presented by China (9.49%), followed by Kazakhstan (5.63%) and Kyrgyzstan (1.21%). Kazakhstan seems to be the only country with the possible expansion of suitable habitats in the case of MRI—a possible temperature increase of approximately two degrees—but the opposite is observed with Kyrgyzstan in the case of the cooling MIROC scenario. With over 93% of the country occupied by mountainous areas [104], Tajikistan is estimated to provide 0.84% of suitable habitats for the lynx, which reduces insignificantly in cases of possible climate change. Afghanistan and Pakistan, with the presence of the Pamir and Hindu Kush mountains, follow with a similar current contribution of 0.77% and 0.74%, respectively.

Overall, in all the countries, it is observed that the lynx habitats might significantly shrink in size in the case of the IPSL scenario (significant warming scenario). Considering the calculations of the most possible MRI scenario, the predicted climate change might result in habitat reduction in all countries except Kazakhstan. This, as presented in Table 3, proves that the lynx moves northwards in the case of the MRI scenario, and we are likely to observe a significant distribution range decline in several lynx habitats with the future climate change. It must be noted that these models are “minimum-estimates”, meaning that with the low number of future predictors available, here we present a possible approximation of the future distribution range using solely seven climate predictors. We encourage future research to develop more future predictors for more robust future Species Distribution Models.

4.5. Habitat Suitability, Protected Areas and Conservation

The high percentage of favorable habitat in Kazakhstan, shown in Table 2, can be explained by a larger area of mountains in the northwest range considered for research. More precisely, the country has the largest part of the Northern Tien Shan and Zhetisu Alatau mountains. As these mountains continue in the east direction into China, the models suggest that the ridges in Xinjiang present a highly favorable environment for the lynx population and their distribution, which can be especially seen in the SDF columns of Table 2. In Kyrgyzstan, despite the lack of data used, the presence of Northern Tien Shan ridges, and a major part of the Western and Inner Tien Shan mountains, allowed the models to estimate that the country provides necessary lynx habitat; this can also be applied for the Western Tien Shan in Uzbekistan, with a high predicted lynx occurrence index in the country.

It is clear that the higher lynx occurrence indices in the study area depend on a number of environmental conditions, such as climate, land cover (coniferous forest and shrubs), access to a water source (the river and lake network), limited human disturbance, terrain (elevation and the level of annual snowfall, as well as the level of snow line), abundance of prey, and the presence of protected areas in these lynx habitats (see Figure A1, Figure A2, Figure A3, Figure A4 and Figure A5 and

Table A2. Percent contribution of 132 environmental predictors to the SDM, all data combined.

Variable	Percent Contribution	Permutation Importance
WorldSoil2	41.9	44.5
WorldSlope1	15.8	0.4
World_MIN_RH_MAR	7.1	0
World_MAX_RH_OCT	5.4	0
World_MAX_RH_SEP	4.1	0
WCaltitude	3.6	5.1
srad6	3.1	0.2
WorldProtectedAreasMerged4	3	4.8
tmax2	1.6	2.3
World_MAX_RH_DEC	1.5	0
tavg11	1.3	0
World_MIN_RH_OCT	1.3	0
Prec09	1.2	5.5
srad8	1.1	1
tmax3	0.8	0.5
BIO1_2_5min	0.7	0
BIO7_2_5min	0.7	0.5
tmax12	0.4	0
tavg3	0.4	0.1
srad12	0.4	0
FAOCC1	0.4	4.1
tmax11	0.4	0.5
tavg10	0.4	0
Prec12	0.4	3
tavg1	0.3	0.3
BIO14_2_5min	0.3	0.6
tmin11	0.3	0
tmax10	0.3	0.1
srad11	0.2	4.8
srad7	0.2	0
GlobalBigRivers11	0.2	0.7
BIO4_2_5min	0.2	0
tmax1	0.1	0.4
World_MIN_RH_JUL	0.1	0.2
World_MIN_RH_JAN	0.1	0.4
World_MAX_RH_APR	0.1	0
BIO5_2_5min	0.1	0
FFNov2020_3	0.1	4.8
VE4	0.1	0.7
BIO2_2_5min	0.1	0.5
GlobalLakes2	0.1	0.5

, Table A3, Table A4, Table A5 and Table A6 of Appendix A). Figure A1, Figure A2, Figure A3, Figure A4 and Figure A5 additionally represent the omission of the models. Table A2, Table A3, Table A4, Table A5 and Table A6 additionally represent the contribution percentage of the corresponding model, representing the “importance” of all included environmental predictors for building the model (retrieved from the HTML MaxEnt output).

The abundance of prey plays an important role in lynx stability, especially during winter [30], while extensive hunting for the prey may also result in a decline in the lynx population. During our research, we observed that in the Northern Tien Shan, during the years where the prey species, in particular, tolai hare, were declining in number, a similar tendency was noted in the lynx. In the habitats where the prey species are abundant, the lynx populations seem to be more stable, as well (

Table 4. Independent captures (IC) and index of the average abundance per 100 camera trap days (IC per season—SIC) of the lynx and its main prey species in the study area, 2013–2022.

Region	Trap Days	Species
		Turkestan Lynx		Tolai Hare		Roe Deer		Siberian Ibex		Red Tree Squirrel		Wild Boar		Red Deer
		IC	SIC	IC	SIC	IC	SIC	IC	SIC	IC	SIC	IC	SIC	IC	SIC
Northern Tien Shan	11,479	67	0.58	105	0.91	260	2.26	539	4.69	101	0.88	162	1.41	228	1.98
Western Tien Shan	353	1	0.28	-	-	3	0.85	75	21.24	-	-	25	7.08	9	2.55
Zhetisu Alatau	182	6	3.29	-	-	-	-	40	21.97	-	-	-	-	-	-
TOTAL	12,014	74	0.61	105	0.87	263	2.19	654	5.44	101	0.84	187	1.55	237	1.97

).

In the Northern Tien Shan, where most of our camera traps were installed, the lynx was more frequently registered by camera traps in the Kungei Alatau range compared to the Ile Alatau (on average, 0.74 and 0.44 lynx individuals per 100 trap days in Kungei Alatau and Ile Alatau ranges, respectively). In Kungei Alatau, during the period of our research, we noted a higher abundance of the tolai hare, the main lynx prey in the Tien Shan (on average, 2.34 hares per 100 trap days). In Ile Alatau, the Siberian roe deer, another key lynx prey, was frequently caught on camera (on average, 3.33 individuals per 100 trap days). In the Western Tien Shan and Zhetisu Alatau, where the camera traps were stationed for a shorter period of time, we noticed the lynx to be more frequent in the habitats with more roe deer, tolai hare (for the Western Tien Shan) and mountain hare (for the Zhetisu Alatau) present. Thus, in the territory with a relatively high abundance of prey and a smaller number of disturbance factors, a higher abundance of lynx can be traced.

Higher lynx occurrences can be noticed in habitats with protected areas within all of our study area. The level of general population stability can especially depend on the level of protection these PAs provide. In order to evaluate the contribution made by PAs for habitat suitability, we calculated the percentage of PAs laid over habitats our models estimated to be suitable for the lynx populations in the south-east of Kazakhstan, as well as the X-fold increase in predicted lynx occurrence indices of the PAs compared to the rest of Kazakhstan (as seen in

Table 5. Lynx habitat distribution percentage (HDP) between the PAs in the south-east of Kazakhstan and Occurrence Presence Factor (OPF) within all of Kazakhstan (X-fold increase in lynx presence compared to the rest of the country).

Regions	Protected Areas	SDM 132 Predictors		SDF Current		SDF IPSL		SDF MRI		SDF MIROC
		HDP	OPF	HDP	OPF	HDP	OPF	HDP	OPF	HDP	OPF
Northern Tien Shan	Ile-Alatau SNNP	15.18	36.38	19.94	15.31	23.25	13.39	18.55	8.63	19.11	16.3
	Almaty SNR	4.57	15.65	5.89	6.66	6.33	5.44	4.79	3.10	6.42	8.33
	Kolsai Kolderi SNNP	12.90	19.98	17.44	8.91	20.07	7.80	16.55	4.84	17.26	10.12
	Sharyn SNNP	7.97	21.20	2.92	2.56	2.81	1.88	3.27	1.64	2.47	2.49
Western Tien Shan	Aksu-Zhabagly SNR	8.26	16.91	7.83	5.29	10.91	5.60	12.40	4.80	10.13	7.85
	Sairam-Ugam SNNP	9.63	18.87	13.58	8.78	18.35	9.02	15.60	5.78	13.65	10.13
	Karatau SNR	0.42	3.05	1.79	4.26	1.16	2.09	1.99	2.71	0.84	2.29
Zhetisu Alatau	Altyn Emel SNNP	12.52	18.97	4.61	2.30	4.40	1.67	7.50	2.15	3.61	2.07
	Zhongar Alatau SNNP	28.54	18.04	25.99	5.42	12.72	2.02	19.36	2.31	26.50	6.35

below). The calculations were made using the “Zonal statistics” tool in QGIS.

In Kazakhstan, the model based on the 132 predictors suggests that the highly favorable habitats can be found in the Northern Tien Shan at the Ile-Alatau and Kolsai Kolderi National Parks and in the Zhetisu Alatau at the Zhongar Alatau and Altyn Emel National Parks. With the calculations of SDF models based on seven predictors, it can be seen that the Ile-Alatau, Kolsai Kolderi, and Zhongar Alatau National Parks and the Sairam-Ugam National Park in the Western Tien Shan play a significant role in conserving the lynx habitats in future climate change scenarios. In the case of Nature Reserves, where the level of protection is the highest, Almaty Nature Reserve provides the key habitats to the lynx populations despite having much less territory compared to the other National Parks in the Northern Tien Shan.

In China, the possible lynx habitats are protected in the Tomur Peak National Nature Reserve, Kalajun Provincal Park, Kuerdening National Nature Reserve, and Bayinbuluke National Nature Park, with our models estimating the most suitable habitats being located in the Kalajun Provincal Park area. In Kyrgyzstan, the models show that the relevant lynx habitats are located at and protected by the Chong-Kemin National Park on the Ile Alatau and Kungei Alatau ranges of Northern Tien Shan, the Besh-Aral Reserve on the Chatkal Valley of Western Tien Shan, and the Issyk-Kul and Naryn Reserves and the Saymaluu-Tash National Park on the Fergana range in the Inner Tien Shan; in the Karatal-Japyryk Nature Reserve in the Inner Tien Shan the lynx number seems to slightly increase over the years [56]. In Uzbekistan, the Ugam Chatkal Nature Park, Chatkal Nature Reserve and Ugam Chatkal Reserve in the Western Tien Shan maintain the habitats for the lynx and its prey species [67]. In Tajikistan, SDM shows the occurrences being mostly near the Komarou Nature Refuge.

As the threats to the lynx populations are still present and increasing in the mountains, the key solution for conservation of the felid and other animal and plant species, as well as their habitats, is to create new PAs in the areas where they are most necessary, as well as expand the territory of the existing PAs. In Kazakhstan, one of key countries for the lynx habitats, we recommend creation of Uzynkara Nature Reserve in the Uzynkara range of Northern Tien Shan, while in this and Terskei Alatau range there is no PA present yet. As the Ile-Alatau National Park is located in the vicinity and partially in the territory of Almaty megapolis, it is also important to increase the size of the park and strengthen its protection system.

5. Conclusions

In summary, our results describe the present distribution and future dispersal for the year 2100 of Turkestan lynx at the northernmost part of their global geographical range—in the Northern and Western Tien Shan and Zhetisu Alatau mountains. The highly predictive distribution models provide the data on suitable habitats in our study area, as well as within all 11 countries of its distribution. According to the SDM models on 132 environmental predictors and our subsequent calculations, we noted the most suitable lynx habitats to be located in Kazakhstan, followed by Kyrgyzstan and China. Based on seven bio-climatic environmental predictors, we estimate a lynx habitat reduction in most predicted future climate scenarios.

With the significant proportion of suitable habitat relevant to the lynx population located in the northwestern part of its global range, it is evident that the Northern and Western Tien Shan and Zhetisu Alatau mountains are crucial for lynx population stability and future survival. For the conservation of the lynx and other rare animals, it is essential to create new and expand existing protected areas, make a collaboration between the countries to study this felid, and create ecological corridors between the mountain regions and countries contributing to lynx distribution in Central and South Asia.