Giants in the landscape: status, genetic diversity, habitat suitability and conservation implications for a fragmented Asian elephant (Elephas maximus) population in Cambodia

Study area

We focused our study on the Asian elephant sub-populations within the Prey Lang Extended Landscape (PLEL), a large expanse of mostly flat land in central-northern Cambodia west of the Mekong river. The PLEL is characterised by a mosaic of lowland evergreen forests and agricultural land, with elevations below 200 m above sea level. Its climate is tropical, characterised by a wet season from May to November and a dry season from December to April. Within PLEL, the study area (103°57′–105°58′ E and 14°26′–12°34′ N) comprises three protected areas where elephant presence has been recorded over the past decade (Prey Lang Wildlife Sanctuary, Preah Roka Wildlife Sanctuary and Chhaeb Wildlife Sanctuary), as well as one protected area (Kulen Promtep Wildlife Sanctuary) where elephant presence has not been recorded over the past decade and was therefore not sampled (Fig. 1).

Elephant presence

Available elephant occurrence records were collated from relevant stakeholders: Cambodia’s Ministry of Environment and conservation NGOs (Fauna & Flora International, Conservation International, Wildlife Conservation Society and Wild Earth Allies). We conducted non-invasive genetic sampling, involving the collection of dung samples, during the 2020–2021 dry season in Prey Lang, Preah Roka and Chhaeb Wildlife Sanctuaries (see below), in areas where elephant presence had been noted by local rangers. These surveys in turn generated additional elephant presence records for the target protected areas. All occurrence records were then used to inform the habitat suitability modelling.

Population survey

Our design followed survey standards for elephant monitoring through capture-recapture sampling (Hedges & Lawson, 2006; Karanth, Nichols & Hedges, 2012). Survey methods were non-invasive, relying on sample collection from elephant dung with the relevant permit from Cambodia’s Ministry of Environment (permit 248-MoE), with no contact made with elephants. Using faeces for DNA collection is in accordance with the 3Rs principle for animal welfare as it is non-invasive (Zemanova, 2020).

We identified focal areas for faecal collection based on data gathered from previous studies in the landscape (personal communication with locally-active conservation NGOs) and on key informant interviews with one staff from each of three conservation NGOs working in the area and with the head Ministry of Environment rangers in charge of patrolling each of the protected areas within the landscape. The study was government-sanctioned, including the interview phase (permit 248-MoE), and the purpose of the interviews was explained in advance to all interviewees, who participated voluntarily. Survey areas were 625 km² in Preah Rokar/Chhaeb Wildlife Sanctuaries and 1,000 km² in Prey Lang Wildlife Sanctuary.

Samples from elephant dung were collected during the 2020–2021 dry season between 29^th November 2020 and 14^th March 2021 by teams of trained field researchers. Five collection visits were conducted in each wildlife sanctuary and faecal samples were collected via two methods. The first method involved using a cotton swab to collect genetic material from the outside of a dung bolus; this has proved successful at preserving DNA collected from African forest elephant (Loxodonta cyclotis) dung samples (Bourgeois et al., 2019). The second method involved collecting exterior parts of a bolus and preserving in absolute ethanol; this was used as a backup, as the technique provides material for repeat extractions. The second method has previously proved successful for DNA collection from Asian elephant dung samples (Fernando et al., 2000). We followed the methods outlined in these aforementioned studies with minor alterations that can be found in Article S1.

DNA extraction and quality check

Samples were transported immediately after each field trip to the Royal University of Phnom Penh (RUPP) Conservation Genetics Laboratory and stored in a −20 °C freezer until DNA extraction. QIAamp® Fast DNA Stool Mini kits were used for DNA extraction with some initial preparation differences depending on whether the sample was a swab or dung preserved in ethanol (see Article S1 for details). A total of 199 dung piles were sampled, with 155 (78%) being sampled via both collection methods. Therefore, a total of 354 sample tubes were provided to the RUPP Conservation Genetics laboratory between Dec 2020–March 2021. DNA preservation success rates for each collection method were determined using 20 dung piles from which both a swab sample and an ethanol sample were obtained. Comparisons were made between the DNA extractions based on performance in one of the methods required for individual identification (Single Nucleotide Polymorphism-SNP-genotyping).

An attempt to genotype the elephant producing each sampled dung pile was made to individually identify the elephants within the PLEL and therefore provide recapture rates. Using faecal samples to estimate population sizes via a capture-recapture approach using genetic material was first performed in 1999 (Kohn et al., 1999) and has become a widely used tool within conservation, including for estimating elephant population size (Eggert, Eggert & Woodruff, 2003; Gray et al., 2014; Bourgeois et al., 2019). Here we use two genetic marker sets for genotyping, one based on microsatellites and the other on SNPs. These are both presumed neutral marker-sets distributed across the nuclear genome.

