Mapping threatened Thai bovids provides opportunities for improved conservation 1 outcomes in Asia 2

36 Wild bovids provide important ecosystem services throughout their native range. In 37 Asia, most are threatened with extinction. Five wild bovids remain in Thailand: gaur 38 ( Bos gaurus ), banteng ( Bos javanicus ), wild water buffalo ( Bubalus arnee ), mainland 39 serow ( Capricornis sumatraensis ) and Chinese goral ( Naemorhedus griseus ). However, 40 their populations and habitats have declined substantially and become fragmented. 41 Here, we identified potentially suitable habitat for these five threatened bovids using 42 ecological niche models, first throughout the species entire distribution and second 43 within Thailand, and quantified how much suitable area remains within protected areas. 44 We combined species occurrence data with 28 environmental variables for modelling 45 and used a spatially-restricted Biotic-Abiotic-Mobility framework for two accessible 46 areas: 1) species-specific accessible areas and 2) a single large accessible area. We 47 applied spatially restricted and weighted average ensembles from eight algorithms 48 when generating maps. For B. gaurus and B. javanicus, the best models predicted 49 suitable habitat was mostly within Southeast Asia , with B. gaurus having predicted large 50 areas in Thailand and India. B. arnee suitable habitat was mostly in India. C. 51 sumatraensis suitable habitat was mostly in Thailand and Myanmar. N. griseus was 52 mainly restricted to China. In Thailand, the highest bovid potential richnesses were in 53 the Northern Forest, Western Forest, Eastern Forest and Dong Phayayen-Khao Yai 54 Forest Complexes. We identified unprotected hotspots with >50% of overall suitable 55 habitat located outside protected areas. B. arnee had the smallest proportion of 56 protected habitat (9%). Suitable areas identified in and out protected areas may guide 57 habitat management and conflict mitigation strategies. 58


Introduction
An important task of wildlife research and conservation is to define the distributional ecology of species, and to understand how they relate to the environment, climate and other organisms (Franklin, 2009).Ecological niche models (ENM) are applied to predict the geographic distribution that is suitable for a species by using ecological niche dimensions combined with species' presence data (Soberon & Peterson, 2005).ENM can be approached using the 'Biotic-Abiotic-Mobility' (BAM) framework which considers the relationship between the species' distribution, geographical and climatic factors and explains the influence of factors on predicted habitat suitability (Peterson & Soberón, 2012).Abiotic (A) factors generally determine the potential distribution (or fundamental niche) of a species, and the intersection of abiotic and biotic (B) factors form the realised niche, or the part of this potential distribution where species actually live (Soberón & Nakamura, 2009).Mobility (M) is the area that is accessible by species related to their distribution over periods of time (the 'accessible area'; Barve et al., 2011).Selecting the extent of species' accessible areas, including buffer zones, impacts model prediction results (Anderson & Raza, 2010;Barve et al., 2011).
Wild Bovidae (Mammalia: Artiodactyla) play significant ecological roles in tropical forests and grasslands (Hassanin, 2014).Bovids are grazers and browsers, modifying plant diversity and abundance within ecosystems (Ripple et al., 2015;Romero et al., 2015).Large wild bovids are also the prey of predators such as tigers (Panthera tigris) and leopards (Panthera pardus) (Simcharoen et al., 2018).Throughout Asia, wild bovid populations are threatened by poaching (Gray et al., 2018) and habitat loss (Nguyen, 2009), especially in South to Southeast Asia (Giam & Wilcove, 2012).Natural habitats have been disturbed by free-grazing livestock, which can lead to interbreeding (e.g. between domestic and wild water buffalo, Kaul et al., 2019), increased competition for food and natural resources (Bhandari et al., 2022), and increased risk of disease transmission between wildlife and livestock (Hassell et al., 2017).
In South and Southeast Asia, there are 27 recognised bovid species (IUCN, 2021), of which seven species are listed as vulnerable, five as endangered and three as critically endangered with extinction.Thailand has five bovids species (Gaur Bos gaurus, banteng Bos javanicus, wild water buffalo Bubalus arnee, mainland serow Capricornis sumatraensis and Chinese goral Naemorhedus griseus) are remaining in their natural habitat.These species are distributed in the other countries from South to Southeast Asia (Fig. 1) and also have different suitable habitats.For example, gaur can be found in evergreen forest or grassland and range from India, Nepal, across Southeast Asia to Peninsula Malaysia (Duckworth et al., 2016).C. sumatraensis also has a wide distribution from Nepal to Indonesia through hill forests to shrubland habitats (Phan et al., 2020).Nevertheless, evaluations of their habitat suitability in Thailand and the other countries have been conducted only in some protected areas (Chaiyarat et al., 2019;Pintana & Lakamavichian, 2013), but not at the regional or national level.
Understanding historical and present distributions of threatened species is fundamental for planning and needs to be constructed at suitable scales, yet are typically lacking in Asia.Regional scale modelling provides an overview of potential habitat and benefits conservation planning (Catullo et al., 2008), such as by informing threatened species reintroduction programmes (Bellis et al., 2021) and management through identifying high-quality habitat connectivity and fragmentation (Crooks et al., 2011), along with global biodiversity trend predictions (Araújo et al., 2019).Several studies have predicted habitat suitability for some of these five wild bovids in local areas, including in Malaysia (Ariffin et al., 2021), China (Ding et al., 2018), Indonesia (Rahman et al., 2019), Thailand (Prayoon et al., 2021) and Nepal (Thapa et al., 2020), but habitat suitability studies for larger extents across their distribution are lacking.
Here, we built ENM for five the Thai bovid species: B. gaurus, B. javanicus, B. arnee, C.sumatraensis and N. griseus at two scales: first, at the regional scale throughout the entire distribution and, second, at the country scale in Thailand.We aim to 1) identify the potential distribution for these five species in South to Southeast Asia, and 2) identify conservation areas in their geographical distribution, with a particular focus on Thailand.

