The Pilbara region of northwestern Australia includes the Hamersley and Chichester ranges and features rugged terrain, drainage channels, and alluvial fans (Department of Water 2016; Rojas et al. 2018). Subterranean habitats are common in mesas and consolidated rock formations, especially where channel iron deposits (CIDs) fill ancient river channels, creating vugs and cavities through iron enrichment and mineral replacement, ideal for stygofaunal habitats (Callan and Gunawardene 2017; Clark et al. 2021).

We selected an 11,560 km² area within the Hamersley Range between the Robe River and Fortescue River drainages (Fig. 1). This area encompasses a high density of sampling sites primarily accessed by exploration bore holes drilled by the mining industry. The number of presences of amphipods in these samples further guided study area location. Sample sites are concentrated along exposed stratigraphy of the Mount Bruce Supergroup (see Fig. 1), including the Fortescue, Hamersley and Turee Creek Groups, overlying the Archaean granite-greenstone terrane of the Pilbara Craton (Rojas et al. 2018).

The study area experiences a semi-arid climate influenced by the Australian monsoon and tropical cyclones, resulting in high summer temperatures and variable rainfall (Rojas et al. 2018). High evaporation rates, approaching 100% under extreme conditions, limit runoff and recharge in this region, with approximately 60% of precipitation lost to evapotranspiration and 40% contributing to groundwater (Skrzypek et al. 2016, 2019).

Recharge occurs mainly through surface runoff and stream flow during heavy rainfall events from cyclones and monsoonal low-pressure systems (Skrzypek et al. 2019; Clark et al. 2021), especially near active drainage systems (Department of Water 2016; Rojas et al. 2018). In the northwest Hamersley, aquifers recharge via infiltration along valley margins, downward leakage via faults and fractures, sheet flow across alluvial fans, and direct rainfall (Department of Water 2016; Winter et al. 1998). Recharge efficiency is influenced by local hydrogeology, including clay layers in alluvium and shallow water tables that can restrict infiltration (Department of Water 2016).

Amphipod sampling in the Pilbara typically occurs near or within sedimentary siliciclastic rocks (Fig. 1), such as CID and alluvial formations, due to reliance on mining boreholes. These formations also support stygofauna habitats (Clark et al. 2021). The majority of sample sites we selected for this study are located in sedimentary rock types within the Brockman Iron Formation or the Wittenoom Formation (Fig. 1). Amphipod occurrence data at these sites was sourced from Framenau et al. (2021). Within our study area, the dataset comprised 5776 samples collected between 1899 and 2016.

Following filtering (see below), each sample site was spatially linked to a set of 122 environmental predictor variables with a pixel size of 0.0025° (Table A1) to characterise local climatic, geological, soil, hydrological and terrain conditions relevant to species distributions. This specific spatial resolution was chosen because it represents the highest and most comprehensive common scale available for the suite of environmental variables covering the study area.

Climate variables (see Table A1), sourced from Williams et al. (2010) and Gibson et al. (2015), were chosen because they are likely to influence amphipod distributions, mainly through their influence on groundwater recharge, temperature, and chemistry, which shape habitat stability and resource availability (Saccò et al. 2019). Terrain variables, extracted from analyses done by Williams et al. (2010) and used by Gibson et al. (2015), are likely to influence hydrological connectivity and habitat heterogeneity, which in turn, shape amphipod distributions (Halse et al. 2014). Soil variables, sourced from The Atlas of Australian Soils (Northcote et al. 1960), available through the Australian Soil Resource Information System (Australian Soil Resource Information System 2024), are likely to influence nutrient availability and substrate properties that regulate above-ground primary productivity (Mokany et al. 2019).

Geological variables and linear structures mapped at 1:2,500,000 were sourced from Department of Energy, Mines, Industry Regulation and Safety (2024). Linear structures, including faults, fractures, and dykes, were included in the analyses as they are geological discontinuities that significantly impact rock permeability and groundwater flow paths. We derived variables of these line features using the “Distance to the Nearest Hub” function in QGIS. These include the minimum distance to the nearest respective line feature in kilometres. A second set of derived variables, “Dyke,” or “Hydro” (Table A1) identifies the categorical description of the nearest dyke or hydrographic feature to each sampling point. Additionally, major streamlines were sourced from Department of Water and Environmental Regulation [DWER] (2007).