SNP genotyping

Twenty SNPs have previously been screened for use in Cambodian elephant populations. All SNPs were biallelic and identified based on high minor-allele frequencies. The genotyping method used KASP assays (LGC Genomics, Teddingto, Herts, UK), with fluorescence after PCR amplification measured on a StepOne™ Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). DNA extracts were diluted 1:5 to decrease inhibition prior to amplification. Included on every 48-well plate were two negative controls and 2–3 positive controls (blood samples) with known variability in genotypes. Genotyping was only conducted after checking that the results of the controls were successful. The PCR protocol can be found in Article S1 and details of the KASP assays can be found in Table S1.

Microsatellite genotyping

DNA extracts were diluted 1:5 and genotyping was attempted at nine microsatellite markers in two panels. These microsatellites had been previously designed for use in elephants and have been used to study genetic diversity in Asian elephants (Kongrit et al., 2008; Comstock, Wasser & Ostrander, 2000; Nyakaana & Arctander, 1998). All microsatellite genotyping was carried out in triplicate and a consensus genotype was created for each sample if the allele was observed ≥2 times. A negative and positive control (blood sample) were included on each plate. See Article S1 and Table S2 for details of the microsatellites used in this study, including primer sequences, fluorescent tags and PCR protocols used for multiplexing.

Population structure

The nuclear markers (SNPs and microsatellites) were used to test for population structure in the Cambodian samples using three methods. Firstly, a Bayesian method using the program STRUCTURE (Pritchard, Stephens & Donnelly, 2000) was performed using an MCMC of 1,000,000, an initial 250,000 were run and discarded as burn-in. An admixture model (infer alpha) and correlated allele frequencies model (Lambda = 1) was used for a K of 1–5. All runs were performed in triplicate to verify the stability of the results. Secondly, a principle component analysis (PCA) was implemented using the dudi.pca command in the ‘adegenet’ package in R Studio (Jombart, 2008). Lastly, the fixation index (FST) was calculated using hierfstat (Goudet, 2005) using the Weir and Cockerham method (Weir & Cockerham, 1984), values were compared to the divergence levels described in the original method development (Wright, 1965; Hall, 2022). All three analyses were performed on the unique genotypes identified via the individual identification methods outlined below, however unique genotypes that were scored at <4 microsatellites were removed (n = 8), to reduce the amount of missingness. Any markers that showed evidence of null alleles, according to the program Cervus v3.0.7 (Kalinowski, Taper & Marshall, 2007), were also removed before analysis. We used three methods to reduce the risk of incorrect assessment, as multiple studies have highlighted assumptions and biases that can result from population genetic approaches (Puechmaille, 2016; Loog, 2021; Elhaik, 2022).

Mitochondrial sequencing

A 152bp mitochondrial sequence was produced using primers designed specifically for use in degraded Asian elephant faecal samples. Due to the degraded nature of DNA extracted from elephant faecal samples and the potential of amplifying non-target nuclear mitochondrial (NUMT) DNA, the mitochondrial sequence was created using two separate primer pairs that were each designed to amplify a small section of an overlapping region of the mitochondrial d-loop. The primers used in the PCR for each fragment were designed using Primer 3 (Untergasser et al., 2012) and a section of the Asian elephant mitochondrial genome sequence (Accession: DQ316068.1). The fragment 1 primers (AEL_dloop) were designed to amplify a 204bp sequence and the fragment two primers (CR_AEL_RZSS) to amplify a 193bp sequence, the overlapping region between them was 111bp. Stringent quality checks and end trimming was conducted in Geneious Prime v 2021.1.1, hence the final alignment being smaller than each target mitochondrial product. This target region was chosen as it is a highly variable section of the mitochondrial genome that exhibits variation between maternal lineages. This section of the control region also overlaps sequencing performed in Asian elephant samples from another eight countries (Vidya, Sukumar & Melnick, 2009), allowing comparisons to be made between the PLEL population and others. Details of the primers and PCR protocols are provided in Article S1 and Table S2. A haplotype network was created in the program PopArt v1.7 using the median joining network algorithm (Bandelt, Forster & Röhl, 1999).

Elephant sexing

The method developed by Ahlering et al. (2011) was used to sex each sample. This technique uses three primer sets to amplify regions of three genes on the sex chromosomes, 1 on the X chromosome (PLP) and 2 on the Y chromosome (AMELY and SRY). As female elephants are homogametic and male elephants are heterogametic, the three genes will only amplify in male elephants, with only one amplifying in females (PLP). The PCR amplification followed the protocol of Ahlering et al. (2011) and products were visualised on a 2% agarose gel stained with GelRed®. A negative and two positive controls (known male and female) were included on every PCR plate. Any samples that showed unclear band separation were tested a second time (n = 8).

