Climatic niche modelling and genetic analyses highlight conservation priorities for the Spotted Softshell Turtle ( Pelodiscus variegatus )

The Spotted Softshell Turtle ( Pelodiscus variegatus ) has been recognised since 2019 from Vietnam and Hainan Island, China, but little information about its population status and distribution range is currently available. The species has been provisionally listed as Critically Endangered by the Turtle and Tortoise Working Group, although the status has not been officially accepted by the IUCN, due to the threats the species is facing, including habitat loss and degradation, overexploitation for food, competition with other non-native softshell turtles and pollution. To identify conservation priority sites for P. variegatus in mainland Indochina, this study combines molecular analyses and species distribution modelling. Our results show that, in Vietnam, Phong Nha-Ke Bang National Park has the largest suitable area and high probability of species occurrence, followed by Vu Quang National Park and Song Thanh and Ke Go Nature Reserves. In addition, the central provinces, from Thanh Hoa to Thua Thien Hue in Vietnam, constitute a key part of the species distribution and should be prioritised for conservation actions. According to the study’s findings, although P. variegatus is possibly found in Laos, the probability decreases sharply at the border between both countries and there is also a gap in the occurrence of wetlands, arguing for strong natural barriers. Unfortunately, to date, only part of the species potential distribution is protected, while no records are known from protected areas, highlighting the need for extended or even new reserves. To recover natural populations of the species and following the IUCN’s One Plan Approach to Conservation, breeding programmes have been established in Vietnam with a potential to expand to other facilities in the country and abroad. Once suitable sites are identified, offspring can be released into the protected ar - eas to improve the current conservation status of this highly-threatened softshell turtle.


Introduction
The Chinese Softshell Turtle, Pelodiscus sinensis, was described nearly 200 years ago and was believed to represent a morphologically variable, geographically widespread taxon (from the Russian Far East through the Korean Peninsula, eastern and central China to Vietnam).The Northern Chinese Softshell Turtle, P. maackii (Brandt, 1857), was described 23 years later, but was then thought to be a synonym of P. sinensis.Only 35 years ago, the populations from the northernmost part of the P. sinensis distribution range were shown to represent a distinct species, based on osteological differences (Chkhikvadze 1987).
Two additional species of this complex from central China were described in the 90s, based on morphological differences: the Hunan Softshell Turtle (P.axenaria) and the Lesser Chinese Softshell Turtle (P.parviformis) (Zhou et al. 1991;Tang 1997).Recent molecular studies confirmed that the genus Pelodiscus constitutes a species complex (Fritz et al. 2010, see also Stuckas and Fritz (2011); Yang et al. (2011); Gong et al. (2018)).Based on genetic and morphological analyses, three new taxa were described: the Spotted Softshell Turtle, P. variegatus from northern Vietnam and Hainan (China), the Horse-hoof Softshell Turtle, P. huangshanensis from southern Anhui Province of China and the Chinese Stone Slap Softshell Turtle, P. shipian from Jiangxi and Hunan Provinces of China (Farkas et al. 2019;Gong et al. 2021;Gong et al. 2022).
Thus, the Pelodiscus sinensis complex at the moment comprises seven species, with six species being distributed in China and four endemic to the country (Gong et al. 2021).The revisions have consequences on the species conservation, as taxonomic splitting implies that the range size and number of individuals decrease for each species.However, these new research results have yet to be reflected in the conservation assessment of the genus.According to the IUCN (2023), only P. sinensis is listed as Vulnerable, but the data are greatly outdated with the last update in 2000 (Asian Turtle Trade Working Group 2000).The Turtle Taxonomy Working Group (TTWG) recently evaluated conservation status of six species in the genus and provisionally categorised three species as Data Deficient (P.axenaria, P. huangshanensis and P. maackii), one as Vulnerable (P.sinensis) and two as Critically Endangered (P.parviformis and P. variegatus; TTWG 2021).The last species, P. shipian, has not been assessed by any previous study.
The exact range of the recently-described P. variegatus is still largely unknown, but historical records suggest that the species occupies lowland areas of central and northern Vietnam and parts of southern China, viz. Hainan Province (Farkas et al. 2019;TTWG 2021).Recently, Ziegler et al. (2020) investigated whether natural populations of P. variegatus still exist, since the description of this species was mostly based on historical museum specimens.To find potential members of P. variegatus, this study focused on surveys of central lowland freshwater habitats, as well as local markets, restaurants and farms in Vietnam.Individuals with the species-specific dark blotched plastron pattern were identified and subsequently genetically screened to confirm their identity.By using this approach, Ziegler et al. (2020) demonstrated that P. variegatus is still extant in the wild in Vietnam, both based on evidence from the trade and on surveys in the natural habitat.The study also showed that P. variegatus occurs in the central provinces of Thanh Hoa, Nghe An, Ha Tinh and Quang Binh, primarily in lakes with flat shores consisting of mud and soft soil, rivers in agricultural landscape and medium-sized streams in secondary forests.
As softshell turtles are common and prized as food where they occur, natural populations are threatened by local hunting, with further threats of habitat loss and competition with introduced softshell turtles (Shi et al. 2008;Le Duc et al. 2020).At the moment, there is little evidence of any viable population of the P. variegatus existing in the wild.It is, therefore, imperative to implement conservation measures in priority areas, where natural populations still likely survive.To establish conservation priorities for the species in Vietnam, the present study performed climatic niche modelling, based on comprehensive distribution data derived from detailed molecular analyses of existing and new samples collected from the country, the Emys system (emys.geo.orst.edu), the Turtle Taxonomy Working Group (TTWG 2021) and the Global Biodiversity Information Facility (GBIF).

