MODELING THE DISTRIBUTION OF THE SOUTHERN YELLOW-CHEEKED GIBBON (NOMASCUS GABRIELLAE) USING MAXENT

The Southern yellow-cheeked gibbon (Nomascus gabriellae) is an endangered species found only in a small region of Indochina, and its populations have declined in most known sites. In this study, we use Maxent, a robust and widely used species distribution modeling approach, to predict the current and future distributions of the Southern yellow-cheeked gibbon over its entire range based on an extensive review of published records. Overall, we compile and provide a comprehensive set of known distribution records of the species from Viet Nam and Cambodia. The model results show that N. gabriellae is likely to occur in many areas around the Central Highlands in both Viet Nam and Cambodia and the southern end of Laos. Our study suggests that protected areas in this region will play a key role in conserving the gibbons. In addition, the distribution of the gibbon under future climate conditions, even in the best case scenario, is likely to shrink significantly, as the species would have to move to higher altitudes. Under such an impact, the populations would become more fragmented and restricted to a few areas of higher elevation. Our study also confirms that the climatic difference across its distribution range may not be entirely responsible for the speciation and biogeography of the N. annamensis/gabriellae complex.


INTRODUCTION
Nomascus is the second most speciose genus of gibbons with a total of seven species, of which six are found in Viet Nam [1]. All six gibbon species in Viet Nam have been severely threatened due to illegal hunting and habitat destruction, and are listed as Endangered or Critically Endangered in the IUCN Red List [1]. Unlike some other Nomascus species that have only been discovered or rediscovered recently, the Southern yellow-cheeked gibbon (N. gabriellae) has an ostensibly longer history of research [1 -3]. N. gabriellae was formerly known to be widely distributed throughout eastern Cambodia, southern Viet Nam and southern Laos, but its status in Laos is controversial [3].
Although there were no specific population data for this species before the 1990s, given the scale of illegal hunting and habitat destruction in Viet Nam over the past 40 years [4], it is widely believed that the population of Nomascus gabriellae has declined in most locations [1], resulting in being listed as endangered on the IUCN Red List [5]. Taxonomic studies of N. gabriellae have proposed a new species, N. annamensis, and that has led to a substantial decrease in the distribution range of N. gabriellae sensulato [3,6]. N. gabriellae and N. annamensis have been hypothesized to be separated by barriers formed by the Ba River in Viet Nam and the Srepok River in Cambodia, with the Southern yellow-cheeked gibbon inhabiting the southern side of the river system [7]. However, N. gabriellae has been reported to occur in the northern side of the Ba River, hence the rivers may not be a tight boundary for the distribution of these two gibbon species [6,8].
Species distribution modeling (SDM) is a relatively new approach that has proven useful in studying biogeography, biodiversity pattern, evolutionary ecology, and environmental issues [9]. In general, SDM explores and uses possible interactions between species occurrence records and defined environmental variables to calculate and generate a probability map for a species in a given area. SDM can help identify the potential distribution range of an endangered species that could be included in protected area planning [10], resolve taxonomic issues for species groups whose other traits may not be clearly distinguishable and conclusive [11], predict areas of possible spread and vulnerability due to risk from invasive species [12], find stable areas of refugia that can be home to many rare and endemic species [13], especially in the context of climate change [14], and examine the processes that affect biodiversity and evolution over a long period of time [15]. The wide range of applications has spurred the development of many different SDM approaches, of which Maximum Entropy (Maxent) is one of the most commonly used [9]. Maxent is a relatively powerful machine learning method, and unlike other approaches that need both presence and absence data, Maxent only requires presence records [16][17][18]. Maxent also allows to predict the distribution of species in future climate scenarios to assess the impact of this threat [16,17].
A previous study (Thinh et al., 2018) attempted to model the current distribution and impacts of climate change on the Southern yellow-cheeked gibbon [18]. However, the authors used the older version of bioclimatic variables from WorldClim version 1, as well as global circulation models of future climate change as guided by the Coupled Model Intercomparison Project 5 (CMIP5). As newer models have been developed and recently released for CMIP6, with more recent data and improved accuracy [20], it is necessary to re-assess the distribution of Nomascus gabriellae under both current and future climate conditions. In this paper, we reviewed and compiled a profile of the Southern yellow-cheeked gibbon across its known distribution range. We then incorporated the local data into the SDM to create a distribution map of the species using Maxent and new climate change models to help improve understanding and conservation measures for this primate.