Genetic diversity

Two approaches to measuring genetic diversity were used in this study, one using the mitochondrial data and the other using the nuclear microsatellite markers. Haplotype diversity (Hd) and nucleotide diversity of the mitochondrial sequences were calculated in DNASP v6.12.03 (Rozas et al., 2017). As this same mitochondrial region has been sequenced in elephants from other countries (Vidya, Sukumar & Melnick, 2009) a 152bp alignment was created using the geneious alignment algorithm in Geneious Prime v 2021.1.1 for diversity comparisons between 15 populations in nine countries. The nuclear diversity was calculated using only the individually identified elephants which were genotyped at ≥4 microsatellites. Only the microsatellites which did not show any evidence of null alleles were included. Observed heterozygosity (Ho), Expected heterozygosity (He), Rarified Allelic richness (Ar) and the Fixation index (FIS) were calculated in the R package hierfstat (Goudet, 2005). Given the notorious difficulty of standardising microsatellite genotyping between studies and laboratories (Hoffman & Amos, 2005; Moran et al., 2006; Ellis et al., 2011) comparisons of the results to other studies were not conducted.

Individual identification

The SNP and microsatellite genotyping data were combined to perform individual identification of each dung sample. Probability of identity scores (PID and PID_SIB) were calculated for the marker-sets in the program Cervus v 3.0.7 (Kalinowski, Taper & Marshall, 2007). Only a marker-set with a PID_SIB score of <0.01 was deemed acceptable for analysis. This led to the inclusion of dung samples that were either i) genotyped at ≥18 SNPs or ii) genotyped at ≥15 SNPs AND ≥2 microsatellites. Individual identification was then performed in the R package Allelematch (Galpern et al., 2012) using an allele mismatch score of 6. Any mismatches were manually checked and groups amalgamated into unique genotypes if they shared multiple matches where >50% of mismatches could be attributed to allelic dropout.

Population size estimates

Population size estimates were made from the unique genotypes generated by the genetic analysis results. The R package Capwire, based on the models presented by Miller, Joyce & Waits (2005), was utilised to estimate population size. Capwire has been used in a range of population estimation studies and is recognised within the field (Lonsinger et al., 2019; Davidson et al., 2014). The Two-Innate Rates Model (TIRM) model was run as it accounts for heterogeneity of capture probabilities for different individuals. The model assumes a closed population and the estimate is made from the capture class (the number of times individuals have been captured i times) and the frequency of each capture class.

Habitat suitability and connectivity modelling

A literature review of available studies identified forest cover (Lin et al., 2008; Gray et al., 2014), water accessibility (Sharma et al., 2020), and anthropogenic disturbance as key variables affecting habitat suitability (Gray & Phan, 2011; Huang et al., 2019; Sharma et al., 2020; Neupane et al., 2019). Data was sourced to represent these variables in suitability models: Global Forest Watch for forest cover and deforestation data (Hansen et al., 2013); Open Development Cambodia dataset for rivers (Open Development Cambodia, 2024); and Open Street Map data for roads and settlements (OpenStreetMap Contributors, 2024). Studies of Asian elephants in Cambodia are underrepresented in the academic literature, with a geographic bias towards China and India (Lin et al., 2008; Liu et al., 2016; Huang et al., 2019). In this study, we represent anthropogenic disturbance by road kernel density, Euclidean distance to villages, and recent and historic deforestation prepared in ArcGIS Pro v2.8 (Esri, 2024). Data sources and preparation can be found in the Supplemental Materials (Table S3).

MaxEnt

MaxEnt models were run to provide variable weighting importance for the weighted fuzzy habitat suitability models (Phillips, Anderson & Schapire, 2006). MaxEnt additionally provides insight into areas of suitable habitat by considering the relationships between variables and occurrence locations. Models were run using the default parameter settings on a 10-fold cross-validation model with a bias file included to minimise overfitting caused by sampling bias and a maximum of 5,000 iterations permitted (Kramer-Schadt et al., 2013; Syfert, Smith & Coomes, 2013). The bias file was created in SDMToolbox V2.4 (Brown, Bennett & French, 2017). Results were then evaluated using the Area Under Curve (AUC) to assess appropriateness for use in conservation management.

Fuzzy suitability model

A fuzzy habitat suitability model was used to incorporate local conservation expert opinions on the relationship, either positive or negative, between variables identified in the literature that affect habitat suitability including elephant relationships to forest cover, areas of deforestation, rivers, villages and road density. A fuzzy approach was taken to account for uncertainties in specific thresholds outlined by local experts for the variables in ranges in this understudied area (Ouellet et al., 2021). The model was run in ArcGIS Pro v2.8 using the Fuzzy Memberships tool (Esri, 2021) to apply linear memberships to all variables and the Weighted Overlay tool to combine layers to identify suitability. The application of linear memberships to variables enabled incorporation of uncertainty within expert opinions. MaxEnt results were then used to weight variables in accordance with their contribution to model performance as proxy for importance for elephant habitat suitability. Details on fuzzy memberships and weightings can be found in the Supplemental Materials (Table S3).