Molecular analyses
To identify new samples collected in the field, we used a molecular approach.In total, 61 newly-collected samples were included in the analyses (Suppl.material 1).Sequences of other Pelodiscus species were downloaded from Gen-Bank.Two species, Dogania subplana and Palea steindachneri, were employed as outgroups (Le et al. 2014).We used the protocols of Le et al. (2006) for DNA extraction and amplification.Two mitochondrial genes, the nearly complete cytochrome b (1110 bp) and partial ND4 (673 bp), were amplified using primers listed in Table 1.Successful amplifications were purified to eliminate PCR components using GeneJET™ PCR Purification Kit, Thermo Fisher Scientific (Vilnius, Lithuania).Purified PCR products were sent to 1 st BASE (Selangor, Malaysia) for sequencing.Sequences generated in this study were aligned using De Novo Assemble function in the program Geneious v.7.1.8(Kearse et al. 2012).
Data were then analysed using Maximum Likelihood (ML) as implemented in IQ-TREE v.1.6.7.1 (Nguyen et al. 2015) and Bayesian Inference analysis (BI), as implemented in MrBayes v.3.2.7 (Ronquist et al. 2012).For ML analysis, we used a single model and 10,000 ultrafast bootstrap replications.The optimal model for nucleotide evolution employed in both methods was determined using jModelTest v.2.1.4(Darriba et al. 2012).The optimal model for nucleotide evolution was set to TIM2+G for ML and single-model Bayesian analyses.For the Bayesian Inference, two independent analyses with four Markov chains (one cold and three heated) were run simultaneously for 10 million generations with a random starting tree and sampled every 1000 generations.Log-likelihood scores of sample points were plotted against generation time to determine stationarity of Markov chains.Trees generated before log-likelihood scores reached stationarity were discarded from the final analyses using the burn-in function.The posterior probability values for all clades in the final majority rule consensus tree were provided.The cut-off point for the burn-in function was set to 71 in the Bayesian analysis, as -lnL scores reached stationarity after 71,000 generations in both runs.Nodal support was also evaluated using ultrafast bootstrap (UFB) in IQ-TREE and posterior probabilities (PP) in MrBayes.PP and UBP ≥ 95% were regarded as strong support for a clade (Ronquist et al. 2012;Nguyen et al. 2015).This study only employed mitochondrial genes to provide taxonomic identification of samples collected from the wild.Although maternally inherited mitochondrial loci cannot help detect hybridisation events, interbreeding between different species of softshell turtles has only been reported in turtle farms (Gong et al. 2018).In addition to genetic screening of captured animals, we morphologically identified the specimens using diagnostic characters.

