The study area is located on the northern and southern of the Tianshan Mountain range in Xinjiang, China, covering the central section of the Tianshan Mountain and the surrounding regions on both sides. This area is a key habitat for Marmota baibacina (Figure 1). The Tianshan Mountain range runs through the central part of Xinjiang, dividing the study area into northern and southern regions, forming an ecological transition zone with distinct dry and wet conditions (Wang et al., 2021). The northern lies at the southern edge of the Tianshan Mountain and the Junggar Basin. The climate is temperate and humid, with high annual precipitation (200–800 mm). Forests, steppe, and alpine meadows are mainly distributed here, making it a region rich in water resources and biodiversity (Zhang et al., 2024). Areas such as the Ili River Valley and the Bortala River Valley are known as the “Jiangnan of the North,” with favorable natural conditions. The southern of the Tianshan Mountain is characterized by a typical arid continental climate, with low annual precipitation, dominated by deserts, gobi, and oasis. The northern of the Tianshan Mountains includes administrative regions such as Ili Kazakh Autonomous Prefecture, Bortala Mongolian Autonomous Prefecture, Urumqi, Changji Hui Autonomous Prefecture, and Turpan, which are economically developed and are key agricultural, pastoral, and industrial areas in Xinjiang (Fang et al., 2024). The southern of the Tianshan Mountains includes parts of Bayingolin Mongolian Autonomous Prefecture and Aksu Prefecture, with an economy based on oasis agriculture and resource development, with Korla City as the core city on the southern. Overall, the northern and southern of the Tianshan Mountains have significant differences in natural geography and ecological environments, making them important regions for studying climate change, biodiversity conservation, and ecosystem services (Wang et al., 2024). Marmota baibacina is widely distributed in the alpine meadows, steppe, and surrounding mountainous areas of the Tianshan Mountains, primarily inhabiting the high-altitude regions between 2,500 and 4,000 meters above sea level (Sun et al., 2018). It relies on the unique ecological and climatic conditions of the area for survival (Du et al., 2022), as shown in Figure 1.

In this study, we collected data on the natural distribution of 190 Marmota baibacina during a field survey conducted during the Third Xinjiang Comprehensive Expedition in 2022-2023, in which the precision, latitude and habitat characteristics of each sample site were recorded in detail. To ensure data accuracy, duplicate points within 1 km were removed using the SDM toolbox in Spatial analysis tool to avoid spatial autocorrelation, and records with ambiguous coordinates were excluded, and subsequent screening yielded a final set of 111 sample points for Marmota baibacina (Figure 1). Spatial thinning was performed to ensure the minimum distance between any two occurrence points was ≥1 km. This threshold was determined based on two considerations: (a) Habitat characteristics of Marmota baibacina: As a colonial rodent with a home range of 0.8–1.2 km² per colony, a 1 km minimum distance avoids over-representing a single colony and reduces sampling bias; (b) Model resolution consistency: The threshold matches the 1 km spatial resolution of environmental variables, ensuring each occurrence point corresponds to an independent environmental grid and avoids pseudo-replication. The thinning algorithm randomly retained one point within each 1 km buffer until no points violated the minimum distance constraint.

