Biogeography Molecular phylogenetic analyses and ecological niche modeling provide new insights into threats to the endangered Crocodile Lizard (Shinisaurus crocodilurus)

USA. *Correspondence: Minh D. Le, le.duc.minh@hus.edu.vn. Abstract The endangered crocodile lizard, Shinisaurus crocodilurus , is seriously imperiled by anthropogenic threats, including habitat loss and degradation and most critically over-collection for the international pet trade. As a result, population sizes of crocodile lizards have sharply declined throughout their range, with only a small number remaining in China and a handful of individuals left in Vietnam. To prioritize conservation measures for the species, in this study, we generate new mitochondrial sequences of important new samples and analyze them with existing data. Our results confirm a new genetically distinct population in China, highlighting cryptic genetic diversity within the species. The assessment of climate change impacts the to locate the natural distribution range of the newly identified molecular clade in China and determine the distribution extension of the Vietnamese population in the border area, especially potential occurrence on the Chinese side. Considering the impacts of climate change on the Vietnamese population, designing a corridor to connect the subpopulation’s habitat in the border area with nature reserves in Yen Tu Mountain Range and/or translocating lizards from the site to more suitable habitats might help secure the subpopulation in the context of climate change. In all recommended conservation measures, close collaboration between Vietnam and China will be crucial to effectively protect this potentially shared subpopulation of the highly threatened species.