Materials and methods
Our workflow consisted of two main processes of data preparation and model building (summarised in Supplementary Fig. S1) that generated habitat suitability maps for all species and accessible areas used.Data preparation consisted of gathering the species occurrence data and environmental data and selecting the accessible areas.Then, the model building consisted of pre-processing, processing and post-processing steps.

Species Occurrence data
We compiled species occurrence data from January 2000 to June 2021 from researchers, government, NGOs (World Wildlife Fund, Wildlife Conservation Society, Freeland, Panthera, Fauna & Flora International, Friends of Wildlife and RIMBA) and open data sources, including GBIF (https://www.gbif.org/)and eMammal (https://emammal.si.edu/).The data coverage by country can be found in Table S2.The occurrence data were collected through observation of animal signs (e.g.footprint and dung) during forest patrols, direct observation during wildlife surveys, camera trapping, and radio-collar signals (Table S3).We excluded occurrence records that fell outside the species-specific accessible area, duplicated records and museum collections.

Accessible areas
The accessible area refers to the parts of the world that have been accessible to species via dispersal over time (Barve et al., 2011).The extent of the accessible area and the inclusion of a buffer zone has an important effect on ENM performance (Anderson & Raza, 2010;Barve et al., 2011).Therefore, we used two accessible area sizes to delimitate our modelling extent (Fig. 1), and compare ENMs with and without spatial restriction (MSDM).To reduce overprediction and make our predictions closer to realised niche estimates, we used an occurrences-based threshold (OBR) method with ensemble models from Mendes et al. (2020) for creating the spatially restricted ENM (hereon MSDM).OBR is an a posteriori method that restricts the suitable areas of our final ensemble models based on presence and the largest nearest neighbour distance among pairs of occurrences.The first larger accessible area (hereon LA) includes most of the Asian continent and its ecoregions, and all species distributions are included as a common extent.We downloaded current IUCN range maps for each species, then overlapped those on ecoregions (Olson et al., 2001), then combined the results with selected ecoregions based on biogeographic knowledge of the species distributions and habitat preference from the literature reviews.For example, gaur (Bos gaurus) habitat typically contains moist evergreen, semi-evergreen, and dry evergreen forests (Steinmetz et al., 2008;Tanasarnpaiboon, 2016), so we included these regions in our accessible areas.The second accessible area was more restricted and cropped based on individual species-specific distributions (hereon SSA) from literature reviews (Table S1), IUCN polygons or 'ranges' (IUCN, 2020) and the terrestrial ecoregions where they occur.Further details on ecoregions included in accessible areas are in Supplementary Table S5.Overall, we built four combinations between two accessible areas with and without MSDM methods for each species, including 1) No MSDM-SSA; 2) No MSDM-LA; 3) MSDM-SSA and 4) MSDM-LA.