We assigned presence (1) if the target taxa was ever found at a given borehole, and absence (0) if it was never found. Thus, overlapping sampling records (i.e., presences or absences recorded at the same site that were repeatedly sampled over different time periods) were removed from this analysis to minimise spatial autocorrelation and sampling bias (Valavi et al. 2021).

Variables that showed clear raster artifacts within the study area, such as abrupt spikes in values, irregular frequency distributions, or systematic striping patterns visible in the raster grid when inspected in QGIS, were excluded, as these features reflect resolution inconsistencies rather than true environmental variation Kuhn and Johnson 2013). Variables with near-zero variance were also excluded to improve model reliability (Kuhn and Johnson 2013).

Using the function “findCorrelation()” from the caret package (Kuhn et al. 2007), we removed highly correlated variables (Pearson correlation coefficient > 0.7 or <-0.7) while retaining one representative from each correlated group, to ensure statistical independence (Aliverti et al. 2021). For each highly correlated pair, the variable with the larger average absolute correlation with all other variables was removed, as it was considered more redundant.

The original database contained 122 variables and 5,776 presence and absence records from multiple sampling events, including repeated observations at the same locations. After applying the filtering process described in Fig. 2, the final modelling dataset was reduced to 999 sampling sites and 31 variables. Of these sites, Amphipoda were present at 180, Paramelitidae at 166, and Paramelitidae sp. B08 at 67.

We first evaluated Generalized Linear Models, Linear Discriminant Analysis, Neural Networks, and Random Forest (RF) models using Accuracy and Kappa Coefficient metrics. A detailed comparison is provided in Suppl. material 1. RF was selected as our classification method based on its superior evaluation metrics. RF algorithm was used as a binary classification model to predict the probability of species presence or absence based on environmental variables.

RF is an ensemble machine learning method that builds many decision trees, each trained on a bootstrap sample, that is, a new dataset created by randomly sampling the original data rows with replacement (Kuhn and Johnson 2013). During training, each tree learns a set of decision rules by repeatedly splitting the data into groups based on the predictor variables, and during prediction, every tree casts a vote for presence or absence; the majority vote determines the final classification (Breiman 2001).

Key RF parameters are the number of trees and predictor variables considered at each split; tuning them manages model complexity and improves predictive performance (Kuhn and Johnson 2013). Bootstrap sampling and random variable selection reduce overfitting and improve stability, enabling RF to model complex, non-linear patterns (Breiman 2001).

To ensure an unbiased evaluation of model performance, the dataset was divided into three parts following the recommendations of Kassambara (2018) and Kuhn and Johnson (2013). The first part is the training set (50% of the dataset). This is the portion of the data used to fit the model, to allow the algorithm to learn the relationships between environmental variables and species presence or absence. The second part is the testing set (30% of the dataset), which is used during model development to assess how well the model performs on data it has not yet seen and to help tune model parameters. Finally, the validation set (20% of the dataset) provides an independent dataset used only at the end of the process to evaluate the model’s ability to generalise to completely new observations.

We developed RF models for the three target taxonomic groups, each trained on an unbalanced dataset (i.e., a dataset where the number of presence records does not match the number of absence records). Such an imbalance can affect how well the model learns to detect the minority class (Kuhn and Johnson 2013). Model performance was evaluated using the receiver operating characteristic (ROC) curve.

ROC analysis assesses the balance between sensitivity (true positive rate) and specificity (true negative rate) across classification thresholds (Sammut and Webb 2011). The threshold is the probability cut-off for classifying sites as “presence” or “absence” (Hoo et al. 2017). Sensitivity is the proportion of true presences correctly identified; specificity, the proportion of true absences (Kuhn and Johnson 2013; Kassambara 2018); see Suppl. material 2.

To improve model performance, we tuned the “mtry” parameter of the RF, which controls how many predictor variables are randomly tested at each split of a decision tree (Kuhn et al. 2007). This procedure does not remove or select variables from the dataset; instead, all predictor variables are used since each tree has many splits, and the forest consists of many trees. Therefore, every variable has the opportunity to be used throughout the forest, even if “mtry” is small. We only adjust how many of them the model is allowed to consider at each split. We tested “mtry” values ranging from 2 to 26 (in steps of two) to identify the value that produced the best predictive performance (Kuhn et al. 2007). Model training was performed by averaging several runs (10) to achieve better performance (Barbet-Massin et al. 2012). The models were built using 1000 trees. These tuning steps were implemented using the “train” function of the “caret” package (v6.0-92) (Kuhn et al. 2007).