The total suitable habitat for Asian elephants was identified by aggregating the fuzzy and Maxent outputs. Fuzzy output areas of >5 were deemed suitable (Aini et al., 2017). Maxent outputs were reclassified using the maximum test specificity plus sensitivity threshold (Liu, Newell & White, 2016). The binary maps were then added to form a comprehensive map of habitat suitability.

Connectivity

Morphological Spatial Pattern Analysis (MSPA) was run using the prepared forest canopy layer to identify core areas of forest. A least-cost path model, adjusted to apply circuitscape theory to identify pinch points in connectivity, was run in Linkage mapper between core areas >25 km² identified in MSPA. A resistance surface was created using Gnarly Landscape Utilities version 0.1.0 (Shirk & McRae, 2013), non-forest areas with <15% canopy cover, settlements and roads, were designated with the highest resistance scores (Table S4) (Suksavate, Duengkae & Chaiyes, 2019; de la Torre et al., 2019). As one of the first studies to analyse Asian elephant connectivity in PLEL the default parameter settings were used including a 200,000 cost weighted distance truncation factor.

Population study

Population structure

The Principal Component Analysis (PCA), the Bayesian STRUCTURE analysis and the fixation index (FST) suggest little to no genetic divergence between the elephants in Preah Roka/Chhaeb and the elephants in the Prey Lang Wildlife sanctuaries. Two markers (SNP023 & FH48) showed evidence of null alleles and FIS values >0.5, and were removed prior to analysis (Table 1). The structure results show no evidence of genetic differentiation between the two areas at K = 3 (Fig. 2 and Fig. S2). The PCA analysis shows a large overlap between the elephant genotypes in the north and south and the low eigen values only explain a small amount of variance, the eigen values produced for PC1 (0.896) and PC2 (0.796) explain 14.2% and 12.6% of the variation, respectively (Fig. 3). The FST value between the elephants found in Prey Lang and those in Preah Roka/Chhaeb was 0.066, indicating a moderate level of divergence. Although these results suggest no, to low divergence, the lack of evidence of any elephants moving between Prey Lang and the other two areas, and the long generation times of elephants, caused us to conduct the remaining analyses separately assuming two populations, Preah Roka/Chhaeb and Prey Lang.

Table 1:

Summary statistics of the genetic markers (SNPs and microsatellites) genotyped in elephant samples collected in Preah Roka/Chhaeb and Prey Lang Wildlife sanctuaries.

Only unique genotypes as identified via Allele-match were included in the analysis and individuals genotyped at <4 microsatellites were excluded. All parameters calculated in Cervus v 3.0.7 except FIS which was calculated in Fstat v 2.9.3.

Locus ID	k	N	HObs	HExp	PIC	F(Null)	FIS
SNP002	2	26	0.5	0.419	0.326	−0.0983	−0.194
SNP004	2	25	0.24	0.216	0.189	−0.0588	−0.096
SNP005	2	25	0.36	0.301	0.252	−0.0969	−0.292
SNP010	2	26	0.615	0.462	0.35	−0.1524	−0.327
SNP011	2	27	0.37	0.492	0.366	0.1318	0.242
SNP013	2	27	0.37	0.307	0.256	−0.1002	−0.214
SNP016	2	26	0.462	0.462	0.35	−0.0097	−0.003
SNP021	2	27	0.111	0.107	0.099	−0.0194	−0.163
SNP023	2	26	0.154	0.462	0.35	0.4927	0.67
SNP025	2	27	0.222	0.307	0.256	0.1518	0.184
SNP026	2	26	0.346	0.382	0.305	0.04	0.081
SNP027	2	27	0.593	0.453	0.346	−0.1429	−0.299
SNP029	2	27	0.111	0.107	0.099	−0.0194	−0.053
SNP031	2	27	0.593	0.453	0.346	−0.1429	−0.308
SNP033	2	26	0.577	0.491	0.366	−0.0901	−0.365
SNP035	2	27	0.148	0.14	0.128	−0.0302	−0.087
SNP036	2	27	0.37	0.352	0.286	−0.0345	−0.073
SNP037	2	27	0.37	0.391	0.31	0.0182	0.03
SNP039	2	27	0.296	0.307	0.256	0.0092	−0.02
SNP040	2	27	0.296	0.307	0.256	0.0092	0.055
Emu04	5	27	0.556	0.655	0.578	0.0797	0.035
FH48	2	27	0.074	0.307	0.256	0.6033	0.657
Emu07	6	26	0.615	0.821	0.777	0.1336	0.213
Emu03	2	17	0.294	0.401	0.314	0.1392	0.279
Emu12	3	27	0.481	0.398	0.352	−0.1271	−0.57
Emu10	5	25	0.6	0.749	0.687	0.0981	0.17
Emu17	5	14	1	0.786	0.719	−0.1481	−0.327
Emu15	6	18	0.778	0.663	0.615	−0.114	−0.222