Model building
We processed the species occurrence files and environmental datasets in R 4.0.1 (R Core Team, 2020).We developed reproducible ecological niche models with optimized processing times using the ENMTML package (Andrade et al., 2020), following three main steps: 1) pre-processing, 2) processing and 3) post-processing.
In pre-processing, we performed occurrence thinning using 2 times the cell-size (1 km 2 ) (Velazco et al., 2019) to reduce clustering of species records and sampling bias.We used principal component (PC) analysis (PCA) to reduce the collinearity of the predictors.We assigned species' accessible areas to determine the species' distributions using a mask function.We used random sampling to create pseudoabsence background points in a 1:1 ratio with presence points (Barbet-Massin et al., 2012).The occurrence and pseudo-absence data was divided into two sets for fitting the model (75%) and evaluating the fitted models (25%), using the bootstrapping partition method with 10 replications for each algorithm.
In the processing step, eight algorithms were used to build the ENMs, namely: Bioclims (Booth et al., 2014), Generalized Linear Models (McCullagh & Nelder, 1989), Generalized Additive Models (Hastie, 2017), Random Forest (Liaw & Wiener, 2002), Support Vector Machine (Karatzoglou et al., 2004), Maximum Entropy default (Phillips et al., 2006), Maximum Likelihood (Royle et al., 2012) and Bayesian Gaussian Process (Golding, 2014).All models used the default settings from the ENMTML package, which included the functions from different packages (e.g.dismo, maxnet) based on the algorithms that use to fit the models.The data type used for each algorithm is in Supplementary Table S6.
In the post-processing step, we created ensemble models using the weighted average (WMEAN) method based on the True Skill Statistic (TSS) values for building final habitat suitability and binary maps.The benefits of ensemble models are 1) robust decision making (Ahmad et al., 2020); 2) reducing uncertainty (Marmion et al., 2009); and 3) a combination of several models into one model prediction (Kindt, 2018).We used TSS to calculate threshold values to convert habitat suitability maps into binary suitability maps (0 = unsuitable and 1 = suitable).We used TSS and area under the curve (AUC) for evaluating model performance.The TSS threshold is calculated using the maximum summed specificity and sensitivity and is not based on prevalence, where an equal TSS score for given models means similar performance (Allouche et al., 2006).Therefore, we selected the final models from the best TSS of weighted average ensemble models.
For calculating total suitable areas, we summed the suitable area of the best TSS binary map models using the zonal function in the raster R package (Hijmans, 2021).Then, we created a species richness map ranging from 0 (no suitable areas) to 5 (suitable areas for all five species) by summing the best TSS binary maps of each species.