Species records and environmental variables
For species records, we carefully checked potential records listed in recent studies, including, Le Duc et al. (2020), Pham et al. (2020), Ducottend et al. (2023), Pham et al. (2023a) and Pham et al. (2023b).However, the papers followed the older version of the Turtle of the World Checklist (TTWG 2018), which omitted information about P. variegatus because the species was only described a year later in 2019.As a result, the records were all assigned to P. sinensis, although the species is now considered only occurring in China and Taiwan according to the new checklist (TTWG 2021).In addition, much of the information came from interviews with local people and is not taxonomically robust.In Pham et al. (2023b), some photos clearly belong to the P. sinensis form, while the others did not show diagnostic characters to allow identification of these individuals.We, therefore, did not add the data to our analyses.In the final dataset, we used 54 unique, georeferenced and genetically confirmed locations located in unique grid cells, as shown in Fig. 2

(see below).
For environmental predictors, we used a combination of weather station-derived precipitation data and remote sensing data.The full list of variables including their interpretation is provided in Table 2. Average annual characteristics of precipitation regimes were obtained from the Worldclim database 2.1 as interpolated elements from different climate conditions collected over a period of 30 years  with a resolution of 30 arc seconds (Hijmans et al. 2005;Fick and Hijmans 2017).In order to characterise seasonal changes in land surface temperatures and in vegetation cover, we used 27 remote sensing predictors derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors of two NASA satellites available through the EDENext project (temporal resolutions: 8-day averages (MOD11A2) and 16-day averages (MCD43B4), spatial resolution 30 arc seconds) (Mu et al. 2007).Availability of wetlands as surrogate of suitable microhabitats was assessed using a recent assessment of tropical wetlands provided by Gumbricht et al. (2017) as a categorical variable.As the wetlands dataset lacks river networks which may provide suitable microhabitats for the turtle, we added a high-resolution water layer as additional category (GRDC 2020).The wetlands dataset was resampled to the spatial resolution of the remote sensing predictors and only grid cells including water bodies were considered.
Multi-collinearity amongst predictors may hamper successful model training and subsequent projection (Brun et al. 2020).Therefore, we computed pairwise Spearman rank correlations and selected each predictor pair with ρ 2 > 0.75 with only one predictor for final model training.The final set of variables comprised both climate and land-cover related variables, which are suitable to characterise the species range limits and microhabitat preferences.As these operate on very different scales, we used a two-step approach for modelling: (1) determining likely range limits with temperature and precipitation related variables and (2) using these range limits to train a second model using land-cover related variables to assess microhabitat preferences.For the first climate related model, the variables comprised the annual mean surface temperature (ED15078_ bio1), maximum surface temperature of the warmest month (ED15078_bio6), surface temperature annual range (ED15078_bio7), annual mean precipitation (bio_12), precipitation of the driest month (bio_14) and precipitation of the warmest quarter (bio_18).Quantification of microhabitat preferences was derived from the annual mean Normalised Difference Vegetation Index (NDVI; ED1514_bio1), annual range of NDVI (ED1514_bio7), annual mean Enhanced Vegetation Index (EVI; ED1515_bio1), annual range of EVI (ED1515_bio7) and the categorical map of wetlands including the categories open water, mangroves, swamps, flood-out swamp, fens, riverine, floodplain, marshes, marshes-dryland/wetland and marshes-wet meadows.