Because our dataset has many more absences than presences, during model training, we used stratified sampling to ensure that presence and absence classes are represented proportionally using the “train” function in the “caret” package (v6.0-92) (Kuhn et al. 2007) and the “sampsize” parameter from the “randomForest” package (Liaw and Wiener 2002).

In addition, we used balanced accuracy and Matthews correlation coefficient (MCC) to assess the performance of the models. The balanced accuracy metric calculates the average of sensitivity (true positive rate) and specificity (true negative rate), making it appropriate for unbalanced datasets where one class (e.g., absence) dominates (Chopra and Dixit 2020; Joloudari et al. 2023), see Suppl. material 2. In general, balanced accuracy values below 0.6 suggest poor performance, while values between 0.6 and 0.7 are considered modest but acceptable in ecological modelling contexts (Allwright 2022).

The MCC metric is a measure of binary classification performance that ranges from −1 to +1. MCC considers all four elements of the confusion matrix (see Suppl. material 2) and produces high scores only when the model performs well across both classes; a positive MCC indicates that predictions are better than random chance (Matthews 1975; Chicco et al. 2021).

To determine the optimal threshold for classifying a case as a predicted presence, we employed two widely used methods based on the ROC curve: Youden’s J statistic and the Closest-Top-Left method. The Youden J statistic (also known as Youden’s index) identifies the threshold that maximises the difference between sensitivity and specificity, thereby balancing the trade-off between true positive and false positive rates (Youden 1969). The Closest-Top-Left method selects the threshold that minimises the Euclidean distance from the point (0,1) on the ROC curve, prioritising an optimal balance by reducing both false positive and false negative errors (Jones and Athanasiou 2005).

We assessed variable importance by ranking each variable’s contribution to the ability of the model to distinguish species presence from absence. Breiman (2001) explains that RF assesses variable importance by measuring the decrease in model accuracy when the values of a variable are randomly shuffled. One common approach described by Breiman (2001) is the Mean Decrease in Accuracy (MDA), which indicates how much the model relies on each variable, with larger drops in accuracy indicating greater importance.

Variable importance for each taxonomic group was calculated using the function “varImp” of the “caret” package (version 7.0.1) (Kuhn et al. 2007), based on the MDA scores computed by the “randomForest” package (see Fig. 3). We defined the argument “scale = TRUE” in the function to scale the importance scores so that the maximum importance score is set to 100, and all other scores are scaled relative to this maximum (Kuhn et al. 2007).

Partial dependence plots (PDPs) help make complex machine-learning models, such as RF, easier to understand by showing how a single variable influences model predictions (Molnar 2022). They display the marginal effect of one variable while averaging out the effects of all others, making it easier to see how changes in that variable affect the predicted outcome (ŷ) (Greenwell 2017). PDPs were performed using the package “pdp” (Greenwell 2017) to examine the contribution of the five top-ranked numerical variables for each taxa in predicting presence. To enhance interpretation, PDPs are presented alongside histograms that display the distribution of the values of each variable.

The histograms highlight the values where each variable is most frequently observed, with higher frequencies indicating regions of greater data coverage. Such regions enable the model to make more confident predictions, as they are supported by robust data (Kuhn and Johnson 2013). We also used the “rug()” function of the “graphics” package (R Core Team 2024) to mark along the x-axis the locations of the data points in the given variable. This gives a visual sense of how different values of the variables are distributed along the axis, helping to interpret the PDP in relation to the actual data distribution.

To create a predicted probability map for our model, we generated a grid of points spaced 750 metres apart using the “Regular Points” function in QGIS version 3.34.3 (QGIS Development Team 2023). We applied the “Point Sampling Tool” to extract variable information at each grid location. We calculated the distance from each grid point to nearby linear structures and identified structure types using the “distance to nearest hub” tool.

Some categories within the hydrology, rock type, and dyke variables are present in the grid of points used to generate the predicted probability map but are not represented in the training subset used to build the model. Since the model cannot make predictions for categories it has not encountered during training, grid points with these unseen categories were excluded. This exclusion resulted in gaps in the predicted probability maps (Fig. 7), corresponding to areas with missing categorical coverage in the training set.