DOI: 10.7717/peerj.18932/table-1

Note:

K, number of alleles; N, individuals genotyped; HObs, Observed heterozygosity; HExp, Expected heterozygosity; PIC, Polymorphic information content; F(Null), Null allele frequency; FIS, fixation index as calculated in Weir & Cockerham (1984) for small f.

Mitochondrial sequencing

A total of 87/163 (53%) of the samples were sequenced at a 152bp sequence of the mitochondrial control region and these sequences were aligned and compared to the 534 available sequences from other range states. A total of 21 haplotypes were observed across the global sample set and clustered into two Asian elephant clades (α and β), and the PLEL elephant samples contained four haplotypes (Fig. 4). One haplotype (I) was in the β-clade and has been observed previously in India, Sri Lanka and Vietnam. This haplotype was only observed in samples collected in the Preah Roka/Chhaeb region. The other three haplotypes were in the α-clade and have previously been observed in other countries (Fig. 4). One of these (IV) was identified in both populations in this study; however, two (II & III) were just observed in samples from the Prey Lang Wildlife sanctuary. The four haplotype sequences identified in the PLEL have been uploaded to NCBI genbank with the following accessions (PP066116–PP066119), but it should be noted that they each align with haplotypes identified from longer mitochondrial sequences previously published as part of a range-wide study.

Genetic diversity

The mitochondrial sequences of the PLEL samples were compared to 13 previously sampled populations in eight countries. They reveal that the two PLEL populations have average levels of haplotype diversity (Prey Lang: 0.55 ± 0.05 and Preah Roka/Chhaeb: 0.46 ± 0.07). However, the two populations in this study differed in their mitochondrial nucleotide diversity. Of note was the nucleotide diversity observed in the Preah Roka/Chhaeb population (0.034 ± 0.005) which was higher than all other sampled populations globally except Sri Lanka (Fig. 5). This is striking given that only two haplotypes were observed in this population; however, they are from highly differing evolutionary histories, one being from the β clade and the other from the α clade (see Fig. 4). This contrasts with the Prey Lang elephants, which only exhibit mitochondrial lineages from the α clade and thus have lower levels of mitochondrial nucleotide diversity (Fig. 4). The nucleotide diversity (microsatellites) did not differ between Prey Lang and Preah Roka/Chhaeb at any of the genetic diversity measures (Table 2).

Table 2:

Genetic diversity calculations using the microsatellite genotyping of individually identified elephants in the Prey Lang Extended Landscape (PLEL).

Only individuals genotyped at ≥4 microsatellites were included. All parameters calculated using hierfstat and standard errors are included in brackets.

Population	N	HObs	HExp	FIS	AR
Prey Lang	17	0.584 (0.074)	0.627 (0.066)	0.069 (0.105)	3.17 (0.376)
Preah Roka/Chhaeb	9	0.571 (0.104)	0.639 (0.063)	0.106 (0.134)	2.95 (0.299)

DOI: 10.7717/peerj.18932/table-2

Note:

N, number of individually identified elephants included; Hobs, Observed heterozygosity; HExp, Expected heterozygosity; FIS, fixation index as calculated in Weir & Cockerham (1984) for small f; AR, Allelic richness.

Habitat models

MaxEnt and fuzzy habitat suitability

Habitat suitability models indicate habitat preferences, mapping potential suitable habitat across the PLEL. The maxent model reached an AUC of 0.941 across the 10 cross-validated models, exceeding the accuracy threshold recommended for application in conservation management. Distance to historic deforestation had the greatest effect on suitability with a 44.7% contribution to model performance, followed by distance to villages (24.4%), and distance to rivers (15.5%). The ensemble Maxent and Fuzzy output identified 276.06 km² of potentially suitable Asian elephant habitat, the majority of which is found within Prey Lang Wildlife Sanctuary and to the north of Preah Roka Wildlife Sanctuary and Chhaeb Wildlife Sanctuary (Fig. 6).

Connectivity