Species distribution modelling
As the algorithm for climatic niche modelling development, we used Maxent v.3.4.0, which is specifically designed to derive potential distributions from presence-pseudoabsence data (Phillips et al. 2006;Phillips et al. 2017) and which can perform well even if the sample size is comparatively small (Hernandez et al. 2006).For the first model, we chose an area defined by a minimum convex polygon buffered with 5 km enclosing all species records as environmental background.Model selection followed the procedure described in Ginal et al. (2022).In Maxent, we allowed only linear, product and quadratic feature classes and used a regularisation multiplier of 0.8, as theses settings had, on average, the minimum delta AICc (689.4 with 8 parameters) and revealed the most realistic response curves.Using these optimal settings for feature classes and regularisation parameter and a bootstrap approach, we computed 100 models, each trained with 80% of the species records and using the remaining 20% for model evaluation via the area under the ROC (Receiver Operating Characteristic) curve [AUC] (Swets 1988).The average across all 100 replicates in cloglog format was used for further processing.
For the second land-cover niche modelling, we reclassified the potential distribution suggested by the first model applying the minimum training presence threshold, which was used as environmental background.Model selection followed again Ginal et al. (2022) and internal settings in Maxent were set to a regularisation parameter of 0.7 and linear, product and quadratic feature classes.The regularisation parameter of the categorical wetland predictor was set to 0.250.The average delta AICc was 760.6 with 10 parameters.Again, the average across all 100 replicates in cloglog format was used for further processing.The joint effects of climatic suitability and microhabitat suitability were estimated by rescaling both average predictions to a scale of 0-1 after applying the minimum training presence threshold and multiplying both.The resulting map highlights areas where both climatic and microhabitat suitability are highest.