Protected area analyses
We used the binary (suitable and unsuitable) predictions from the best performing model SSA and LA ensembles for calculating the suitable areas inside and outside protected areas.The source for our protected areas map was the World Database of Protected Area (WDPA) (UNEP-WCMC and IUCN, 2021).We classified protected areas based on IUCN protected areas from WPDA into 8 categories including: categories 1 to 6 as IUCN management categories I to VI; category 7 as 'not applicable' which is included 'not reported', 'not applicable' and 'not assigned' protected areas; and category 8 as non-protected areas, which are the remaining areas that have not been classified as IUCN categories 1 to 7 (UNEP-WCMC and IUCN, 2021).Then, we used the zonal function again to calculate overlapping areas between binary map predictions and protected areas for each species.We used 1 km resolution and WGS 1984 projection for all the calculations.
We calculated the percentage of suitable habitat in WDPA polygons using the threshold values from the best TSS binary suitability maps to separate suitable and unsuitable areas.The areas that had suitability values higher than each species' threshold were classified as suitable and those lower as unsuitable.We used the exact_extract function in exactextractr package (Baston et al., 2021) for extracting the suitable areas (values =1) from binary map rasters in each WDPA polygon.Then, we classified each PA into 5 different suitability categories, based on the percentage of suitable habitat in the PA: low suitability (0 -20%); low -medium suitability (>20 -40%); medium suitability (>40 -60%); high suitability (>60 -80%) and very high suitability (≥80%), and selected only the PAs that have the proportion of suitable area larger than species home range in the result.We have provided the code for creating the models in a GitHub repository.
We found that the PCA reduced the 28 environmental variables into 12 PCs that explained 95% of the environmental variance in the variables for the LA models for all species.The PCs for SSA models explained more than 96% of the total variance and the PC number varied by species, comprising 13 PCs (B.arnee), 11 PCs (B.gaurus, C. sumatraensis), and 10 PCs (B.javanicus, N. griseus).The bioclimatic variables were the important variables in all species models.For LA models, the first two axes (PC1 and PC2) have high contributions from the annual mean temperature (bio01), mean temperature of the coldest month (bio06), mean temperature of the driest quarter (bio09) and mean temperature of the warmest quarter (bio10).The first two axes of SSA models showed high positive contributions from mean temperature of the coldest month (B.gaurus), minimum temperature of coldest month (B.javanicus, C. sumatraensis), annual mean temperature (B.arnee, C. sumatraensis), and precipitation of the wettest quarter (N.griseus).We also found that NDVI, elevation, slope and human population density have less effect on explaining the variability for the first two PCs for all species.
The correlations between PCs and individual environmental variables, PC biplots and percentage of explained variance are summarised in Table S7 and Fig. S2.

Ecological niche models
Overall, all ensemble models showed high performance both for TSS and the area under the curve (AUC) with the highest performing models over 0.8 for all species (Table S8).Models with species-specific accessible areas were not always the best performing models, but most ensemble models performed above 0.7 TSS.The habitat suitability prediction maps using the best model ensembles are in Supplementary Fig. sumatraensis is No MSDM-SSA.We found that all species have small predicted suitable habitats.Moreover, all species models predicted less than 50% of the suitable areas inside PAs.
For B. gaurus, the best model predicted the total suitable area was 612,548 km 2 , of which Thailand has the largest area (29%; 177,216 km 2 ).B. javanicus has a similar predicted distribution to B. gaurus, but is more restricted to Southeast Asia.The total suitable area for B. javanicus predicted was 446,075 km 2 (MSDM-LA) and most of the suitable areas were located in mainland Southeast Asia (96%, 429,759 km 2 ).B. arnee has a patchy distribution of suitable habitat in the southern tip of India, Sri Lanka, through the south of Thailand.The best model (MSDM-SSA) predicted the greatest suitable area in India (35%, 213,247 km 2 ), followed by Thailand (135,552 km 2 ) and Myanmar (99,844 km 2 ).C. sumatraensis suitable habitat is predicted throughout the region modelled, down through the south of Thailand -Myanmar border forest patches and in the southwestern parts of Cambodia in the species accessible area.The total predicted suitable area for C. sumatraensis for the best performing model (no MSDM-SSA) was around 290,000 km 2 .Myanmar had the greatest suitable area of around 23% (66,892 km 2 ).N. griseus suitable habitat predictions ranged through northeastern India, Myanmar and northern Thailand to parts of China, with a total suitable area of 169,882 km 2 with the MSDM-LA model.China has the greatest suitable area, with 96% (162,929 km 2 ) of the total suitable area.Thailand has no suitable areas predicted by the best LA models but around 4% (14,847 km 2 from 331,704 km 2 ) of the total suitable area was predicted in Thailand by the best SSA models.Total of the suitable area in km 2 for each species and countries is shown in Fig. 2, and the prediction maps is in Fig. 3. Fig. 2. Total of the suitable area in km 2 for each species and countries.Blue is the species-specific accessible area (SSA) and grey is the large accessible area models (LA) (see details in Supplementary Table S9).presents the map of species numbers for five species using the binary suitability results, ranging from 0 (not suitable area for the five species) to 4 (suitable area for four species), which is the highest value range for the best model predictions.Interactive maps are provided in the supplementary material (link).