Connectivity models highlight areas of potentially constricted movement across the landscape. High-cost areas in which movement is restricted are found to the north-west of Prey Lang, south-west of Preah Roka and south-east of Chhaeb wildlife sanctuaries. Similarly, to the south and east of Prey Lang Wildlife Sanctuary. Pinch points demonstrate where landscape connectivity would be disproportionately affected if these pathways were lost are also demonstrated connecting all focus protected areas (Fig. 7). Pinch points are highlighted by corridors of medium connectivity north of Prey Lang Wildlife Sanctuary towards Chhaeb Wildlife Sanctuary. Pinch points are also found within Preah Roka Wildlife Sanctuary as movement is restricted by forest degradation within the sanctuary.

Population size and fragmentation

Using genetic data, we estimated the Asian elephant population in the Prey Lang Extended Landscape, adding to the body of evidence demonstrating the value of non-invasive capture-recapture methods for estimating the size of relatively small populations of elusive wildlife. Understanding population size is key to assess the conservation status of species and to identify management actions needed to secure long-term population persistence (Chaudhary & Oli, 2020). This has been highlighted as increasingly important for Asian elephants, with most population estimates based on educated guesses (Blake & Hedges, 2004). In contrast, our study provides a reliable estimate that can inform management decisions in Cambodia.

The Prey Lang Wildlife Sanctuary elephant population is estimated to number 31 individuals, with a 95% confidence interval of 24 to 41 individuals. The Preah Roka/Chhaeb Wildlife Sanctuaries population is estimated to number 20 individuals, with a 95% confidence interval of 13 to 22 individuals. These estimated population sizes are relatively small but somewhat larger than had been anticipated based on local knowledge gathered to identify survey areas. Minimum viable population sizes in vertebrates depend on factors such as population growth rate (Reed et al., 2003). In the case of elephants, which are severely impacted by habitat fragmentation (Leimgruber et al., 2003; Ripple et al., 2017), it has been proposed that for a population to be viable, it should be >500 elephants (Santiapillai, 1997). This has been suggested to avoid the most severe consequences of inbreeding but does not eliminate them and is unlikely to retain genetic diversity needed for the population to adapt to future challenges, such as climate change or disease (Frankman, Bradshaw & Brook, 2014). Over 1,000 individuals may be necessary on timescales longer than 100 years (Sukumar, 1993) and viability will also depend on the population’s vital rates such as female mortality (De Silva & Leimgruber, 2019). Therefore, the population estimates produced in this study show that the elephants in this landscape are in an extremely precarious situation, with numbers far lower than needed for continued population survival.

The genetic results from this study suggest that the elephants in the landscape are fragmented into two populations. Evidence for the movement of six elephants between Preah Roka and Chhaeb wildlife sanctuaries was produced from the genetic survey, with unique genotypes from the same elephants being found in both protected areas. However, none of the 35 unique genotypes were shared between Prey Lang and the Preah Roka/Chhaeb sanctuaries. This result, combined with the genetic differences observed in the mitochondrial lineages, the population structure analysis and the large barriers to movement observed to the northwest of Prey Lang, suggests that fragmentation into two extremely small populations has occurred recently. These results add to the body of evidence (e.g., Lowe & Allendorf, 2010) on how genetic data, in combination with habitat analyses, can be used to assess fragmentation in wildlife populations. It is therefore of vital importance to (i) increase population sizes as rapidly as possible and (ii) reconnect these populations with each other and/or other neighbouring populations (Johnsingh & Williams, 1999).

This study has identified 27,606 hectares of suitable elephant habitat within the PLEL, the majority of which is within Prey Lang Wildlife Sanctuary. However, as elephants are a long-lived species, with long gestation times, and single calves and can take 13 years to reach sexual maturity (De Silva et al., 2013), there are many constraints that limit rapid population growth. Additionally, potential carrying capacity of this area is difficult to predict as Asian elephant densities vary considerably between previous populations studied, from 0.053 elephants per km² in China (Sun et al., 2021) to 1.7 elephants per km² in India (Kumara et al., 2012). Studies performed in populations geographically closest to Cambodia, such as Thailand, have produced densities between 0.1 and 0.7 elephants per km² (Htet et al., 2021). This would predict that the carrying capacity for the remaining suitable elephant habitat (27,606 ha) could roughly lie between 27–193 elephants. Given this uncertainty we advise that it is a priority to identify any factors that may be limiting population growth (and remove them if possible) and identify whether the area can support a larger population. Even if the population cannot grow rapidly, ensuring that it can grow is vital to retain the genetic diversity that remains.

Another way of maintaining, and potentially increasing, genetic diversity is by reconnecting small fragmented populations with others. This can consist of restoring natural corridors between populations to facilitate movement or by meta-population management via regular human-mediated translocations between populations (Hoffmann, Miller & Weeks, 2020). Given the difficulties and risks (to both elephants and humans) involved in elephant translocations (Santiapillai, 1997) creating natural corridors is often the most suitable. The large separation that now exists between Prey Lang and Preah Roka/Chhaeb also makes this option challenging, and while elephants can thrive in degraded habitat, a significant concerted effort by government and local stakeholders would be required to establish passages across infrastructure barriers and to manage potential human-elephant conflict.