Defining conservation priority sites
We merged the occurrence data with existing protected areas (Reserves, National Parks etc.) in the country.Information of protected areas was obtained from the world dictionary of protected areas/protected planet (https://www.protectedplanet.net).We selected the targeted area in north-central and central Vietnam because most of distribution records were reported from the region.In total, there were 42 potential conservation units within the general area and, for each Reserve, we computed the number of suitable grid cells, the sum of probabilities and the mean probability in QGIS 3.14.Map resolution ca. 1 km (30 arc sec).

Results
The molecular matrix contained 1921 aligned characters.Both BI and ML analyses showed that new samples belong to Pelodiscus sinensis and P. variegatus.The former species was only moderately supported (PP = 89%, UBF = 90%), while the latter received strong support from both BI and ML (PP = 100%, UBF = 99%).In total, 61 newly-collected and four GenBank samples were identified as P. variegatus (see Suppl.material 1: fig.S1, for the full tree).The localities of the samples were used for the species distribution modelling.
When integrating the probabilities of occurrence derived from climatic and land-cover variables, the most suitable habitats for P. variegatus are near the Vietnam coastline, where extensive freshwater wetlands exist (Fig. 2C).Rivers and other water bodies in mountainous areas are partly suitable, but they cover much less area.
Only part of the potential distribution for P. variegatus in Vietnam is protected and no known occurrence is directly located within reserves, although they are close, such as nearby Ke Go Nature Reserve and Vu Quang National Park (Fig. 2D).Table 3 has information on the relative ranking of protected areas, based on their IUCN status, the total area suitable for P. variegatus, the sum of potentially suitable sites and the average potential for the species.The Table is sorted with descending total suitable areas per Reserve.

Discussion
Farkas et al. ( 2019) stated that, in Vietnam, most records of P. variegatus fall within the "Northeast Lowlands Subregion" of Bain and Hurley (2011).The zoogeographical affinities of Hainan are closely related to this area as well as mainland south-western China, while the southern portion of the purported range forms part of the "Central-South Vietnam Lowlands Subregion" as defined by Bain and Hurley (2011).However, our molecular data show that most distribution localities of the species occur in the north-central region of the country, except for records from Dong Mo Lake (Fig. 1).It is likely that a majority of extant populations of the species is restricted to this region in Vietnam, although it is possibly found in Laos, based on climatic niche modelling results.However, the probability based on climate decreases sharply at the border between both countries and there is also a gap in the occurrence of wetlands.Thus, the border area may represent a strong natural barrier.
Our habitat suitability analysis predicts that two most important protected areas for P. variegatus include Phong Nha-Ke Bang and Vu Quang National Parks in Vietnam.Other protected sites with the largest suitable sizes consist of Song Thanh and Ke Go Nature Reserves and Ben En National Park (Table 3, Fig. 2).In terms of Average Probability, Hue Saola Nature Reserve receives the highest value (0.422), followed by Phong Dien Nature Reserve, Bach Ma National Park and Ke Go Nature Reserve (Table 3).The findings show that the central provinces from Thanh Hoa to Thua Thien Hue form an important part of the species distribution.Although one record in our study suggests that the species occurs in Laos (Fig. 2), its approximate field coordinates could not be used to absolutely confirm the species' presence in the country.It is, therefore, essential to conduct field surveys at suitable sites in Laos to verify its presence or absence.
Prior to its discovery, P. variegatus was considered part of P. parviformis.The latter species was already assessed as threatened and included on appendix II of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES).As its southern population has been split into a different species,  P. variegatus, the overall population of each species becomes even smaller than previously thought.Rhodin et al. (2018) regarded P. parviformis as "Critically Endangered (CR)".Consequently, TTWG (2021) provisionally categorised P. variegatus as CR, although the assessment has not been officially accepted by the IUCN Red List.The species was recently proposed to list as Vulnerable (VU) in the Viet-nam Red Data Book, based on estimation of population reduction approximately of over 30% in the past of 30 years (Nguyen, per. comm. 2023).
Ziegler et al. ( 2020) found during their market and trade surveys in Vietnam that human of softshell turtles appears massive.In addition, interbreeding events were observed to occur amongst softshell turtle species in farms.Most of the inhabited freshwater bodies and their surroundings showed signs of human encroachment, such as fishing, vegetation transformation, conversion and pollution.Thus, both in situ and ex situ conservation measures seem essential for protecting P. variegatus from extinction.Some of the investigated freshwater habitats are already located inside protected areas, such as Ben En, Phong Nha-Ke Bang and Vu Quang National Parks and Ke Go Nature Reserve in central Vietnam.Some of the protected areas are well known internationally as either special bird areas (Ke Go Nature Reserve) or with spectacular mammals, including Saola (Pseudoryx nghetinhensis) and the Large-antlered Muntjac (Muntiacus vuquangensis) (Vu Quang National Park).Improving conservation in and around the protected areas will benefit a suite of critically-endangered and endemic species, including the P. variegatus.As a first measure, based on the genetically-identified individuals, a conservation breeding programme has been established.This is following IUCN's One Plan Approach to Conservation, developed by the Conservation Planning Specialist Group (CPSG), which combines in situ with ex situ conservation measures for the optimum protection of a given species (Byers et al. 2013).For the build-up of the conservation breeding programme, the individuals identified as P. variegatus were transferred to the Melinh Station for Biodiversity of the Institute of Ecology and Biological Resources (IEBR), Hanoi.In the Station, located in Vinh Phuc Province, northern Vietnam, besides existing outdoor tank facilities, an exclusive softshell turtle breeding facility was constructed recently.To maximise positive outcomes and for security reasons, viz., to extend the conservation breeding network, another group of the genetically-identified P. variegatus was provided to another softshell breeding facility in northern Vietnam.Successful breeding has already been observed in the colony and offspring are ready for release to the original habitat sites.To extend the conservation breeding programme and, thus, contribute to the build-up of a stable assurance colony and conservation breeding network, a plan has been developed to transfer a limited number of surplus offspring to other facilities in Vietnam and overseas.In addition to these ex situ conservation measures already being in place, focus should now be directed to improving in situ conservation of this beautiful, but threatened softshell turtle species.

Note added in proof
In the late 2023, 50 young and healthy spotted softshell turtles from the in-country breeding program initiated by the Institute of Ecology and Biological Resources (IEBR), Vietnam, together with the Cologne Zoo, Germany, were successfully released to a site in northern Vietnam.

Figure 1 .
Figure 1.Trimmed phylogram, based on the Bayesian analysis.Number above and below branches of major nodes are Bayesian posterior probabilities and ML ultrafast bootstrap values, respectively.Sample highlighted in bold and orange is the paratype of Pelodiscus variegatus.

Figure 2 .
Figure 2. Potential distribution of Pelodiscus variegatus in Vietnam based on A climate B land cover C climate and land cover and D coverage with protected areas as number of suitable grid cells.

Table 1 .
Primers used in this study.

Table 2 .
Variables used for climatic niche modeling computation.NDVI = Normalized Difference Vegetation Index, EVI = Enhanced Vegetation Index.

Table 3 .
Protected areas in Vietnam with projected proper climatic conditions and land-cover (Average Probability > 0.1) for Pelodiscus variegatus, sorted by size of suitable area.A full list of protected areas in China, Laos and Vietnam is presented in the Suppl.material 1.