Literature review
We reviewed the available records of the Southern yellow-cheeked gibbon by searching the Web of Science, Google Scholar, Research Gate, GBIF, and iNaturalist using the following queries: "Nomascus gabriellae", "Nomascus", "gabriellae", "Red-cheeked Gibbon", "Yellowcheeked Gibbon", "Buff-cheeked Gibbon", "Crested Gibbon", and their equivalent keywords in Vietnamese such as "Vượn má vàng", "Vượn đen má vàng", "Vượn đen má hung", "Vượn". In addition, library archives, reports, and specimens from related institutions were also examined. For all records mentioning the Southern yellow-cheeked gibbon but coming from the distribution range of the recently described species N. annamensis, we still reported the localities but excluded them from modeling. The collected records were then evaluated, checked, and filtered to avoid erroneous or unreliable locality information, and then the final set was used to train the SDM for the Southern yellow-cheeked gibbon.

Data preprocessing and model tuning
From the collected records, to avoid spatial autocorrelation, we used the spThin package [21] in R [22] to thinly split the localities over a distance of 10 km [23], resulting in the final set of 98 localities. We constructed the SDM using 19 bioclimatic variables at 2.5-minute resolution available at WorldClim 2.1 database [20], and limited the extent using one-degree buffer around minimal convex polygon around the occurrence localities [24]. We ran all analyses in Maxent version 3.4.1 [17]. However, since Maxent tends to produce overfitting models, we performed the following set of tuning steps to minimize overfitting and maximize discriminatory ability using ENMeval package in R [25].
We performed the tuning process using all combinations of feature classes and tested the models with regularization multipliers ranging from 1.0 to 10.0 in increments of 0.5. Other model parameters, e.g. convergence threshold and feature selection, followed recommendations from model developers [17]. We then used the spatial fourfold cross-validation method recommended for models that need to be projected under novel climatic conditions [25] to build the model. This method randomly splits the occurrence data into four equally sized partitioned folds based on the spatial distribution of the input localities. Four models were then created, leaving out one-fold each time as test data for model evaluation.

Model selection and future projections
To evaluate model performance and select the best fit, we used a 10 % omission rate threshold to select models that showed the least overfitting. In this subset, we then chose the model with the highest AUC value. The final models were then compared using the Akaike information criterion (AIC), which balances complexity with model fit [26]. For the final model, we used a 10 % training presence threshold to classify between suitable and unsuitable areas for the Southern yellow-cheeked gibbon.
The best model was then projected to two future periods, 2041 -2060 and 2061 -2080, using data for two global climate models, BCC-CSM2-MR [27] and MIROC6 [28]. The two models have been used in similar studies with comparable geographic locations, and are only distantly related to each other, hence they may better reveal uncertainties in the Maxent projections [29 -31]. All four scenarios, including SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5, were used for each model, ranging from SSP1-2.6, which assumes there will be an increasing shift towards sustainable practices in the global economy, to SSP5-8.5, which assumes the world will adjust towards an energy-intensive, fossil fuel-based economy [32]. The individual binary maps for each scenario of the same model and of the same duration were then stacked and combined to generate presence probability maps for the species.

Distribution records review for the Southern yellow-cheeked gibbon
Our review of available literatures, including published papers, reports, and unpublished field records yielded a comprehensive set of locality records of the Southern yellow-cheeked gibbon from Cambodia, Laos, and Viet Nam. We have explicitly specified every possible record of Nomascus annamensis. All records of the species are illustrated in Supplementary material and Figure 1.