From a methodological standpoint, we encountered high amounts of DNA degradation that needed to be overcome for accurate genotyping, similar to previous attempts to conduct genetic population estimates using faecal samples (Costa et al., 2017; Taberlet et al., 1996). To increase our ability to determine unique genotypes we have combined microsatellite and SNP datasets. Alongside (Crouthers, 2024), this is the first time this has been conducted in an Asian elephant study. We also tested two collection/extraction techniques and showed that samples collected via the swabbing technique produced more suitable DNA for genotyping than the material stored in ethanol. The sample success rate for the SNP genotyping was 65%, which is a substantial improvement on the previous sampling of elephants in a dung survey in the Cardamom Mountains during 2015–2016, where the sample success rate for microsatellite genotyping was 30% (unpublished data). However, this improvement is likely related to the markers being genotyped rather than the collection technique. As in Taberlet, Waits & Luikart (1999), we advise that investigation is conducted to determine the most appropriate sampling, extraction and genotyping technique for a given study, as variation between regions and studies can be high (Piggott Maxine, 2004; Broquet, Ménard & Petit, 2007; Bourgeois et al., 2019).

Genetic diversity and structure

A positive result for both populations considering they are so small was that mitochondrial haplotype diversity were still average when compared to 13 populations in eight other range states (Vidya, Sukumar & Melnick, 2009). Additionally, the nucleotide diversity for Preah Roka/Chhaeb was unusually high, driven by the presence of two highly divergent mitochondrial lineages, one from the α-clade and one from the β-clade (Fernando et al., 2000). These two clades are thought to have evolved independently in the north and south of the Asian elephant range respectively (Fernando et al., 2000; Vidya, Sukumar & Melnick, 2009) but have been found to exist together in some populations, such as those in Sri Lanka and Myanmar (Vidya, Sukumar & Melnick, 2009). Several hypotheses have been proposed to explain their co-occurrence including natural migration and human-mediated transport of elephants, as Sri Lanka was historically a trading hub for captive elephants (Fleischer et al., 2001; Sukumar, 2003). The presence of both haplotype clades in Preah Roka/Chhaeb and not in Prey Lang could be due to random stochasticity in these small populations or increased historical connectedness of the former to other elephant populations (natural or human mediated). Other populations that they could have been naturally connected to would have been across the Thai border or potentially more recently to the remaining strongholds of elephants in Cambodia. Conservation efforts have focused on these strongholds in southern (Cardamom mountains) and eastern (Eastern Plains Landscape) Cambodia (Maltby & Bourchier, 2011). Currently only one scientific article documenting the use of a genetic faecal population estimate has been produced in Cambodia (Gray et al., 2014). It was conducted in 2009 and estimated that Phnom Prich Wildlife Sanctuary, in the Eastern Plains Landscape, contained 136 ±18 (SE) elephants (Gray et al., 2014). In future, it would be useful to gain more regional estimates to put our results within a national context for country-wide conservation planning of this iconic species. Genetic diversity is one of the pillars of biodiversity and the use of genetic diversity in establishing conservation priorities has been strongly supported (DeWoody et al., 2021). We believe our study presents a useful case study where genetic diversity is considered alongside population size and habitat suitability and fragmentation to provide a wholistic understanding of the status of a threatened vertebrate as a basis to inform management and conservation.

Although in this study genetic diversity comparisons with published studies were based on mitochondrial sequences, this study also genotyped nuclear microsatellites that allowed nuclear diversity comparison between the Preah Roka/Chaeb and Prey Lang elephants. No difference in genetic diversity was found, however the ability to compare nuclear diversity across populations at a greater scale throughout Asia would be extremely beneficial, as there are several studies in other range countries (e.g., Zhang et al., 2015; Parida et al., 2022). Microsatellite datasets are notoriously difficult to standardise across studies (Hoffman & Amos, 2005; Moran et al., 2006; Ellis et al., 2011); however, there are multiple strategies that can be used to gain more robust nuclear diversity comparisons in the future. One is to pursue a genotyping method that can be more easily standardised, such as SNP competitive assays (Vignal et al., 2002; Semagn et al., 2014; von Thaden et al., 2017), an alternative approach which alleviates some but not all problems with microsatellite comparisons is the ‘yardstick’ approach (Skrbinšek et al., 2012). This has recently been attempted for Asian elephants using a population in Laos (Budd et al., 2023). This ‘yardstick’ population has been genotyped at the same seven microsatellites analysed in this study. Using the dataset of this population and removing individuals genotyped at <4 of these microsatellites we calculated the rarefied Allelic Richness for the PLEL and the Laos ‘yardstick’ population. The diversity ratio was 0.840 (±0.139 SE). This suggests that the elephants in the PLEL are less diverse than the ‘yardstick’ population, however 88% (n = 28) of populations compared to the ‘yardstick’ have had lower diversity ratios than this one (Budd et al., 2023). This level of nuclear genetic diversity in the PLEL seems to corroborate (Budd et al., 2023)’s findings that elephant populations in Southeast Asia, such as the PLEL, retain unusually high levels of diversity compared to populations elsewhere.