Asia
Most suitable habitats in protected areas are located in IUCN category Ia (Strict nature reserve), Ib (Wilderness area) and II (National park) areas for the best TSS models for all species, while IUCN category V (Protected landscape or seascape) has the least.
Overall, more than half of the species' suitable habitat is not under any form of protection defined by the WDPA (Supplementary Table S10, Fig. S6, Fig. S7).The proportion of the suitable area in each WDPA of the best models is shown in Fig. 4; and the proportion of the best models from SSA and LA for each species are presented in Supplementary Fig. S8 and Fig. S9.
Our results predicted B. gaurus has 27% (166,092 km 2 ) of the total suitable area inside PAs with the best model (No MSDM-LA).We identified 179 PAs (from 2132 PAs; Fig. 4) that had a high proportion of suitable area (≥80% Only 2% (3,904 km 2 ) of the suitable area predicted for N. griseus is currently protected in a PA, and 12% were in 'not applicable' PAs, according to the best model (MSDM-LA).
Our best model had highly suitable areas in only 10 of 2132 PAs (Fig. 4).The highest suitable areas were found in Jin Zhai Goa scenic area (75%) in China and Emae bum NP (75%) in Myanmar, followed by Sichuan Giant Panda Sanctuaries (70%), Three Parallel Rivers of Yunnan (60%) and Hkakaborazi National Park in Myanmar (40%).Wiang La and Mae Ping (Fig. 5 and 6).The hotspots for each species can be found in the supplementary material (Fig. S10).We found that the highest percentage of suitable area was mixed deciduous forest for all species, followed by evergreen forest for B.  Western, Dong Payayen-Khao Yai and Eastern forests have suitable areas for 3 -4 species inside PAs and also in the surrounding areas.Overall, most protected areas had two to three species in the south, such as Kaeng Krachan and Khlong Saeng-Khao Sok forest complex.N. griseus suitable areas are restricted only to the northern forests (Lum Nam Pai-Salawim, Mae Ping-Om koi and Srilanna-Khun Tan), which make these PAs suitable for all species (Fig. 6 B1 and B2), but only in a small area compared to total PAs areas.The best model by TSS for N. griseus (MSDM-LA) showed no suitable areas in Thailand, though the SSA model did.

Discussion
We modelled the potential distribution for the five threatened wild bovid species, distributed in South and Southeast Asia, providing high resolution habitat suitability.Our aim was to build predictive models to identify conservation areas, and potential species richness maps in their entire geographical distribution and then with a particular focus on Thailand.Therefore, we used ensembles of models that make better predictions and achieve better performance than any single contributing model, reducing dispersion of predictions.We identified that suitable areas were fragmented and often (50%) located outside PAs.Those suitable areas outside PAs could possibly be managed as corridors or buffer zones to connect currently fragmented bovid populations, thereby enhancing long-term wild bovid conservation success (Karanth, 2016;Penjor et al., 2021), which requires further investigation.
Our study found most suitable areas for gaur were similar to IUCN range assessments (Duckworth et al., 2016) and consistent with studies that have confirmed species presences, such as in Thailand's PAs (Prayoon et al., 2021), Myanmar (Hein et al., 2020) and Western Ghats in southwestern India and Manas WS in the Himalayan foothills (Choudhury, 2002).Our study predicted larger gaur suitable habitats in Thailand inside (~82,400 km 2 ) and outside (95,000 km 2 ) PAs than Prayoon et al (2021), who predicted 39,508 km 2 of total suitable habitat.Our predictions used NDVI and land coverage fractions (Table S4) for predicting greenness, which may be useful for predicting the vegetation quality and availability for ungulates (Borowik et al., 2013).
However, NDVI is difficult to differentiate vegetation variations (Didan et al., 2015;Martinez & Labib, 2023), such as between specific agricultural areas, grassland, and dense forest canopy.This may include other vegetation types other than the species' habitat in suitable areas and estimated larger suitable areas predicted in non-forest areas and non-PAs to be identified in our study, compared to Prayoon et al's study.
Other studies suggest B. gaurus does use crop plantations or man-made grasslands, which may increase the suitable areas in our prediction, even if these are not their natural habitats and lead to conflict between human and gaur (Chaiyarat et al., 2021).
We found a high percentage of predicted suitable areas in Eastern Plains Landscape (ELP) and Chhaeb WS in Cambodia; the former supports the likely largest B. javanicus population globally (Gray et al., 2012).However, our results showed low habitat suitability in Sundaic Southeast Asia, with just 2% of the total suitable area in Indonesia (mainly in Alas Purwo NP, Java) and 2% of the total suitable area in Malaysia.B.
javanicus populations and habitats in Southeast Asian islands (Borneo, Kalimantan, Java, and Bali) are threatened due to hunting for horn and meat consumption and habitat loss (Dewi et al. 2020).In Thailand, we found high suitability similar to previous studies in Eastern (Menkham et al., 2019) and Western forest complexes (Jornburom et al., 2020), including reintroduction areas in Salak Pra WS (Chaiyarat et al., 2019) and where recent recolonization by natural population movement has occurred in Mae Wong NP (Phoonjampa et al., 2021).
Bubalus arnee has been domesticated and bred as livestock, making it hard to distinguish between the free-grazing domestic buffalo and B. arnee as domesticated animals may replace wild animals in suitable habitats and cause the high suitable area prediction outside PAs, especially in overlapping habitats (Zhang et al., 2020).We estimate the highest percentages of suitable areas at Kaziranga NP in India, currently with the largest population of wild water buffalo (Kaul et al., 2019).Grasslands and flood plain areas of Manas NP (500 km 2 ) and Kaziranga NP (>850 km 2 ) in India contain the most suitable habitat and are the main population strongholds for B. arnee (Choudhury, 2014).Moreover, we identified suitable areas for B. arnee in Huai Kha Kheang WS in Thailand, located close to streams and plains, but it is around 43% of Huai Kha Kheang WS.Thus, we recommend protecting these habitats to ensure the protection of wild cattle, such as active patrolling to reduce illegal intrusions, snare removal and habitat management based on their diet diversity (McShea et al., 2019).
We estimated suitable C. sumatraensis habitats in 12 countries, all consistent with IUCN polygons.However, our models predicted 293,883 km 2 of suitable area, substantially smaller than the IUCN polygons (around 3,500,000 km 2 ) (IUCN, 2021).We identified suitable C. sumatraensis and B. gaurus habitats in Dawna-Tanintharyi Landscape (DTL), Western Forest Complex, Kaeng Krachan and Kuiburi, which are transboundary parks connecting Myanmar and Thailand.These areas were identified as important conservation areas to build habitat connectivity for maintaining large felids, their prey, and other animal populations (Greenspan et al., 2021).However, our models predicted low suitability areas where C. sumatraensis has been recorded in Myanmar, such as Hkakaborazi NP and Hponkanrazi (Lwin et al., 2021;Rao et al., 2010).Fortyseven percent (185/388) of the data occurrences were collected from Thailand (Table S2), and only 13% were collected in Myanmar, so we suggest collecting more occurrence data from these PAs in other regions for validating and improving the model predictions.
We predicted N. griseus to have the lowest total suitable habitat area (around 170,000 km 2 ) among the five species, using its best model (MSDM-LA), with suitable areas restricted to China and Myanmar.This area is smaller than the IUCN polygon (Duckworth et al., 2008).We identified suitable habitat along the border of northern Myanmar (Hkakaborazi NP), China (Three Parallel Rivers of Yunnan PA), and part of the far-eastern Himalayan Landscapes (FHL), from which N. griseus as been observed in field surveys (Li et al., 2014).In Thailand, we predicted N. griseus may only be distributed in the northern part, with 13 PAs having medium to very high (40 -90%) suitable areas for the MSDM-SSA models (Supplementary Fig. S8) at higher elevations, close to a recent study that reported the species in 11 PAs (Jangtarwan et al., 2020).
Pintana & Lakamavichian (2013) report potential locations outside Qinling PAs, supporting our findings that more than 50% of suitable areas are found outside PAs.
We acknowledge sampling deficiencies across the regions.We had fewer occurrences in Vietnam, Laos, Myanmar and Indonesia compared to Thailand, from which a large number of our data points came (30,512 points in Thailand,3,152  griseus from which these species have been reported.This would likely be improved if more spatial data were available for these species.
In this study, we included all subspecies data points in our model ensembles as we aim to extrapolate and predict the entire range of species' habitat suitability, but this may increase uncertainty (Dormann, 2007).These five bovids have multiple subspecies, including 3 subspecies of B. gaurus (Duckworth et al., 2016), B. javanicus (Gardner et al., 2016) B. arnee (Kaul et al., 2019) and C. sumatraensis (Mori et al., 2019), and 2 subspecies of N. griseus (Duckworth et al., 2008).Subspecies may vary in niche, climate and biological interactions that could affect the model predictions.For example, C. s. thar live in Himalayan areas, and our models predicted the only suitable area was in Khaptad NP, whereas the Annapurna Conservation Areas, where the species has been detected (Aryal, 2009), were unsuitable.The low habitat suitability of our study in Borneo for B. javanicus could be because climatic and geographic conditions differ for B.j. lowi and those in mainland Asia, affecting model transferability across different regions (Zhu et al., 2021).Equally, Mori et al. (2019) suggest that N. griseus should be reclassified as a single species with Himalayan (N.bedfordi) and Burmese goral (N. evansi).Future analyses must consider these taxonomic reclassifications.However, we modeled species level habitat suitability, rather than the subspecies, as we assume that there is less likely to be habitat and environmental condition variation at the subspecies level for these bovids (Smith et al., 2019).
We found that using the MSDM OBR technique made the predicted ecological niche closer to a realised niche for species with more restricted ranges like B. javanicus, B. arnee and N. griseus, with higher performance TSS values compared to No MSDM models.We recommend restricting the accessible area for predicting B. arnee potential habitat to reduce overprediction caused by overlapping areas with domestic water buffalo.However, OBR can be sensitive to the distribution of occurrence data, because it keeps predicted suitable areas close to the occurrence locations.This may lead to the exclusion of potential suitable areas driven by a lack of occurrence data in those areas.
For example, the B. arnee No MSDM predicted potentially suitable habitat around the Sre Pok Wildlife Sanctuary in Cambodia where the species is distributed (Gray et al., 2012;IUCN, 2008), but after the MSDM, this potential habitat was excluded as we lack occurrence data in Cambodia.Although our study showed slightly different TSS values between two different accessible area extents, we encourage testing the different accessible areas as it affects the model results (Anderson & Raza, 2010).Moreover, model performance varied with accessible area sizes and spatial restrictions, emphasising the need for careful accessible area definition in ecological modelling (Barve et al., 2011).Still, we suggest further studies assess habitat quality and finescale validation to assess high suitability PAs and their capacity to maintain populations.

Conclusion
Our study provided an overview of the suitable remaining habitat for threatened bovid species at a regional scale using high-resolution environmental variables and species occurrence data from multiple observation methods.Our predictions showed that the suitable areas are small and fragmented for all five species, and more than 50% of

Fig. 1 .
Fig. 1.Species occurrence data (yellow circles), IUCN polygons (blue areas) and study areas (grey areas) used in model building for five wild bovid species.First, a common large 'accessible area' (A) was used for all species for model building, and then species-specific accessible areas (B-F) for individual species.
S3 (SSA) and Fig. S4 (LA), and the binary maps which were used for calculating the suitable area in Supplementary Fig. S5.The performance of spatially restricted ensembles was higher in comparison with the No MSDM models, as the TSS was improved for B. javanicus, N. griseus and B. arnee.The lowest performing model for B. arnee was the No MSDM-SSA (TSS = 0.57).The best model for B. gaurus was no MSDM-LA, B. javanicus and N. griseus is MSDM-LA, B. arnee is MSDM-SSA, and C.

Fig. 3 .
Fig. 3. Habitat suitability prediction maps for five species (A-E) using the best model from either accessible areas and the weighted average ensemble.The value ranges from 0-1: yellow represents low suitability and dark brown represents high suitability.(F)

Fig. 5
Fig.5Habitat suitability for five species in Thailand.The value ranges from 0-1: yellow represents low suitability and dark brown represents high suitability.(F) presents the map of species numbers for five species using the binary suitability results.
suitable areas are outside of protected areas.Those suitable areas outside PAs could possibly become efficient conservation areas, such as forest corridors or buffer zones to connect fragmented bovid populations and enhance long-term habitat conservation.Our predictions may inform conservation actions to avoid further defaunation of wild bovidae such as the management of human-wildlife conflicts and habitat quality for long-term species survival.Caesar Kleberg Wildlife Research Institute, Texas A&M University, USA), Freeland Foundation, WCS Thailand, Smithsonian Institution, PANTHERA USA in Thailand.WWF Thailand would like to thank: WWF Germany, WWF Sweden, WWF US, B.Grimm, WWF Japan, and WWF Switzerland for wonderful support of field projects, and Department of National Parks, Wildlife, and Plant Conservation for kind permission and collaboration.WCS Cambodia, Friends of Wildlife, Wildlife Alliance and Ministry of Environment, Royal Government of Cambodia for the camera trap data in Cambodia.CarBi Project of WWF Lao and WWF Greater Mekong for the camera trap data in Laos.Camera trap in Myanmar: Friends of Wildlife and Data collecting in Tanintharyi, Myanmar was partially funded by the European Union, Helmsley Charitable Trust and mainly funded by Integrated Tiger Habitat Conservation Project, through the Fauna & Flora International (FFI).Open source databases: https://www.gbif.org,camera trap data from eMammals (https://emammal.si.edu/) by William J. McShea (Conservation Ecology Center, Smithsonian Conservation Biology Institute) and Megan Baker-Whatton (Smithsonian Conservation Biology Institute, and George Mason University) under the project: Qionglai Mountains Project, Liangshan Mountains Project, Habitat Connectivity in Minshan Mountains, Qinling Project, HKK ForestGEO Project, and Carnivore Intraguild Interactions in Select Thailand Reserves.Also, the other researchers who share the species occurrence data with us.We thank the IUCN specialists for commenting on the early version of the habitat suitability maps.The authors wish to acknowledge the use of New Zealand eScience Infrastructure (NeSI: https://www.nesi.org.nz)high-performance computing facilities, consulting support and/or training services as part of this research.
For C. sumatraensis, our findings showed 41% (110,535 km 2 ) of the predicted suitable area inside PAs by the best model (No MSDM-SSA).A high percentage of suitable area (≥80%) was identified within 64 PAs (from 835 PAs, Fig.4).Most suitable PAs are found in Southeast Asia, including Thailand, Cambodia, Myanmar, Laos and Malaysia.
suitable PAs were identified in Sri Lanka (128 PAs from 146 PAs, Fig.4).India had the highest proportion of suitable area in PAs (≥80%), identified in Kaziranga NP and Manas WS.
Kaengkrachan and Huai Kha Khaeng, for C. sumatraensis Ton Pariwat, Khlong Yan and Khlong Naka, for B. arnee Phu Wua and Dong Yai, and for N. griseus Mae Wa, Doi points outside Thailand, Supplementary TableS2).Occurrence data based on data accessibility can have sampling bias, particularly clustered points for B. gaurus, B. javanicus, and C. sumatraensis.We minimised these biases through spatial thinning.Since we found large amounts of the suitable area outside of Thailand, the imbalance has not likely impacted our results, however, we suggest that future studies should focus on ensuring monitoring of bovid populations in other countries, especially in India and Myanmar.However, missing data has impacted some results.The model TSS values for endangered B. javanicus and N. griseus are over 0.8, yet our models predict unsuitable areas in part of Indonesia (east and central Kalimantan; Dewi et al. 2020) for B. javanicus and China (e.g.Beijing and northeast Inner Mongolia; Yang et al. 2019) for N.