Species distribution modelling results for current climate conditions
For the SDM, the Maxent models showed reasonable predictive ability for the distribution of the Southern yellow-cheeked gibbon with mean AUC values greater than 0.82. The best model had a regularization multiplier value of 6.0 and a combination of linear, quadratic, and hinge feature classes, and an AUC value of 0.839. All the final models were quite similar in their ability to predict the overall distribution of Nomascus gabbrielae and differed only slightly in exact locations and total suitable areas.
However, a regularization multiplier value of 6.0 for the best model means that the final model was quite flexible, which can lead to a broader and larger distribution. The final prediction should therefore be carefully interpreted as "potential zones", not "core zones". Hence, it may include regions that are not highly likely to be suitable for the Southern yellowcheeked gibbon.
The current potential distribution range of the Southern yellow-cheeked gibbon can be broadly divided into two main regions, and the southern part contains the largest area and covers most of its currently recognized range. In particular, the suitable area covers parts of Kon Tum, Gia Lai, Dak Lak, Dak Nong, Lam Dong, Khanh Hoa, Ninh Thuan, Binh Thuan, Binh Phuoc, Binh Duong, Dong Nai, and Ba Ria -Vung Tau in Viet Nam as well as Kratie, Mondulkiri, and Ratanakiri in Cambodia. The total suitable area of the gibbons is estimated at about 77,000 km 2 . Compared to a previous study [19], our model predicted a larger suitable area, and also included more locations in the fringe regions of the map, such as Gia Lai and Ninh Thuan in Viet Nam, where the gibbons were previously documented in historical records.

Species distribution modelling results for future climate scenarios
Based on future model projections, it is clear that Southern yellow-cheeked gibbon's distribution would shift in future climate change scenarios. Changes in the distribution are illustrated in Figure 3, which is a probability density map, ranked on a graduated color scale, showing the likelihood of the gibbon presence calculated from a large ensemble [33]. Each ensemble is a combination of all four scenarios for one particular timeframe from one global circulation model. The map shows the highest value (100 %) which means that the area was estimated to be suitable for the gibbons in all scenarios. Although the predictions of those models differ slightly in terms of the overall location of the changes, there appear two trends for Nomascus gabriellae's range shift in the future. Firstly, the distribution is likely to shrink significantly compared to the current predicted range, with an average decrease of approximately 41.3 % (at standard deviation ± 12.8 %) in terms of total area, even in the best-case scenario ( Table 1). The second trend shows that the species would move upwards to higher elevations, 500 m or more, in the Central Highlands. As shown in Figure 3, the northern part of the distribution range will concentrate around the Kon Tum plateau, and the southern part will be situated around the Dak Lak, Da Lat, and Di Linh plateaus. The populations in lower areas, such as parts of Dong Nai -Binh Phuoc in Viet Nam or Mondulkiri -Kratie region in Cambodia will become restricted and occupy a few higher elevation areas. In such case, the landscape of Dak Lak -Da Lat -Di Linh plateaus may serve as the main refuges for the species in future. Furthermore, while two global circulation models used in our research are only distantly related, they showed quite similar results in terms of future suitable areas and suitable area changes in the same timeframe. Also, for different timeframes, both global models suggested that the range contractions of the gibbons in the 2061 -2080 period, approximately 38,000 km 2 , are likely to be larger than those in the 2041 -2060 period, which means that about 27,000 km 2 will be lost. Therefore, in the long term, climate change will pose an ever-increasing threat to the Southern yellow-cheeked gibbon.