Introduction
Although the endangered Crocodile Lizard (Shinisaurus crocodilurus) is the sole living representative of the monotypic family Shinisauridae, it has been poorly studied until recently. Until 2003, more than 70 years after its description by Ahl (1930), the species was known only from China. Subsequent discoveries of populations in Vietnam extended its known distribution more than 500 km to the south (Le & Ziegler 2003, Hecht et al. 2013, van Schingen et al. 2014a, 2016a. Initial evaluation of genetic and morphological variability of geographically isolated populations based on a limited dataset did not provide a conclusive distinction between the populations (Ziegler et al. 2008). However, a recent study based on more comprehensive molecular, morphological, and ecological data along with ecological niche modeling showed that populations from Vietnam and China represent two separate conservation and taxonomic units, resulting in the description of the Vietnamese form as a distinct subspecies, Shinisaurus crocodilurus vietnamensis (van Schingen et al. 2016b) (Fig. 1). Due to their distinction, the study also recommended the two taxa be managed separately in order to maintain their genetic integrity.
Although morphological and ecological data revealed differences between Chinese and Vietnamese crocodile lizards, molecular analyses with limited data and taxonomic sampling did not clearly resolve the phylogenetic relationships of respective populations (van Schingen et al. 2016b). Specifically, although samples from Vietnam unambiguously formed a reciprocally monophyletic and distinct clade, those from China were split in two independently evolving lineages with one receiving low statistical support values from both Bayesian and maximum likelihood analyses. The results resemble those reported in Huang et al. (2014), which recovered two separate groups of mitochondrial haplotypes in China. Nonetheless, their analysis of microsatellite data and a recent study using genomic data (Xie et al. in press) showed a total of three, instead of two distinct genetic clades. More recent genetic screening of captive individuals in European zoos, however, revealed three mitochondrial clades from China (Ngo et al. 2020), suggesting that earlier analyses could have missed the third clade due to the lack of samples.
Potential distribution of the species was first modeled in Vietnam in 2014 to support field surveys and identify priority areas for conservation measures of this poorly studied species (van Schingen et al. 2014). The prediction of suitable habitat led to the discovery of another population of crocodile lizards in the border area between Vietnam and China (van Schingen et al. 2016a) and highlighted the importance of conducting further surveys in the habitat block, especially in the Chinese side, as the species has never been recorded there. The results also emphasize the need for cross-border collaboration in conserving this highly threatened species in the context of rampant poaching across its range in both countries, which has already led to extirpation in several locations in Vietnam (van Schingen et al. 2015a, 2016a. Li et al. (2013)  In terms of its conservation status, the Crocodile Lizard is facing an increasingly uncertain future. Anthropogenic threats such as habitat destruction through coal mining, dam construction, poaching, and electro-fishing have imperiled populations of the species in both countries (CITES 1990, Huang et al. 2008, van Schingen et al. 2015a, Auliya et al. 2016. Latest estimates showed that populations from China and Vietnam have significantly declined over the last decade with roughly 950 remaining in China and fewer than 150 mature individuals in Vietnam (Huang et al. 2008, van Schingen et al.2015a, 2016a. In addition, populations from China and Vietnam are further subdivided into small-isolated subpopulations (Huang et al. 2008, 2012, van Schingen et al. 2015, 2016a. As a consequence of such small population sizes and increasing threats, the species has been listed as Endangered in the IUCN Red List and included in Appendix I of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) since 2017 , van Schingen et al. 2015a, 2016, CITES 2017.
To better prioritize conservation measures for this critically threatened species, it is important to further clarify its genetic diversity, especially the level of molecular and ecological distinction between populations from Vietnam and China. In order to address these questions, we employed genetic analyses of three mitochondrial genes, cytochrome b, partial ND6, and tRNA-Glu. Ten additional samples from all three known Vietnamese subpopulations and three samples from an unknown clade were sequenced and analyzed to shed light on the unresolved issues related to its population genetics. Moreover, improved species distribution models (SDMs) were applied to assess potential ecological differentiations between different genetically distinct populations that were identified via molecular analyses. We also assessed the impacts of climate change on the populations and identified potential climate refugia, where habitat suitability will likely be maintained under the future climate change scenarios.

Taxonomic sampling and molecular analyses
Thirteen new samples were incorporated in this study, including ten from Vietnam and three with unknown origin from an unidentified clade. Nineteen other sequences published in previous studies were downloaded from GenBank and incorporated in the analyses. New samples from Vietnam were collected during previous field surveys between May and July 2014(van Schingen et al. 2014b, 2015b, 2016a at two sites in Tay Yen Tu Nature Reserve, Bac Giang Province, and Dong Son -Ky Thuong Nature Reserve, Quang Ninh Province, northeastern Vietnam. Crocodile lizards were captured by hand and photographed. Small tissue samples were taken from the tail tip or shading skin, the respective skin parts treated with antiseptic, before animals were released at the same spots. Three taxa, Ctenotus piankai, C. schomburgkii, and Heloderma suspectum were selected as outgroups based on the results of an earlier study (Li et al. 2012) (Table 1).
Total genomic DNA was isolated using Dneasy Blood and Tissue Kit (Qiagen, Germany) for saliva samples and GeneJET Genomic DNA Purification Kit (Thermo Fisher Scientific, Lithuania) for tissue samples. Three fragments of mitochondrial genes, including cytochrome b, partial ND6, and tRNA-Glu, were sequenced using primer 1 and primer 2 from Huang et al. (2014). Tissue samples were extracted using DNeasy blood and tissue kit, Qiagen (Germany). The PCR volume contained 21 μl (10 μl of mastermix, 5 μl of water, 2 μl of each primer at 10 pmol/μl and 2 μl of DNA or higher depending on the quantity of DNA in the final extraction solution). PCR condition was: 95 °C for 5 min to activate the taq; with 40 cycles at 95 °C for 30 s, 50 °C for 45 s, 72 °C for 60 s; and the final extension at 72 °C for 6 min. PCR products were subjected to electrophoresis through 1% agarose gels (Serva, Germany). Gels were stained for 10 minutes in 1X TBE buffer at 2 pg/ml of ethidium-bromide and visualized under UV light. Successful amplifications were purified to eliminate PCR components using GeneJET™ PCR Purification kit (Thermo Fisher Scientific, Lithuania). Purified PCR products were sent to First Base (Malaysia) for sequencing.

Mitochondrial data analyses
Sequences generated in this study were edited using the program Geneious v.7.1.8 (Kearse et al. 2012). The sequences were ClustalX v2.1 (Thompson et al.1997) with default settings. Data were analyzed using the Bayesian inference (BI) as implemented in MrBayes v3.2.1 (Ronquist et al. 2012). We used the optimal model for nucleotide evolution determined using Jmodeltest v2.1.4 (Posada & Crandall 1998) with parameters estimated by MrBayes v3.2.1. Two simultaneous analyses with four Markov chains (one cold and three heated) were run 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. Two independent analyses were run simultaneously. The posterior probability values for major clades in the final majority rule consensus tree were provided.
We also analyzed data using maximum parsimony (MP) as implemented in PAUP*4.0b10 (Swofford 2001) and maximum likelihood (ML) as implemented IQ-TREE v1.6.7.1 ) with a single model and 10,000 ultrafast bootstrap replications. The optimal model for nucleotide evolution was set to TPM2uf+G as selected by Jmodeltest v2.1.4 and applied to both ML and BI analyses. Nodal support was evaluated using Bootstrap replication (BP) as estimated in PAUP*4.0b10 and ultrafast BP in IQ-TREE v1.6.7.1 and posterior probability (PP) in MrBayes v3.2. BP ≥ 70 (Hillis & Bull 1993) and ultrafast BP (UBP) and PP ≥ 95% are regarded as strong support for a clade (Ronquist et al. 2012. Nodes with PP and UBP values between 90 and 94 were considered well-supported. The uncorrected pairwise distance (p) were calculated in PAUP*4.0b10.

Ecological niche modeling
Based on the molecular results supported by previous studies, four genetically distinct populations, three from China and one from Vietnam, were selected for modeling habitat suitability. Although mitochondrial data only recovered two independently evolving lineages in China, microsatellite markers and genomic data clearly revealed three genetically distinct and geographically isolated groups (Huang et al. 2014, van Schingen et al. 2016. On the other hand, the population from Vietnam showed little structure (Nguyen et al. 2017). A total of 41 occurrence records was first split into four genetically distinct clusters for niche suitability modeling. These included three groups from China, consisting of group 1 coded in blue (12 localities), group 2 in yellow (6 localities), and group 3 in green (7 localities) (Fig. 3, Huang et al. 2014) and one group from Vietnam, group 4 in red (16 locations) (Fig. 3, (Fick & Hijmans, 2017). Additionally, we added 12 precipitations, 12 maximum monthly average temperatures, and 12 minimum monthly average temperature layers and other environmental variables, including elevation, aspect, and slope, to the input dataset. The Digital Elevation Model (DEM) was obtained from the global Shuttle Radar Topography Mission (SRTM) (Farr et al. 2007). Values for elevation, aspect, and slope were extracted and interpolated from the Digital Elevation Model (DEM) by ArcGIS 9.2 (ESRI, California).
We performed the tuning process using four feature class combinations (linear, linear and quadratic, linear and quadratic and hinge, and hinge features), and we tested a range of regularization multiplier ranging from 1.0 to 10.0 by increments of 0.5. To identify the best model settings, we adjusted feature types and the regularization multiplier value and obtained the settings that generated models for four populations with lowest average 10% training omission rate. We used ENMeval package in R program to performed this analysis (Muscarella et al., 2014). The optimal model had a regularization multiplier of 3.5 and hinge feature class. Other model parameters, e.g., convergence threshold and background samples, followed re commendations from model developers (Phillips et al. 2006). We employed a geographic partitioning (block method; Muscarella et al. 2014) validation approach for withholding testing data (Shcheglovitova and Anderson 2013). We calibrated mode performance using the area under the receiver operating characteristic curve (AUC). In the context of SDM, the AUC value indicates the ability of the model to discriminate the presence localities from background (Swets & Swets 1988). An AUC value between 0.7 and 0.8, between 0.8 and 0.9, and higher than 0.9 indicates a fair model, a good model, and an excellent model, respectively (Elith et al. 2006). The importance of each variable was assessed based on percent contribution values reported in Maxtent's output files (Bradie & Leung 2017).
For the future climates, we obtained the data layers of the future projections for 2050 and 2080 also from WorldClim 2.1 database (Fick & Hijmans 2017). An ensemble of three general circulation models (GCMs) including BCC-CSM2-MR, CNRM-CM6-1, and MIROC6 under Shared Socio-economic Pathways 245 (SSP 245) was employed. We then applied the same model for 2050 and 2080 but replaced the current climate (2020) with future climate to predict future habitat suitability of the species. The Maxent outputs were classified into four levels of habitat suitability from low to very high. We evaluated the impact of climate change on each category by comparing areas of habitat suitability in the corresponding categories between 2020, 2050, and 2080.

Mitochondrial data analyses
We obtained a final matrix of 1368 aligned characters. In total, 290 characters were found to be Our phylogenetic analyses recovered at least four reciprocally monophyletic clades with one from Vietnam and three from China (Fig. 2). Unlike the results reported in van Schingen et al. (2016b), three clades, except for China Clade 1, received well or strong support from all analyses. China Clade 1 was well corroborated only by the Bayesian analysis (Node 3). However, Node 2, representing all samples of the clade except for the one with Genbank accession number HQ008865, was strongly supported by all three analyses. Sample HQ008865 clustered with Si85 in China New Clade from Reherp with strong support in Ngo et al. (2020), but our phylogenetic analyses showed that all samples recovered in China New Clade originated from Reherp (Fig. 2). The genetic divergence between sequences derived from Vietnamese samples and those from Chinese ones was about 2.2-3.4% based on mitochondrial genes. Three clades from China (Clade 1, 2, and the new clade) were about 2.0-3.7% genetically divergent from each other (Table 2).

Ecological niche modeling
Niche modeling revealed distinct areas of suitable habitats between Vietnamese and Chinese populations (Fig. 3). Maxent models showed reasonable prediction power, and the optimal model for the current potential distribution of the modeled populations had a regularization multiplier of 3.5 and hinge feature class (average test omission rate at the 10% training presence threshold = 0.057, and average AUC value for Vietnam group = 0.997; 3 China group Yellow = 0.998; Green = 0.713; and Blue = 0.956). For each population, environmental parameters contributed differently to its habitat suitability. For the Chinese red population, the aspect accounted for 55.3%, slope 18%, and January precipitation 13.4%; for the Chinese yellow population March precipitation 40.7%, aspect 21%, and slope 15%; for the Chinese green population aspect 38.4%, March precipitation 30%, and February precipitation 13.4%; for the Vietnamese population, aspect 58.3%, slope 13.4%, minimum temperature of the coldest month 9.1% (Table 3).
Under predicted climate change scenarios, the Chinese blue population will apparently become most affected. The distribution range of the population will shrink under all scenarios, except for low habitat suitability. For the three other populations, the potential distribution range will decrease by 2050 and 2080 compared to that of 2020 only under the low habitat suitability category. On the contrary, under other categories of habitat suitability, nearly all populations will gain in range size (Table 4, Fig. 3-5). Distribution ranges of the Vietnamese and China blue populations will become fragmented and/or will shift to other yet unsuitable areas (Fig. 3-5).

Discussion
Our phylogenetic analyses strongly support the presence of four genetically distinct populations with three from China and one from Vietnam. All clades are well supported regardless of the analyses, except for China Clade 1. However, if the sample with GenBank accession number HQ008865 is removed, the clade is also strongly supported in all analyses. The sample was grouped with China new clade in a previous study (Ngo et al. 2020), which used a short gene fragment. Nonetheless, it is distantly related to the samples within China new clade in this study (Fig. 2). The results confirm that data from the new clade were not included in the analyses undertaken by Huang et al. (2014) and van Schingen et al. (2016), as suggested by Ngo et al. (2020).  Although there are some inconsistencies between mitochondrial and microsatellite results, both types of analyses robustly corroborate distinct genetic divergence between populations from Vietnam and China. In addition, morphological, ecological, and behavioral differences between Vietnamese and  (Hu et al.1984, Zhao et al. 1999, Zhu et al. 2002, Ning et al. 2006, Zollweg 2012, Zollweg & Kühne 2013, van Schingen 2014, van Schingen et al. 2015b, Werner 2015. Similar morphological and ecological studies need to be conducted for the lineages in China because they show an equivalent level of genetic divergence ( Table 2). The confirmation of the new clade suggests that genetic diversity of the species has not been  (Ngo et al. 2020). As it represents an independent conservation unit within Shinisaurus crocodilurus, the adult individual is kept separately. The niche model analyses corroborate genetic results based on microsatellite and genomic data, as habitat suitability clearly differs between the Vietnamese and three Chinese populations (Fig. 3). By adding additional environmental variables, our study reveals even clearer niche separation between the populations than that supported by van Schingen et al. (2016b). Their habitat suitability is also increasingly differentiated in future climate scenarios ( Fig. 4 and Fig. 5). While the parameters, slope and aspect, which were not included in the previous study, contribute most significantly to the habitat suitability of all populations, annual mean and range temperatures and precipitation of the driest and warmest months were also found to account for a large proportion of environmental variables within suitable niche as reported in van Schingen et al. (2016).
The assessment of climate change impacts on the crocodile lizard shows that the China Blue population will be most affected with its distribution range becoming highly fragmented and decreased in size in almost all habitat suitability categories from 2050 to 2080. The situation may not be as severe as suggested by Li et al. (2013) and van Schingen et al. (2016a), which predicted that habitat suitability of the species in China will largely vanish by 2080. Based on the present analyses, however, projected distribution ranges of the other two populations in China will remain stable and even expand under impacts of climate change in 2050 and 2080. The marked difference in predicted distribution models between our and previous studies can be explained by the inclusion of other environmental variables, such as slope and aspect, as input data in our modeling analyses. In addition, in this study, we modeled the distribution of genetically distinct populations separately instead of regarding them as a single entity. It is, however, important to note that the models did not incorporate the effect of deforestation, which will likely lead to a further decrease in suitable habitat among all populations of this forest-dwelling species. Three populations in China are not currently well covered by existing protected areas. According to our analyses, only three protected areas, Gu Po Mountain, Da Yao Mountain, and Da Ping Mountain nature reserves, fall within the peripheral portion of the Blue population's range while Yellow and Green populations are not protected by any nature reserve. The future projected scenarios will likely get worse for the Blue population as the highly suitable habitat will move away from the protected areas. The projected distribution ranges under different climate change scenarios continue to reveal that Yen Tu Mountain Range will serve as a refuge for the Vietnamese population in the future, as demonstrated in van Schingen et al. (2016a). Three protected areas, Tay Yen Tu, Yen Tu, and Dong Son-Ky Thuong nature reserves, will play an important role in safeguarding the population long into the future.
Nonetheless, our models suggest that the site near the border between China and Vietnam will become unhabitable under the changing climate. The distribution area, which might extend to China, has been considered a priority region for conservation in the context of climate change by previous studies (van Schingen et al. 2014b, van Schingen et al. 2016a). As the forest on the Chinese side has not been well studied, we recommend that additional field surveys be undertaken in the area to determine the distribution range of the population. International collaborations, involving both Vietnamese and Chinese scientists, will have the highest chance of success in identifying a distribution extension of the population. In addition, as suggested by our results, since the area will likely become unsuitable for the population, it is also crucial to design a corridor connecting the site with Yen Tu Mountain Range, especially with the nearest Dong Son-Ky Thuong Nature Reserve, to facilitate the movement of the subpopulation. Developing new biological corridors in Vietnam is a conservation priority as specified in the incoming National Biodiversity Strategy and Action Plan. This possibility should therefore be further investigated to identify suitable area for the potential corridor. Moreover, translocating the lizards to one of the three nature reserves is another potential conservation measure to safeguard them from impacts of climate change. In both cases, transboundary cooperation between the two countries will be needed to ensure a successful outcome.
Viet Nam Academy of Science and Technology (VAST, Project Code: ĐLSĐ00.01/20-23), the Cologne Zoo, the European Union of Aquarium Curators (EUAC) Cologne Zoo is partner of the World Association of Zoos and Aquariums (WAZA): Conservation Projects 07011, 07012 (Herpetodiversity Research, Amphibian and Reptilian Breeding and Rescue Stations). This study was funded by the Prince Albert II of Monaco Foundation and the Ministry of Science and Technology's Program 562 (grant no. ĐTĐL. ). We thank the participants of the workshop entitled "Consultation and training workshop: Collaborative transboundary conservation of vulnerable species and habitats under climate change" at Hanoi University of Science, Vietnam National University in Hanoi, Vietnam (July 11-12, 2018) and the symposium entitled "Crossboundary Cooperation for Biodiversity Conservation in Asia under Global Change" at Henan University in Kaifeng City, China (29-31 July, 2019) for their comments and suggestions on the study. Comments from two anonymous reviewers and editors helped improve the manuscript.