To show spatial variation in prediction consistency across taxonomic levels, we calculated the standard deviation of predicted probabilities at each grid cell for species-, family-, and order-level models. We then plotted this grid together with the true-positive locations predicted by each model using the validation dataset (Fig. 8). Higher standard deviation values indicate greater disagreement among taxonomic resolutions predicted probabilities.

The R code and tables used for building SDMs in this study are available on Figshare (https://doi.org/10.6084/m9.figshare.29192723). Table A2 presents the ODMAP (Overview, Data, Model, Assessment, and Prediction) protocol (Zurell et al. 2020).

For clarity, the variables are grouped into five categories: climate, geology, hydrology, soil, and terrain, with their respective contributions from each category to the respective models shown in Fig. 3.

Terrain variables contributed 52.2%, 44%, and 41.4% to the Amphipoda, Paramelitidae, and Paramelitidae sp. B08 models, respectively (Fig. 3), with the standard deviation of the percent slope ranking high for all taxonomic groups, and slope ranking second for Amphipoda, and first for Paramelitidae sp. B08 (Fig. 3).

The contribution of geological variables was higher for Paramelitidae sp. B08 (22%) and lower in the Amphipoda and Paramelitidae model (15.6% v. 21%). Minimum distance to a dyke ranked among the top five variables across all taxonomic levels. In the Paramelitidae sp. B08 model, distance to the nearest fold ranked third, while distance to the nearest fault was not present as in the other two models (Fig. 3).

Contribution of climate variables is similar across models, with precipitation of the coldest quarter ranking high in all (Fig. 3). Proximity to major tributaries ranks low but appears in the top 20 in all models. Proximity to minor rivers ranks within the top 20 for the Paramelitidae and Paramelitidae sp. B08 models (Fig. 3). The low relative importance of soil variables in influencing the probability of presence (Fig. 3) is not surprising given that amphipods live in groundwater (Gibson 2018).

Figs 4, 5, 6 indicate that areas with less variation in percent slope tend to have higher predicted outcomes. This suggests that amphipods occur in areas with relatively flat terrain with minimal variation in slope. This trend is particularly apparent for all amphipoda, with the PDP for terrain ruggedness (Fig. 4) showing higher predicted values at lower ruggedness and as percent slope decreases. Predicted outcomes peak at a distance of 10 km from a dyke for the Amphipoda and Paramelitidae models (Figs 4, 5).

The coefficient of variation for monthly diurnal temperature range reaches its highest predicted outcome at low values (around 9.2) for the Amphipoda and Paramelitidae models (Figs 4, 5), suggesting a preference for stable diurnal temperatures. All models show higher predicted outcome for less than 57 mm of precipitation during the coldest quarter (Figs 5, 6, A1), but this is also where the data are most concentrated.

Analysis of observed and predicted species presence, along with key predictors for each taxonomic group, indicates that most confirmed and predicted occurrences of amphipods are located in areas characterized by low slope values and low standard deviation of percent slope (Fig. 7).

The northwestern part of the predicted probability maps (Fig. 7) corresponds to the sedimentary rock areas of the Marra Mamba and Wittenoom Formations (Fig. 1), which form a narrow, diagonally oriented band extending toward the southeastern part of the study area and continuing northeastward. This region, characterised by low variability of the monthly diurnal temperature and percent slope, coincides with valley fill alluvium, CID and paleovalley areas (Rojas et al. 2018).

In the northwestern area (Fig. 7), the Amphipoda and Paramelitidae models predict a high probability of presence (>0.5), while Paramelitidae sp. B08 model predicts values near zero. The Paramelitidae sp. B08 model shows a slightly higher probability in the southwestern portion of the map, overlapping with high values from the Amphipoda and Paramelitidae models. Low standard deviation values in this region (Fig. 8) indicate strong agreement among the models’ predicted probabilities.

In the eastern side of the maps (Fig. 7), high predicted probability of presence for all taxonomic groups align with the Fortescue River catchment (Fig. 1). The predicted probabilities of presence increase from the eastern side of the map toward the bottom-centre, where most true positives are concentrated. The lower standard deviation in this region suggests it meets the shared ecological requirements of each taxon (Fig. 8). This area is aligned with minor streamlines near the Queens Valley, Firetail, and Kings Valley mines, extending south toward Hamersley Road.

Terrain variables accounted for roughly half of the total importance in all models, making topography the strongest predictor for all taxonomic groups of amphipods. For Amphipoda, climate variables were the second-largest contributor after terrain, whereas for Paramelitidae, and Paramelitidae sp. B08, it was the geological variables. Hydrological and soil variables made only minor contributions in all cases. Perina et al. (2024) suggests that stygofauna distributions are shaped more by paleohydrological patterns rather than current river catchments.

Mouron et al. (2022) reported comparable results, with slope making a secondary contribution to explained variance in amphipod models, while the highest contributions were associated with water properties, and geological variables not included in our study due to substantial missing data.

Our quantitative results align with qualitative observations that stygofauna are commonly associated with paleochannels or paleodrainage valleys in the Pilbara (Humphreys 2008; Rojas et al. 2018). The higher probability of amphipod occurrence in low-relief terrains also aligns with findings from Mokany et al. (2019), who reported species richness of stygofauna to be highest in the valleys where groundwater is close to the surface. These relatively low-relief areas may act as groundwater recharge hotspots or provide stable aquifer structures that support amphipod populations (Halse et al. 2014; Hose et al. 2015; Perina et al. 2024). Flatter terrains may allow more water to percolate, recharge the aquifer, and increase the input of nutrients to subterranean ecosystems (Perina et al. 2024), whereas steeper slopes may lead to faster runoff, reducing the amount of water that infiltrates the ground. Valley bottoms accumulate sediments, forming floodplains or fans that act as hydrological buffers and improve the connectivity of the catchment (Gallant and Dowling 2003).

Distance to the geological linear features of dykes, folds and faults are all included in the top 20 ranked variables influencing occurrence of amphipods. Intrusive dykes compartmentalise the fractured-rock aquifers of the Hamersley Range by cutting across the stratigraphy and disrupting groundwater connectivity, creating semi-isolated hydrogeological blocks that can restrict stygofauna dispersal and generate distinct species assemblages across adjacent gullies (Perina et al. 2024).

However, a dyke may also function as an aquifer if it meets certain criteria (Gupta et al. 2012). It must consist of deeply weathered rock (e.g., dolerite), contain dense fracture networks to boost permeability, align with the hydraulic gradient, and possess structures that allow groundwater storage and transmission.

In folded terrains, aquifers become interlayered with variable dips, while large-displacement faults juxtapose units of contrasting hydraulic properties, generating both permeable conduits and barrier zones (Ball et al. 2014). Faults, which often occur together with folds, enhance the downward movement of water and promote more intense weathering (Perina et al. 2024). These linear structures modify aquifer connectivity, recharge pathways, and the scale of groundwater flow systems. They frequently separate valley aquifers from the recharge that comes from areas of high elevation and transfer flow control from topographic gradients to systems regulated by permeability (Ball et al. 2014).

As such, it may be that these geological features influence the occurrence of amphipods in complex ways by affecting the influx and flow of water from above and below the ground (Rojas et al. 2018) as well as nutrient and oxygen levels (Maurice and Bloomfield 2012; Marmonier et al. 2023). As mentioned above, this level of complexity in the three dimensional context requires further investigation using geophysical techniques to assess the dynamics of subsurface water. For example, Kirsch (2006) explains that electrical resistivity measurements from geophysics can be used to estimate water content and porosity, which then allows calculation of percolation time.

Boucher et al. (2009) have also used electrical resistivity tomography, resistivity logging, time-domain electromagnetic, and magnetic resonance sounding to assess how water moves through an unconfined sandstone aquifer in semi-arid southwestern Niger. Their study provided direct estimates of water content, regional groundwater recharge, and highlighted zones of enhanced infiltration and focused drainage beneath gullies and alluvial fans. However, this level of detail is only likely to be achievable in a more local context and remains a challenge for understanding the distribution of stygofgauna at regional scales.

It is important to note that the modelled climate variables used in this study do not accurately capture local underground conditions. As Roghanchi et al. (2017) points out, subterranean climates are influenced by site-specific factors such as geology, rock properties, groundwater flow, and elevation, factors that can vary significantly even over short distances. While Mammola and Leroy (2018) also highlights the limitations of bioclimatic variables in representing subterranean systems, they note that the temperature of surrounding rock tends to equilibrate with the mean temperature of infiltrating water. This suggests that external climatic predictors can still offer a useful approximation of underground microclimates.

The negative relationship between amphipod presence and increasing variation in the monthly diurnal temperature range may indicate a narrow thermal tolerance consistent with the stable thermal regime of groundwater systems (Hose et al. 2015). Because temperatures below ≈ 10 m remain effectively constant in the absence of sunlight, stygofauna evolve under highly stable thermal conditions and may be particularly sensitive to environments where subsurface temperatures fluctuate (Hose et al. 2015). This pattern may also reflect their preference for flatter terrain (Eberhard et al. 2005), where surface and subsurface temperatures tend to vary less than on steep slopes that are shaded for part of the day.

The negative relationship of fauna presence with increasing winter rainfall may indicate that a period of low groundwater recharge is needed to stabilise trophic conditions following the summer months when rainfall events are likely to be highly sporadic and extreme (Rojas et al. 2018; Saccò et al. 2019).

In our spatial models, the highest predicted probabilities for all fauna followed major valleys in the Marra Mamba and Wittenoom formations and in the catchment area of Serpentine Creek, Duck Creek, Caves Creek, the Robe River, and the Fortescue River. The southeastern portion of the study area within the Brockman Iron Formation shows low standard deviation and consistently high predicted probabilities across taxonomic levels, likely indicating environmentally optimal conditions and a core area of habitat suitability.

These low-relief valley floors and gullies, where recharge is focused, and groundwater is shallower, receive greater inputs of oxygen and nutrients than adjacent elevated ridges, and likely provide more favourable habitat for stygofauna (Halse et al. 2014; Hose et al. 2015; Perina et al. 2024).

The high standard deviation in the northwestern and northeastern areas within the Marra Mamba and Wittenoom formations likely results from missing presence records for Paramelitidae sp. B08. While presence is consistently predicted at the family and order levels, this pattern may indicate either a sampling gap or the presence of another paramelitid species with distinct ecological requirements; this discrepancy warrants further investigation.

The high predicted probability areas also match the geological and hydrological context described by Rojas et al. (2018), where high-permeability, mineralised banded iron formations in the Marra Mamba Iron Formation host spatially restricted aquifers sustained by local groundwater recharge, but capable of producing significant yields under dewatering conditions. In such hydraulically constrained systems, groundwater abstraction and dewatering may disproportionately affect subterranean habitats and fauna (Perina et al. 2024).

Rojas et al. (2018) describes groundwater in the Wittenoom Formation as primarily recharged by downward leakage from overlying inland alluvial and CID/pisolite aquifers associated with paleovalley systems in the Hamersley Basin. Pisolite deposits often provide suitable subterranean fauna habitat due to their extensive cavities and secondary porosity (Perina et al. 2024). These sediments are hydraulically connected to the underlying fractured and karstified dolomites, facilitating vertical groundwater transfer. As a result, the Wittenoom Formation is known to host significant groundwater resources, with discharge largely driven by mining dewatering and water supply operations (Rojas et al. 2018).

A key limitation of this study is that our models do not account for physical or geographical barriers. This leads to an assumption of unlimited dispersal for the target species, overestimating accessible habitat. Future work should include data on physical dispersal barriers or simulate spatial dispersal limitations (Leroy et al. 2016).

Another key limitation of this study is the clustering of sample locations and the coarse taxonomic resolution of specimen identifications (Halse et al. 2014; Mokany et al. 2019). Order-level models in our study may be affected by uneven family representation, as Paramelitidae is over-represented in the sample. This can bias model predictions toward the environmental associations of this over-represented family and may not accurately reflect the ecological diversity of the entire taxonomic group (Kuhn and Johnson 2013). Incorporating higher-resolution taxonomic data enhance model precision and ecological inference.

Key aquifer properties like water chemistry, oxygen levels, and salinity influence species distribution at fine scales (Humphreys 2008), but are excluded here due to limited informative data at sampled boreholes. While these coarse-resolution surface variables successfully predicted stygofauna distribution, they may not fully capture all local-scale sub-surface attributes driving these communities.

Future refinements to improve model performance could include integrating high taxonomic resolution through eDNA-based identification (Saccò et al. 2022b) and expanding the environmental variable set to increase discriminatory power (Kuhn and Johnson 2013). Our findings are consistent with those of Seeland et al. (2019), who noted that while higher taxonomic levels may share characteristics that facilitate classification, the loss of species-specific information can negatively affect predictive performance. Similarly, the species-level model in our analysis achieved higher sensitivity and balanced accuracy than the family- and order-level models. By addressing current limitations, future SDMs can provide more robust tools for understanding and conserving subterranean biodiversity in the face of environmental change.