Discussion
Based on these results, we suggest that protected areas (sanctuaries) located in continuous and climate suitable habitats for the gibbons should be prioritized for future species conservation measures. They are presented in Figure 2 and include: A Yun Pa Nature Reserve (Viet Nam): The northern limit of Nomascus gabriellae is at A Yun Pa Nature Reserve. Although the total size of the gibbon population in this area has been unknown, recent records indicate that this area still harbors an important gibbon population [8].
Chu Yang Sin National Park, Bidoup -Nui Ba National Park, and Phuoc Binh National Park landscape (Viet Nam): The adjacent protected areas form a relatively intact and continuous forest region of about 200,000 ha. As a large population of the Southern yellow-cheeked gibbon has been recorded in this landscape [1,29,30], it should be a priority for any species conservation initiative. Nam Nung and Ta Dung Nature Reserves (Viet Nam): The small size of both reserves, and their relative isolation from other sites, could be problematic for long-term conservation measures. However, the sites are known to support a significant population of the Southern yellow-cheeked gibbon, and are still among suitable habitats [31,32].
Bu Gia Map National Park (Viet Nam), and Phnom Prich -Seima -Snoul Wildlife Sanctuaries (Cambodia) landscape: Together they form a vast network of transboundary protected areas, probably home to the largest population of N. gabriellae [1,33]. They will undoubtedly play a key role in any global conservation initiative for the species.
Cat Tien National Park, and Vinh Cuu Nature Reserve (Viet Nam): Cat Tien is considered a key conservation site for the Southern yellow-cheeked gibbon with the largest population of the gibbon [1,34]. Therefore, Vinh Cuu and surrounding areas should be prioritized for the conservation of N. gabriellae in Viet Nam.
In this paper, we used bioclimatic variables with 2.5-minute resolution, which are the finest data available to date for future climate scenarios from WorldClim version 2. While they are indeed coarse scale variables, and 30-sec resolution data may be available soon, we argue that SDM, even at a coarse scale, does help approximate the locations for species range, and where we should focus our analysis with fine-scale data in the future. But indeed, future studies should perform modelling with finer-scale analyses, as higher resolution data may improve the results shown here for more targeted conservation initiatives. Including other environmental variables in future analyses will undoubtedly enhance the model performance as well.
Based on the model results, the current potential distribution range of the Southern yellowcheeked gibbon can be broadly divided into two main regions, and the southern part contains the largest area and covers most of its currently reported range. Two regions are connected by an eastern small "corridor", which is coincidentally located in the A Yun Pa Nature Reserve. Interestingly, the whole northern part falls totally outside the current reported range of the N. gabriellae, and well within the southern distribution limit of the N. annamensis, a subspecies of the N. gabriellae until recently [3,6]. The "corridor", A Yun Pa, has been considered as the northern limit of the Southern yellow-cheeked gibbon, and even the reserve itself is located in the northern side of the Ba river, a purported barrier for two gibbon species [6,8]. Like previous studies [6 -8], which have raised concerns regarding the unclear natural barrier limit, our study reveals that environmental variables alone may not be responsible for the separation between the two species of the N. annamensis/gabriellae complex. Moreover, Rawson et al. (2011) noted that in Cambodia, the biogeographic barrier between N. gabriellae/annamensis could be the dry deciduous dipterocarp forests, located on the eastern side of the Mekong River, extending from Banlung to around Phnom Prich [1]. The lack of suitable food sources in the region is perhaps a bigger problem for species dispersion, especially for a highly frugivorous and folivorous taxa such as Nomascus [40].

CONCLUSIONS
In this study, the final Maxent model showed reasonable predictive power with mean AUC values greater than 0.82, and the best model had an AUC value of 0.839. All final SDMs were similar in terms of overall distribution pattern of Nomascus gabriellae. The results showed that N. gabriellae inhabits areas along Dak Lak -Da Lat -Di Linh plateaus in both Viet Nam and Cambodia sides, and protected areas in this region, including A Yun Pa, Chu Yang Sin, Bidoup -Nui Ba, Phuoc Binh, Nam Nung, Ta Dung, Bu Gia Map, Cat Tien, and Vinh Cuu in Viet Nam, as well as Phnom Prich, Seima, and Snoul in Cambodia, will play a key role in any conservation action for the gibbons, especially in the context of climate change where the species is thought to be restricted to higher elevations. The climate change will also pose a serious threat to the species in the long term, as the gibbons' distribution is likely to shrink significantly with an average decrease of approximately 41.3%. Our study suggests that environmental variables alone can't explain the current distribution of N. gabriellae and N. annamensis. Other biological factors, such as competition and isolation by distance or both, should be examined in future studies to determine the role of evolutionary processes in separating the two sister species.
Credit authorship contribution statement. CTHN: Data analysis, model running, and paper writing. LDM: Theoretical development and paper writing. NTA: Theoretical development, data analysis, model running, and paper writing.

Declaration of competing interest.
We have no competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.