Considering previous studies (Koshkina et al., 2020; An et al., 2023) and the habitat characteristics of Marmota baibacina, 29 environmental variables were selected potentially influencing its distribution (see Table 1). These variables include historical climate data obtained from Worldclim (https://worldclim.org), providing 30 arc-second (approximately 1 km) data for the period 1970-2000. While field surveys were conducted in 2022-2023, historical climate data (1970-2000) were used to characterize long-term habitat suitability, as short-term climate fluctuations (e.g., recent warming) may not yet fully reflect species distribution shifts. Future projections (2041-2100) account for ongoing climate changes. Specifically, the environmental variables considered are as follows:

(1) Terrain variables: slope, slope direction, and elevation data sourced from the geospatial data cloud (http://www.gscloud.cn/). (2) Soil variables: two soil variables crucial for Marmota baibacina habitat—soil pH value (Soil_ph) and soil symbol (Soil_symbol)—extracted from global soil pH data and China’s 1:400,000 soil symbol map compiled by the Nanjing Institute of Soil Science, Chinese Academy of Sciences. The soil map includes 72 soil classes and 247 subclasses, with Xinjiang covering 23 soil symbols, such as thin layer soil, glacier, alluvial soil, water body, calcareous gray soil, calcareous soil, chestnut soil, dune quicksand, sandy soil, impinged soil, anthropogenic soil, black soil, gley soil, embryonic soil, saline soil, denatured soil, gypsum soil, alkaline earth, salt works, urban industrial and mining areas, loose lithologic soil, and leaky rock. (3) Geomor: geomorphic data specific to Xinjiang (https://www.geodata.cn/). (4) River: distance from rivers, centered on river channels, with a 1 km buffer (https://ngcc.cn/). (5) NDVI (Normalized Difference Vegetation Index): maximum NDVI dataset for China from 2000 to 2020, processed on the Google Earth Engine platform using Landsat5/7/8 remote sensing data. This dataset, with a spatial resolution of 30 meters and annual temporal resolution, involved cloud and shadow removal, followed by NDVI extraction using linear interpolation and S-G smoothing methods. (6) Vegetation: Xinjiang vegetation symbol data and grass symbol data (https://www.geodata.cn/). These variables were selected based on their relevance to the habitat preferences and ecological requirements of Marmota baibacina, with the aim of providing a comprehensive analysis of its distribution patterns in the Tianshan Mountains of Xinjiang.

Future climate variable data for this study were derived from CHELSA CMIP6 scenario data, a high-resolution climate database for global land surface areas (https://chelsa-climate.org/) (Karger et al., 2017). Furthermore, the CMIP6 scenarios represent the most advanced generation of climate projections, offering improvements over the previous CMIP5 framework. The SSP scenarios integrate socioeconomic narratives with emission pathways, providing a more realistic and comprehensive basis for modeling future climate impacts compared to the Representative Concentration Pathways (RCPs) used in CMIP5. This makes CMIP6 the current state-of-the-art for assessing biodiversity responses to climate change. The dataset has a spatial resolution of 30 arc-seconds (approximately 1 km × 1 km) and it includes bioclimatic variables projected for two future time periods: 2041–2070 and 2071-2100.These future climate scenarios are based on the Shared Socio-Economic Pathways (SSP), which represent different socio-economic development trajectories and their implications for greenhouse gas emissions and climate change impacts: SSP126: represents a sustainable development pathway with low greenhouse gas emissions, aiming for a radiative forcing of 2.6 W/m² by 2100. SSP370: represents a middle to high-end emission scenario, leading to moderate levels of warming with a radiative forcing of 7.0 W/m² by 2100. SSP585: represents a high-emission scenario with extensive fossil fuel use, resulting in high greenhouse gas concentrations and a radiative forcing of 8.5 W/m² by 2100 (Li et al., 2025). These climate scenarios provide essential data for assessing potential future impacts on ecological systems, including habitat suitability for species like Marmota baibacina in Xinjiang, considering varying levels of climate change and greenhouse gas emissions.

Spatial analysis tool was used to process the original data, including the removal of invalid data, filling in missing values, standardizing the data, and unifying its resolution and projection (WGS1984) for subsequent analysis and modeling (Guisan et al., 2017). For the environmental raster data, the preprocessing included: (1) masking all layers to the unified study area extent; (2) converting all layers to a consistent spatial resolution of 30 arc-seconds (1 km) and the WGS 1984 geographic coordinate system using the bilinear resampling method (for continuous variables) or the nearest neighbor method (for categorical variables); and (3) ensuring no cells contained NoData values within the study area mask. This process resulted in a harmonized and analysis-ready dataset. The resampled environmental variables were then converted to ASCII format using the SDM Toolbox v2.5 extension tool (An et al., 2023).

Given the differences in input sample data between machine learning models and the MaxEnt model, and to avoid the influence of variable correlation on prediction results, the screening process for the 29 environmental variables was systematically optimized: all models first performed basic variable screening, where for machine learning models, the 29 environmental variables in species presence-absence data were extracted, with variables of low importance removed after calculating their importance scores (An et al., 2023); for the MaxEnt model, the 29 environmental variables and sample points were input into the model for preliminary computation, with variables and sample points showing zero contribution eliminated, ensuring that variables entering subsequent analyses had basic ecological relevance (Koshkina et al., 2020). Following initial screening, machine learning models further tested for multicollinearity using the Variance Inflation Factor (VIF) - multicollinearity refers to high linear correlation between two or more independent variables in a regression model, which may cause unstable regression coefficients and reduced model interpretability (Guisan et al., 2017). With a threshold of VIF = 10, variables with VIF>10 were excluded, leaving 10 core variables.

For the MaxEnt model, after initial screening, variables were refined through Principal Component Analysis (PCA) and importance evaluation: the PCA tool analyzed remaining variables, and if the absolute correlation coefficient between two variables exceeded 0.8, the one with lower contribution was removed to reduce information redundancy; meanwhile, variable importance was comprehensively assessed using contribution rates, permutation importance, and the Jackknife test, ultimately retaining 16 variables to ensure the model preserved ecologically significant drivers while eliminating redundancy (Phillips et al., 2006).

This study used the R software platform to implement four commonly applied machine learning models (XGBoost, RF, SVM, and LogBoost) and one species distribution model (MaxEnt) to investigate the suitable habitat distribution of Marmota baibacina. The five models selected for this study each have unique advantages in classification and prediction tasks, with different theoretical foundations. RF and XGBoost are powerful ensemble learning algorithms adept at capturing complex nonlinear relationships and handling high-dimensional data (Zhang et al., 2011). SVM is effective in high-dimensional spaces and for cases where the number of dimensions exceeds the number of samples (Li et al., 2017). LogBoost (LogitBoost) is a boosting algorithm designed for classification. MaxEnt is particularly robust for presence-only data, making it a standard in SDM studies (Phillips et al., 2006; Elith et al., 2008). To ensure optimal performance of both the individual machine learning models and the MaxEnt model, the presence and absence point datasets of Marmota baibacina were randomly divided, with 25% allocated for validation and 75% for training (Li et al., 2019). Ten repeated experiments (Logistic format) were conducted, which enhanced the model’s ability to accurately predict the species’ potential range. Finally, all models were evaluated for accuracy using the Receiver Operating Characteristic (ROC) curve and the Area under the Curve (AUC) (Koshkina et al., 2020). The ROC curve is a graphical tool for assessing the performance of binary classification models, while the AUC represents the area under the ROC curve, with values ranging from 0 to 1 (Luo et al., 2017). A higher AUC value indicates better model performance. The impact of environmental variables was comprehensively assessed using the percentage contribution, permutation importance, and Jackknife test from the MaxEnt model (An et al., 2023). The percentage contribution represents the contribution of each climatic variable to the geographic distribution of Marmota baibacina during the training process, while the permutation importance quantifies the decrease in the model’s AUC value when the climatic variables in the training points are randomly replaced (Araújo and New, 2007). The Jackknife method is similar to cross-validation; it involves excluding one or more sample points at a time and calculating a corresponding statistic using the remaining sample points. This method analyzes the importance of each individual variable in constructing the distribution model (Dai et al., 2022).

The classification of model simulation results was performed using the reclassification tool in spatial analysis, employing the Jenks natural breaks classification method (He et al., 2023). The Jenks natural breaks classification method is a technique that minimizes within-class variance and maximizes between-class variance when categorizing data (Zhai and Li, 2012; An et al., 2023). This method ensures that the differences within each category are as small as possible, while the differences between categories are as large as possible. The simulation results were divided into four categories: unsuitable (0-0.2), low suitability (0.2-0.3), moderate suitability (0.3-0.5), and high suitability (>0.5) (Jiang et al., 2022), to determine the potential geographic distribution of Marmota baibacina in the Tianshan region. After reclassification, the number of grids in each category was calculated, and the area of suitable habitat for Marmota baibacina under different climate scenarios was computed.

The AUC results of the models used in this study were all above 0.9, indicating high accuracy (Figure 2). Machine learning models identified 10 key variables for simulations, while MaxEnt selected 16 variables based on feature importance and correlation metrics. To ensure comparability between MaxEnt and machine learning approaches, MaxEnt was executed using both the 10-variable subset and its native 16-variable set. All models achieved AUC values exceeding 0.9, with Random Forest exhibiting the highest AUC, followed by the 16 variables MaxEnt configuration. LogBoost demonstrated the lowest accuracy among the evaluated models.

The potential distribution area of Marmota baibacina is shown in Figure 9. Under the current climate scenario, its potential suitable habitat is primarily distributed across Bayingolin Mongolian Autonomous Prefecture, Hejing County, Ili Kazakh Autonomous Prefecture, and the western part of Bortala Mongolian Autonomous Prefecture. The total suitable habitat area for Marmota baibacina in the Tianshan Mountains of Xinjiang ranges from 2.75×10⁴ km² to 13.59×10⁴ km². The area predicted by the MaxEnt model is the smallest, while the area predicted by RF is the largest, mainly composed of low suitability and moderate suitability areas, which cover 3.55×10⁴ km² and 6.40×10⁴ km², respectively. These areas are mainly distributed in the northern and southern regions of the Tianshan Mountains in Xinjiang, including the Bogda Peak in the northern Tianshan. The area of high suitability predicted by the XGBoost model is the largest, covering 3.88×10⁴ km², and is mainly distributed in Bayingolin Mongolian Autonomous Prefecture HeJing County. Among high-suitability zones predicted by five machine learning models after variable screening, MaxEnt produced the smallest area covering 0.62×10⁴ km²; this value further decreased to 0.42×10⁴ km² when using MaxEnt’s native 16-variable set.

Based on three common socio-economic pathways proposed by the IPCC (SSP1~RCP2.6, SSP3~RCP7.0, SSP5~RCP8.5), the geographic distribution of Marmota baibacina was predicted under three future climate change scenarios for the periods of 2041–2070 and 2071–2100, as shown in Figure 10. A comparison in Table 2 indicates that, under all five models, the area of future potential suitable habitats for Marmota baibacina continues to decrease, although the overall spatial patterns remain highly consistent with the current period.

Under the three future climate scenarios, in the SSP126 scenario, the XGBoost model predicts an increase in low suitability areas and a decrease in high and moderate suitability areas, with the total suitable habitat area reducing by 0.7×10⁴ km² by 2041–2070. The RF model shows an increase in low and moderate suitability areas, with a decrease in high suitability areas, resulting in an overall increase of 2.87×10⁴ km² in suitable habitat area. In contrast, the SVM and MaxEnt models predict a decrease in low, moderate, and high suitability areas, with total suitable habitat areas decreasing by 5.21×10⁴ km², 0.58×10⁴ km², and 0.35×10⁴ km², respectively. The LogBoost model predicts a reduction in low and moderate suitability areas but an increase in high suitability areas, leading to a total increase of 0.53×10⁴ km² in suitable habitat area. By 2071–2100, these trends are somewhat mitigated, and the total suitable habitat area slightly increases.

Under the SSP370 and SSP585 scenarios, these trends intensify further. The XGBoost, SVM, LogBoost, and two MaxEnt models predict a reduction in the total suitable habitat area for Marmota baibacina, with the decrease being particularly pronounced under the high emission scenario (SSP585). High suitability areas shrink further to 0.35×10⁴–3.02×10⁴ km², and the regions where suitable habitats decrease are mainly concentrated in Bayingolin Mongolian Autonomous Prefecture HeJing County. Spatially, the contraction of suitable habitats is most pronounced in southern Tianshan, this migration may compress available grazing lands for livestock, intensify competition between marmots and domestic herbivores, and further degrade already fragile alpine meadows. The reduced habitat area coupled with concentration in specific regions could elevate human exposure to plague-infected fleas, particularly in pastoral zones where herders and livestock frequently interact with marmot burrows.

Based on the predictions of the five models, the suitable distribution area of the Tianshan marmot under current climatic conditions and different future climate scenarios shows significant differences, as illustrated in Table 2. Under current climatic conditions, the RF model predicts the largest total suitable habitat area, reaching 13.59×10⁴ km². In contrast, the MaxEnt model predicts the smallest total suitable habitat area, with only 2.75×10⁴ km², and the high-suitability area is limited to 0.42×10⁴ km². Under future climate scenarios, the predictions of suitable habitat areas exhibit diverse trends across models. In the SSP126 scenario for 2041–2070, the RF model predicts the largest total suitable habitat area, reaching 16.46×10⁴ km², the highest among the five models, including a high-suitability area of 2.54×10⁴ km². Compared to the current climate scenario, the high-suitability area decreases slightly. Conversely, the MaxEnt model predicts the smallest total suitable habitat area at 2.40×10⁴ km², with the high-suitability area further reduced to 0.39×10⁴ km². As emission intensity increases under the SSP585 scenario, the high-suitability area predicted by all models shrinks further. For example, the high-suitability area predicted by the XGBoost model decreases from 3.88×10⁴ km² under current conditions to 2.71×10⁴ km².

In terms of regional distribution, the high-suitability areas under current climatic conditions are mainly concentrated in parts of Bayingolin Mongolian Autonomous Prefecture, Ili Kazakh Autonomous Prefecture, and Bortala Mongolian Autonomous Prefecture. These regions are also where the reduction in suitable habitats is most pronounced under future climate conditions. Particularly under the SSP585 scenario, the high-suitability areas in these regions decrease sharply.

The AUC values of the models used in this study were all above 0.9, indicating a high level of accuracy, but significant performance differences were observed across the models. The RF model predicted the largest total suitable habitat area, attributable to its ensemble structure that aggregates predictions from multiple decision trees through bootstrap aggregation (bagging). This approach effectively captures complex nonlinear interactions among environmental variables while mitigating overfitting, thereby maximizing habitat inclusivity. In contrast, the XGBoost model, being highly sensitive to parameter settings and suitable for complex variable structures, performs slightly worse than the RF model. The SVM model, which performs well on high-dimensional data, requires extensive parameter tuning, which led to its slightly inferior performance in this study. The LogBoost model performed relatively poorly, possibly due to its dependence on the distribution of training samples, which caused overfitting and reduced its ability to fit nonlinear relationships. The models predicted the distribution of Marmota baibacina under current and future conditions, and the results indicate that climate change will have a significant impact on the species’ habitat area and this finding were consistent with (An et al., 2023). A comparative analysis revealed that using the 16 environmental variables screened by the MaxEnt model alone yielded slightly higher prediction accuracy than the RF model, followed by the XGBoost and SVM models. The LogBoost model showed the poorest performance among all (Zhang et al., 2011; Luo et al., 2017). This hierarchy in performance can be briefly attributed to the respective algorithms’ strengths: MaxEnt’s probabilistic framework is particularly adept at handling presence-only data and capturing niche boundaries, which may explain its superior accuracy with its own variable set (Elith et al., 2008). RF, while robust and capable of modeling complex interactions, might be slightly less optimized for this specific data structure than the specially tuned MaxEnt (Zhang et al., 2011). The poorer performance of LogBoost could stem from its higher sensitivity to noisy data and potential overfitting in this ecological modeling context. The LogBoost model showed the poorest performance among all (Li et al., 2019). In terms of predicted suitable habitat area, the RF model estimated the largest total suitable area, while the two MaxEnt models predicted the smallest total suitable area. However, further validation indicated that despite its superior statistical accuracy, the MaxEnt model yielded conservative predictions of suitable habitats, deviating from field survey data (Guo et al., 2023). In contrast, the RF model demonstrated excellent agreement with actual sampling data, particularly in accurately delineating highly suitable areas on the northern of the Tianshan Mountains (Koshkina et al., 2020; An et al., 2023). Therefore, considering both model accuracy and ecological plausibility, the RF model was identified as the optimal predictive tool.

Ecologically, MaxEnt’s focus on slope (2°–8°) and NDVI (0.1–0.55) aligns with field observations of marmots’ preference for gentle terrain and moderate vegetation cover—critical for burrow stability and predator detection (Koshkina et al., 2020). In contrast, machine learning models emphasized climatic variables (Bio18, Bio2, Bio11, Bio15) contributing >70% to predictions, reflecting the dominance of seasonal water and temperature dynamics in shaping broad-scale marmot distribution. This divergence is not a limitation but a strength: MaxEnt captures fine-scale microhabitat features that machine learning models may overlook, while machine learning models identify macroclimatic constraints that presence-only models cannot fully resolve (Araújo and New, 2007).

The study considered 29 environmental variables, and after selection, 10 key variables with higher contributions to the five models were chosen, as well as 16 key environmental variables for the MaxEnt model. Machine learning models (RF, XGBoost) prioritized climatic variables (Bio18, Bio2, Bio11, Bio15), collectively contributing >70% to predictions, reflecting the dominance of seasonal water and temperature dynamics in alpine ecosystems (Zhang et al., 2011). In contrast, MaxEnt emphasized slope, NDVI, and Bio18, indicating that presence-only models rely more on proximate habitat features like vegetation productivity and microtopography (An et al., 2023). This divergence underscores the impact of data type: models requiring presence-absence data excel at detecting broad climatic constraints, while presence-only models depend on fine-scale environmental proxies (Araújo and New, 2007). Notably, Bio18 emerged as a universal driver, with optimal values at 150–180 mm, aligning with field observations of Marmota baibacina’s preference for mesic alpine meadows (Koshkina et al., 2020). Soil symbol and slope (2°- 8°) were critical in MaxEnt, suggesting that edaphic properties and terrain gradient are essential for modeling species presence in data-scarce environments (Wang et al., 2024). In contrast, the study by Liao et al. (2020) used fewer key variables, focusing primarily on the effects of temperature and precipitation. This study not only included these climatic variables but also comprehensively considered factors such as symbol, vegetation, and topography, making the model more comprehensive and the predictions more precise (Zhao et al., 2023). Furthermore, although Li et al. (2022) also considered different climate scenarios (e.g., RCP2.6 and RCP6.0), their time frame only extended to 2080. In this study, more climate scenarios (e.g., SSP126, SSP370, SSP585) and a longer time span (up to 2100) were included, providing a deeper insight into the far-reaching effects of climate change on the future distribution of Marmota baibacina (Wang et al., 2021). This detailed approach not only improved the spatial accuracy of the model but also further demonstrated the long-term impact of climate change on species survival over time, ensuring the comprehensiveness and objectivity of the analysis (Li et al., 2025). Existing studies on the potential distribution of Marmota baibacina show that future climate change will significantly affect the species’ suitable habitat area, with a decreasing trend (An et al., 2023). This study’s findings are consistent with previous research results (Wang et al., 2024). However, an increasing number of studies are now focusing on ensemble modeling of individual models. Thuiller et al. (2009) proposed the first computational platform framework, BIOMOD, which integrates multiple species distribution models to improve classification accuracy and precision. Ensemble models, through weighted averaging or voting mechanisms, mitigate single-model biases and enhance robustness (Araújo and New, 2007; Marmion et al., 2009). For instance, Banda et al. (2024) significantly improved the accuracy of endangered Apalis flavigularis distribution predictions by integrating MaxEnt, RF, and GLM. Our focus on standalone models aimed to rigorously assess algorithm-specific suitability for alpine species habitat modeling, rather than optimizing predictive accuracy. Future studies could build on our multi-model comparison to explore Bayesian averaging or dynamic weighting strategies, reconciling discrepancies to enhance ecological plausibility.

The machine learning models and the MaxEnt model differ in the selection of variables. The MaxEnt model assumes the independence of variables and directly analyzes the contribution of environmental factors (Phillips et al., 2006), while machine learning models are based on data training and assess the interactions between variables through feature importance (Peters et al., 2007; Zhang et al., 2011). Future research could further integrate the advantages of both approaches to improve prediction accuracy and spatial adaptability. The dominant variables driving Marmota baibacina’s distribution—identified via machine learning models (10 variables) and MaxEnt (16 variables)—are ecologically meaningful, supported by both habitat requirements of the species and field-validated evidence, rather than mere repetition of model results. Bio18 (150–180 mm optimal) is the universal primary driver (42.8% contribution in MaxEnt; >25% cumulative in machine learning models) because it directly controls soil moisture—an essential factor for alpine grass growth (the core food source of Marmota baibacina) (Koshkina et al., 2020). Field surveys confirm that mesic meadows (Bio18: 150–180 mm) support 2–3 times higher marmot burrow density than arid steppes (Bio18 < 100 mm) (An et al., 2023). Bio11 (coldest quarter temperature < -15°C) ensures sufficient hibernation duration (6–8 months at high elevations), a physiological adaptation to conserve energy in alpine cold (Wang et al., 2024). Bio15 > 55 in MaxEnt; >80 in machine learning reflects stable moisture supply, avoiding extreme droughts/floods that disrupt foraging and burrow stability (Stenseth et al., 2008).

Field observations revealed that, in addition to climate variables, soil also play a role in species distribution. Therefore, these environmental variables were included in the model construction for this study. Vegetation provides essential food and habitat for Marmota baibacina, influencing their nutritional status and ability to avoid predators (Wang et al., 2024). As vegetation height and coverage decrease, Marmota baibacina become better able to detect predators, thus reducing predation risks (Chen et al., 2017). Soil type (alpine meadow/steppe soils preferred) and slope (2°–8° optimal) are critical in MaxEnt due to their roles in burrow construction and foraging efficiency. Alpine meadow soils have high organic matter content (10–15% higher than desert soils) and good drainage, reducing burrow collapse risk (An et al., 2023). Gentle slopes (2°–8°) balance two needs: steep slopes increase burrow instability, while flat areas have dense shrub cover that obstructs predator detection. NDVI (0.1–0.55 optimal) serves as a proxy for vegetation productivity—values <0.1 indicate sparse forage, and >0.55 indicate excessive shrubs, both unsuitable for marmots. Compared with other studies Liao et al. (2020), this study includes more environmental variables, such as Soil symbol and topography, making the models more comprehensive and the predictions more accurate. From the results of all five models, it is evident that Bio18 is the primary factor affecting the distribution of Marmota baibacina. The species thrives in areas with dense vegetation, such as alpine steppe and alpine meadows, where precipitation directly impacts soil moisture and vegetation growth. According to the response curves, Bio18 values between 150 and 180 mm are more suitable for Marmota baibacina’s survival. This is consistent with the study’s findings that Marmota baibacina is most likely to be found in alpine cold desert soils, alpine meadow soils, chestnut-calcareous soils, and steppe soils, with the highest distribution probability at slopes between 2° and 8° and NDVI values between 0.1 and 0.55. According to Li et al. (2022), vegetation and Soil symbol also play key roles in marmot habitat selection, consistent with the findings in this study. The potential distribution of Marmota baibacina is not only influenced by climate, topography, soil, and vegetation but also by human activities. Human infrastructure, such as villages and roads, causes wildlife to avoid these areas, potentially negatively impacting their behavior, reproduction, or survival (Bergström and Skarpe, 1999). The construction of buildings and roads can lead to soil compaction and reduced vegetation coverage (Chen et al., 2017). Therefore, the predicted suitable habitats in this study may be overestimated. Future research should integrate climate, biological and human activity factors to improve the accuracy of habitat predictions.

The suitable habitat of Marmota baibacina is primarily concentrated in mid- to high-altitude areas, ranging from 1,500 meters to above 3,500 meters, which is consistent with the distribution point coordinates collected in the field, indicating a good model simulation. Under both current and future climate scenarios, the potential suitable habitat of Marmota baibacina is mainly distributed in the central part of the Tianshan Mountains in Xinjiang, including regions such as Bayingolin Mongolian Autonomous Prefecture HeJing County, Ili Kazakh Autonomous Prefecture, and the western part of Bortala Mongolian Autonomous Prefecture. These areas show a higher probability of Marmota baibacina’s presence.

According to the predictions under future climate scenarios, climate change will lead to significant ecological changes. Rising temperatures and changes in precipitation patterns will cause both the expansion and contraction of habitats (Jiang et al., 2023; Kang et al., 2023). Previous studies have shown that, by the mid-21st century, precipitation in Xinjiang is expected to increase by 10% to 25%, with temperatures rising by 1.5°C to 2°C; by the end of the century, precipitation could increase by more than 25%, and temperatures may rise by 4°C to 6°C (Wang et al., 2021). This warming shortens Marmota baibacina’s hibernation period by 10–14 days (An et al., 2023), reducing energy storage for reproduction and lowering juvenile survival by 22%. Additionally, advanced plant phenology reduces overlap between peak forage quality (grasses, sedges) and the marmot’s active period, limiting fat accumulation before hibernation (Parmesan and Yohe, 2003; Post et al., 2008). While total precipitation in Xinjiang may increase by 10%-25% by mid-century, the intensification of precipitation seasonality poses significant threats (Stenseth et al., 2008). More frequent summer droughts reduce grass cover by 10%-15%, while heavy winter snowfall blocks burrow entrances, causing substantial juvenile mortality (Linné Kausrud et al., 2007). Under the high-emission SSP585 scenario, precipitation in the warmest quarter decreases by 15%-20% in the southern Tianshan, transforming crucial mesic meadows into arid steppes with insufficient forage (Wang et al., 2021).

Climate warming drives progressive upward shifts of alpine vegetation belts at 30–50 meters per decade, compressing Marmota baibacina habitat into increasingly isolated high-elevation patches (Chen et al., 2017). This fragmentation is most severe in the arid southern Tianshan (e.g., Hejing County), where habitat contraction exceeds northern regions by 2–3 times (An et al., 2023). These isolated patches curtail gene flow and intensify competition with livestock for dwindling forage resources, further exacerbating alpine meadow degradation (Prakash and Ghosh, 2012). Our findings, particularly those derived from the MaxEnt model, resonate with other studies on steppe rodents in the region (An et al., 2023), reinforcing that climate change poses substantial threats to species habitats through significant range reduction. The variation in predictions between models reflects their different sensitivities to environmental variable interactions and habitat definitions (Araújo and New, 2007; Zhang et al., 2011). Future research should prioritize ensemble modeling approaches to reconcile these differences and yield more robust forecasts.

This study has three primary limitations. First, the pseudo-absence points generated through random sampling may introduce spatial bias, and future studies could optimize this process using environmental stratification. While our occurrence data were collected during a single-year expedition, the use of long-term climate normals (1970-2000) for modeling mitigates the direct impact of interannual climatic variability on our projections. Nevertheless, multi-year survey data would be beneficial to account for potential population fluctuations and further validate model stability over time. Second, the exclusion of human disturbances (e.g., roads, grazing intensity) from the models may lead to overestimations of high-suitability areas. Third, the static time frame fails to account for lag effects in species migration and vegetation succession, potentially underestimating habitat fragmentation. While ensemble modeling was not adopted in this work, the study establishes a methodological foundation for multi-algorithm comparisons in alpine species research. Future efforts could combine ensemble predictions with dispersal models to dynamically simulate distribution shifts (Beaumont et al., 2016). Additionally, incorporating participatory GIS (PGIS) to merge herder traditional knowledge with spatial data could significantly enhance the practical applicability of predictions.