We found an even sex ratio in the samples as a whole; however more males than females were recorded in Prey Lang but a greater number of females compared to male elephants was recorded in Preah Roka/Chhaeb. This could be due to stochasticity in these small populations, but we advise that sex ratios are monitored in any future surveys. A severely skewed female-bias sex ratio can be indicative of elephant populations that have been subjected to high levels of ivory poaching (Sukumar, Ramakrishnan & Santosh, 1998) as males are targeted for their tusks (Mondol, Mailand & Wasser, 2014). This result could therefore indicate a higher poaching threat in the Preah Roka/Chhaeb population compared to the Prey Lang region, however the sample size is too small and would require further investigation before any conclusions are made.

Identification of the relatively high mitochondrial genetic diversity provides some positivity for the elephant population in the PLEL. Typically, populations that have undergone extreme population declines lose substantial genetic diversity, indicating that diversity levels may have been even higher in the past. The remaining high levels of genetic diversity warrant conservation and could aid the prioritisation of the landscape as a national stronghold for this flagship species. However, as genetic diversity is lost more quickly from small populations (Frankham, Ballou & Briscoe, 2010), the conservation need is immediate and requires substantial effort to rescue this remaining diversity.

Habitat suitability and connectivity

Habitat suitability results show extensive areas of suitable Asian elephant habitat in Prey Lang, Preah Roka, and Chhaeb Wildlife Sanctuaries, providing a spatial indication of key areas to focus conservation effort for Asian elephants in the PLEL. However, limited or no recent elephant records exist from some areas modelled as highly suitable habitat, particularly in northern Chhaeb Wildlife Sanctuary on the border with Laos PDR and in central-western Prey Lang Wildlife Sanctuary (Maltby & Bourchier, 2011). Observed discrepancies may be a result of uneven survey effort resulting in an incomplete picture of actual occurrence, and/or limitations in the modelling (Araújo et al., 2019), such as factors not fully accounted for or weighted for, e.g., levels of current or historical human disturbance, access to water sources or seasonal movements. Further research would be required to determine the reasons why occurrence records are concentrated in particular areas of the landscape but apparently absent in other areas with seemingly suitable habitat, so population recovery plans can be better informed. Connectivity results indicate low cost to movement within central regions of protected areas similarly indicating the critical role they play in facilitating species connectivity, and therefore metapopulation sustainability, within the Prey Lang extended landscape. Pinch points in connectivity from northern Prey Lang to Chhaeb Wildlife Sanctuary further demonstrate the importance of what little connectivity there is between the north and south of the PLEL. Currently, this route does require a road crossing. Whilst elephants have been documented to cross roads (Wadey et al., 2018), efforts to decrease resistance to movement in this area through restoration of secondary forests, open habitats and food resources (Wadey et al., 2018) could enhance the viability of road crossing to maintain connectivity. Restoration of forest cover between the two sanctuaries could play a critical role in facilitating the movement of individuals from the two sub-populations identified from the genetic analysis. Despite the areas between Prey Lang and Chhaeb not being designated as protected areas, adequate management could support the re-establishment of vital connectivity (Huang et al., 2019).

This study provides useful insights into areas of suitable habitat and possible connectivity for Asian elephants in the PLEL that will help prioritize conservation efforts. Following on from results of this study, future research should consider the key factors identified here, forest cover, deforestation and proximity to rivers and villages, which align with key variables found in other regions (Huang et al., 2019; Sharma et al., 2020), to develop more locally specific models to improve upon this initial assessment. It is recommended that the relationship between habitat suitability and distance to rivers is explored in more depth to take into account fluctuations during wet and dry seasons which may influence seasonal patterns of suitability and connectivity. This study has shown the importance of forest cover for Asian elephant suitability, and future research could consider generating a bespoke landcover map to overcome the biases and limitations faced in this study using a global forest cover dataset (Shimizu, Ota & Mizoue, 2020). More broadly, our study may serve as a model to demonstrate an integrated approach to threatened species assessments to inform management, combining habitat suitability and connectivity analyses with non-invasive methods to determine population size, genetic diversity and fragmentation.