Conservation prioritisation through genomic reconstruction of demographic histories applied to two endangered suids in the Malay Archipelago

The biodiversity of the Malay Archipelago is the product of the region's rich biogeographical history with periods of island connectivity and isolation during the Pleistocene glacial cycles. Here, the case of two endemic suid species, the Javan (Sus verrucosus) and Bawean (S. blouchi) warty pigs, was used to illustrate how biogeographic processes and recent anthropogenic pressures can shape demographic histories with significant implications for species conservation.


| INTRODUC TI ON
The Malay Archipelago, also known as the Indo-Australian Archipelago, continues to serve as a cornerstone for the field of biogeography (Lohman et al., 2011). Ever since the seminal work by A.R. Wallace (1860Wallace ( , 1863, the distribution and dynamics of biota across Sundaland (comprising the Malay Peninsula, Borneo, Sumatra, Java, Bali and smaller islands) and Wallacea (Sulawesi and smaller western Indonesian islands; Figure 1a) inspired the formulation of theories on speciation, dispersal and vicariance. The regional biodiversity is the product of its rich geological history and the drastic sea level changes of the Pleistocene glacial cycle that affected island connectivity (Van Den Bergh et al., 2001). The vast land bridges of the Sunda shelf during glacial maxima facilitated the dispersal of biota across Sundaland (Husson et al., 2020;Sathiamurthy & Voris, 2006;Voris, 2000). Conversely, rising sea levels during interglacial periods enforced the isolation of islands thereby disrupting gene flow and fostering allopatric speciation through vicariance (Brown et al., 2013).
Sundaland now features among the world's leading biodiversity hotspots, with 15,000 plant species and 700 vertebrate species documented as endemic (Myers et al., 2000). With the advent of genetic and genomic approaches, these numbers will stand to rise, as even in well-studied taxa like mammals, previously undescribed or cryptic species are being identified (e.g. Wilting et al., 2007).
The region is facing immense anthropogenic pressures from deforestation, landscape fragmentation, introduction of nonnative species and over-exploitation of natural resources (Myers et al., 2000). An estimated 70% of Sundaland has already been heavily modified by humans, with a significant increase of human pressures reported even across protected areas (Verma et al., 2020).
Restricted geographic ranges and small effective population sizes make insular endemic species particularly prone to extinction risk (Frankham, 1998). In insular systems, most extinctions have been caused by anthropogenic impacts, rather than by environmental or stochastic changes (Wood et al., 2017). However, the impact of recent or ongoing pressures on species persistence can be difficult to interpret without a good understanding of their biogeographical history (Whittaker et al., 2005).
While past cyclical climate fluctuations have led to the largest diversity of wild pig species (suids) on Sundaland and Wallacea (Frantz et al., 2013), most suid species are currently threatened by increasing anthropogenic pressures (Blouch, 1995). Among the five wild pig species of the genus Sus in the Indonesian Archipelago (Figure 1b), only the Indonesian subspecies of the Eurasian wild boar, the banded pig (Sus scrofa vittatus), is thriving and widely distributed across mainland Southeast Asia and Indonesian islands, including Java and Sumatra. Meanwhile, the bearded pig (Sus barbatus) and Sulawesi warty pig (S. celebensis) are listed as 'Vulnerable' (Luskin et al., 2017) and 'Near Threatened' (Burton et al., 2020), respectively, on the International Union for Conservation of Nature (IUCN) Red List of Threatened Species. The Javan warty pig (S. verrucosus ssp. verrucosus) is currently listed as 'Endangered' according to IUCN Red List criteria . S. verrucosus is endemic to Java where it is threatened by hunting, poisoning, and habitat destruction (Semiadi & Meijaard, 2006). An interview-based survey found a 53% population decline since 1982 and only highly fragmented populations persist (Semiadi & Meijaard, 2006, Figure 1b). The species was formerly also present on Madura, but deforestation is thought to have precipitated its extinction there (Blouch, 1988). The genetic resilience of the Javan warty pig may be further compromised due to introgressive hybridisation with the sympatric banded pig (Blouch & Groves, 1990;Drygala et al., 2020).
The small remnant volcanic island of Bawean (197 km 2 ) in the Java Sea is home to the fifth Indonesian suid, the Bawean warty pig (S. verrucosus ssp. blouchi). Groves and Grubb (2011) upgraded this suid to full species status with the name of Sus blouchi. The Bawean warty pig distinguishes itself from its sister species S. verrucosus by morphological (smaller skull and body size) and phenotypical (different by a recent bottleneck that reduced the effective population size to less than 20. The genomic assessment supports the single species status of S. blouchi, as was previously proposed based on morphometrics. The demographic history of S. verrucosus showed evidence of secondary contact with the sympatric banded pig (S. scrofa vittatus) that colonised Java 70 k years ago.
Main Conclusions: While the Javan and Bawean warty pigs have persisted throughout the Pleistocene climatic oscillations, contemporary pressures from human activities threaten their survival and immediate action should be taken to grant legal protection to both S. verrucosus and S. blouchi. This study highlighted the use of demographic history modelling using genomic data to identify evolutionary significant units and inform conservation.
Henceforth, we will follow the nomenclature proposed by Groves and Grubb (2011) to refer to the Bawean warty pig as S. blouchi. Owing to its small population size (172-377, Rademaker et al. (2016);234-467, Rode-Margono et al. (2020)) and its small geographic range, the Bawean warty pig was listed as 'Endangered' on the IUCN Red List (Rademaker, 2016;Rode-Margono et al., 2020). As its Javan relative, it is not legally protected by Indonesian law, and threatened by illegal logging and retaliatory hunting in response to crop-raiding (Nijman, 2003;Rode-Margono et al., 2016. Based on camera traps, S. blouchi is the only suid species on Bawean Rode-Margono et al., 2020). There are currently no Bawean warty pigs in captivity as part of a breeding programme to safeguard the wild population .
Here, the case of both the Javan and Bawean warty pigs was used to illustrate how long-term biogeographic processes and recent anthropogenic pressures can shape the demographic history of species, with significant implications for their conservation. The long-term persistence of the endangered S. verrucosus on Java and S. blouchi on Bawean depends on the timely implementation and compliance to conservation measures; yet significant knowledge gaps on the biogeographical context and the contemporary genetic diversity of these species may hinder full protection status.
Regarding the Bawean warty pig, this study aimed to conduct the first genomic assessment of S. blouchi, in relation to other suid species in the Indonesian archipelago. Specifically, we tested the hypothesis that S. blouchi diverged from S. verrucosus as recently as 10,000 years ago, when Bawean became isolated from mainland Java as suggested by Meijaard (2003). We thus expected genomic differentiation between S. blouchi and S. verrucosus to be much lower than differentiation with S. barbatus and S. celebensis that were estimated to have diverged ca. 2.8 million years ago (mya; Frantz et al., 2013). We hypothesised that the founder effect and ongoing anthropogenic pressures have led to depauperate genomic diversity and a small effective population size in the Bawean warty pig. As to S. verrucosus, it was hypothesised that ongoing hybridisation with the sympatric S. s. vittatus has resulted in significant levels of introgression from the wild boar into the warty pig gene pool. To gain a better understanding on the long-term co-occurrence of both suid species on Java, we compared alternative hypotheses of population divergence and connectivity using demographic model selection.
F I G U R E 1 (a) Present-day and past geography of the Indonesian Archipelago. During glacial maxima, vast areas of the Sunda shelf were exposed; the light grey area corresponds to the sea level at -116 m below present-day levels as described by Sathiamurthy and Voris (2006) with the Sunda river systems (Voris, 2000). (b) The distribution of the five Indonesian Sus species is shown as shaded areas: banded pig (Sus scrofa vittatus, purple), Javan warty pig (S. verrucosus, green), bearded pig (S. barbatus, teal), Bawean warty pig (S. blouchi, yellow), Sulawesi warty pig (S. celebensis, blue). S. verrucosus is endemic to Java where it only occurs in highly fragmented and small populations. S. blouchi is endemic to the small island of Bawean in the Java Sea. Distributional data for S. verrucosus, S. barbatus and S. celebensis were obtained from the IUCN Red List website; distribution of S. s. vittatus followed the description by Groves (1983). Tissue samples from trapped warty pigs were obtained using tissuesampling ear tags (CAISLEY International). No animals were killed for the purpose of this study. Biospecimen samples were stored in 96% ethanol and/or frozen at −20°C prior to laboratory analysis. Genomic DNA was extracted at the Indonesian Institute of Science (LIPI), Indonesia, using an ammonium acetate-based salting out procedure (Miller et al., 1988) or the spin column method using Qiagen DNeasy Blood and Tissue Kits (Qiagen, US).
Ten S. s. vittatus, 14 S. blouchi and 20 S. verrucosus were genotyped using the Porcine SNP60 v2 BeadChip (Illumina Inc.) following manufacturer's instructions. GenomeStudio 2.0 software (Illumina Inc.) was employed to call genotypes. The standard cluster file was used instead of a custom cluster file because the GenomeStudio clustering algorithms require ~100 samples to produce representative clustering positions (Illumina, 2014). Only autosomal SNPs, mapping to chromosomes 1-18 on the reference genome Sscrofa build 11.1, were retained for analysis. The resulting 50 K genotypes were merged with publicly available data of five S. verrucosus, six S. celebensis, nine S. barbatus and four common warthogs (Phacochoerus africanus; Yang et al., 2017) using PLINK 1.9 (Purcell et al., 2007).
Only loci present in both libraries were retained for analysis to reduce a potential batch effect. The presence of a batch effect was investigated by comparing clustering of S. verrucosus samples with respect to their library attribution.
Quality control filtering (call rate > 90% and missing genotypes <10%) of the merged dataset was carried out in PLINK 1.9. As the Porcine SNP60 V2 BeadChip was developed largely based on European domestic pig breeds, the application of the SNP panel to related species may lead to ascertainment bias (Barbato et al., 2020).
To reduce the effect of ascertainment bias, we applied a minor allele frequency filter (MAF >0. 01) and pruned SNPs with high linkage disequilibrium from the reduced SNP panel (using-indep-pairwise 2000 kb 10 0.5), as suggested by Barbato et al. (2020). These stringent filters reduced the available SNP panel from 51,265 to 7364 SNPs, hereafter referred to as the 50 K and 7 K SNP panels respectively.

| Molecular diversity and population structuring
Observed H o and expected heterozygosity H e and the inbreeding coefficient F IS were calculated for each partition with the dartR R package (v.2.0.3; Gruber et al., 2018) using both the 7 K and 50 K SNP panels. H e was adjusted for sample size. The degree of genetic divergence among species groups and partitions was estimated using Weir & Cockerham's (1984;hereafter Weir's F ST ) as implementing in the StAMPP package (v. 1.6.1; Pembleton et al., 2013) in R (v.3.6.0; R Core Team, 2019). Ninety-five percent confidence intervals (95% CI) were estimated based on 100 bootstraps across loci.
The Bayesian clustering approach as implemented in STRUCTURE (v. 2.3.4.;Pritchard et al., 2000) was employed with the 7 K SNP panel set to determine the most likely number of distinct genetic clusters K. The analysis was limited to the 7 K SNP panel set, because STRUCTURE assumes independence among loci (Pritchard et al., 2000). Both the no-admixture model with independent allele frequencies and the admixture model with correlated allele frequencies were tested because the most appropriate model, dependent on the amount of shared ancestry between S. verrucosus and S. blouchi, was not known a priori. An alternative ancestry prior of = 1 ∕ K as starting value was employed as recommended by Wang (2017).
Each estimation comprised an initial 70,000 iterations as burn-in, followed by an additional 200,000 iterations. Ten replicate estimations were conducted for values of K ranging from one to eight. The most likely number of K was inferred from Pr X| K , where X denotes the data, as described by Pritchard et al. (2000), and the ad hoc statistic Δ K developed by Evanno et al. (2005). For each K, the replicate with the highest Pr X| K was plotted with the pophelper R package An unrooted neighbour-joining (NJ) tree estimation was generated using the ape R package (Paradis & Schliep, 2019). For comparison, the PCA and NJ tree estimation were also conducted with the 50 K SNP panel set.

| Phylogenetic structure and divergence dating
To assess the support for splitting S. verrucosus and S. blouchi into separate taxonomic groups, we estimated bootstrap support for the population tree in TREEMIX (Pickrell & Pritchard, 2012) and employed Bayes Factor Delimitation in SNAPP (Leaché et al., 2014).
The topology of group splits was inferred from allele frequency variations of the 7 K SNP panel set with the TREEMIX algorithm (v. 1.13; Pickrell & Pritchard, 2012). P. africanus genotypes were used as outgroup to root the tree (Frantz et al., 2013). The TREEMIX input file was generated using the gl2treemix conversion function in dartR package. The maximum likelihood tree was estimated assuming a block size of 5 SNPs over three independent replicate runs.
Four species delimitation models ( Figure S1) were compared using Bayes Factor as model selection tool (Grummer et al., 2014). Leaché et al. (2014) modified the approach to work with SNP data by combining the Bayesian tree inference framework of SNAPP (v.1.5.2; Bryant et al., 2012), implemented as plug-in in BEAST2 (v. 2.6.7; Bouckaert et al., 2014), with methods for estimating marginal likelihoods. The 7 K SNP panel set was used because SNAPP assumes that all SNPs are unlinked (Bryant et al., 2012). The dataset was further limited to four samples from each group to reduce computational time. A detailed description of the models and parameter setting is provided in Supplementary Material S1. Bayes factors were The best-supported species tree was subsequently used to estimate divergence times. Due to ascertainment bias in the SNP data, branch lengths, which SNAPP estimates in coalescent units, cannot be readily converted to units of time. Instead, the time calibration approach described by Stange et al. (2018) was employed to calibrate the strict molecular clock model using divergence times from previous studies. The root of the species tree was constrained to (1) 4.2 mya (SD = 0.1) following Frantz et al. (2013) and (2)

| Demographic modelling
The effective population size N e of S. blouchi was estimated using the linkage disequilibrium method (Hill, 1981). NeEstimator (v.2.1;Do et al., 2014) was employed to estimate contemporary N e with both the 7 K and 50 K SNP datasets. The random mixing and monogamy models were both considered, given the mostly polygynous mating system in the Sus genus (Meijaard et al., 2011). Ninety-five per cent confidence intervals were estimated with the jackknife method (Do et al., 2014). A correction factor was applied based on linkage map length (~2000 cM, Tortereau et al., 2012) to reduce bias caused by the inclusion of potentially physically linked loci (Waples, 2006).
Trends in N e trajectories for S. blouchi were further investigated with SNeP (v.1.1; Barbato et al., 2015) using r 2 adjustments for limited sample size, a minimum of 100 items per bin, and Sved and Feldman (1973) as recombination rate modifier.
The occurrence of a recent population bottleneck for S. blouchi was tested using the programme BOTTLENECK (v.1.2.02; Cornuet & Luikart, 1996) with the infinite alleles model. The significance of heterozygosity excess compared to levels expected under mutationdrift equilibrium was tested using the standardised differences test (Cornuet & Luikart, 1996).

| Introgressive hybridisation
The advent of admixture between species was assessed using Patterson's D (Patterson et al., 2012) as implemented in Dsuite (Malinsky et al., 2021) to test for deviations from a strict bifurcating evolutionary history. The 7 K SNP data in VCF format were used as input file.
The extent of potential introgression from S. s. vittatus into S. verrucosus was investigated further using the two-layer hidden Markov model implemented in ELAI (V.1.01; Guan, 2014) using the 7 K and 50 K SNP panel set. The method assigns local ancestry with reference to ancestral populations that contributed to the current genomic composition in the admixed individuals. As no historical samples for S. verrucosus were available for this study, we used six S. verrucosus samples that originated from captivity (controlled breeding) and that were previously identified as pure Javan warty pigs (Drygala et al., 2020).   Table 2).
Both the no-admixture model with independent allele frequencies and the admixture model with correlated allele frequencies yielded comparable clustering solutions in STRUCTURE, except for higher admixture levels at K = 2 in S. barbatus and S. celebensis ( Figure S4-S7). The uppermost hierarchical level of structure, as identified by Evanno et al.'s (2005) Δ K method at K = 2 ( Figure S4 and S6), grouped S. verrucosus and S. blouchi into a separate cluster from S. s. vittatus, S. barbatus and S. celebensis ( Figure S5). However, log-likelihood estimates increased with higher Ks, levelling off at K = 4, thus supporting the presence of further clustering ( Figure S4 and S6). At K = 4, S. verrucosus and S. blouchi were partitioned into separate clusters, while S. barbatus and S. celebensis still formed part of the same cluster and S. s. vittatus accounted for the fourth cluster ( Figure S5). S. barbatus and S. celebensis only split into separate clusters at K = 5.
To test the hypothesis that the heterozygote deficit observed among S. verrucosus samples was caused by the Wahlund effect, the presence of population structuring was assessed using Bayesian clustering in STRUCTURE. The highest log-likelihood was observed at K = 3 ( Figure S8), which partitioned samples broadly into a cluster with the 10 individuals from the Yang et al. (2017) study, a western Java cluster (N = 5) and a central-eastern Java cluster (N = 14). Samples from captive individuals were part of the latter cluster.
The PCA (Figure 2a) and NJ tree estimation (Figure 2b) showed congruent patterns of hierarchical clustering and supported the differentiation into separate clusters between S. verrucosus and S. blouchi. In comparison, the 50 K SNP panel set yielded comparable clustering, but with S. s. vittatus showing stronger differentiation ( Figure S9).

| Concurrent divergence among insular Sus species
Tree topology inferred from TreeMix suggested high (100%) bootstrap support for the node split between S. verrucosus and S. blouchi The species delimitation model with a single group for S. verrucosus and S. blouchi was rejected in favour of a five-species model that separates S. blouchi as a distinct species. Bayes factor support for this model was decisive and other models that lumped species together ranked significantly lower ( Table 3).     Figure S1.  (Table S1, Figure 4b). This model was six times more likely to be the best model than the next-best model without the change in population size.

| Demographic histories marked by secondary contact and bottlenecks
For both JSFS, model residuals were overall low and goodnessof-fit testing suggested an adequate fit, despite some features of the empirical SFS not being fully captured by the simplified demographic histories.

| Low levels of genome-wide introgression from S. s. vittatus into S. verrucosus
Patterson's D estimates indicated significant deviations from expected frequencies of ABBA-BABA patterns within several groups ( clouded leopard (Neofelis diardi), Wilting et al., 2007), amphibians and reptiles (Inger & Voris, 2001). Geographical (e.g. river) and ecological barriers (e.g. climate, forest types) to dispersal have been proposed F I G U R E 3 Phylogenetic relationships among Sus species as inferred from SNP data in SNAPP. The tree root was time-calibrated (Stange et al., 2018) to estimate divergence times in millions of years (mya). Root calibration using (a) Zhang et al. (2022)  to explain the phylogeography of these species (Brown et al., 2013;Meijaard, 2003). In the case of Bawean, the East Sunda river system (Voris, 2000) or unsuitable xeric intervening habitat could have acted as a barrier to gene flow between Java and Bawean. Based on our results, a small number of warty pigs succeeded to disperse to Bawean during the Chibanian. The founder event was followed by a population expansion and divergence from the Javan warty pigs.
S. blouchi was characterised by low heterozygosity. This was in line with the expectation of lower genetic variation in island species (Frankham, 1997), which often experience genetic bottleneck effects during foundation, followed by loss or fixation of alleles caused by finite population sizes (Jaenike, 1973). Genetic variation is generally related to island size, with smaller islands, such as Bawean, supporting smaller populations and thus lower genetic diversity (Frankham, 1996). There was some evidence for a recent unidirectional admixture event from Java to Bawean, suggesting that the barrier between both islands was not absolute during some periods of lower seawater levels. Such gains from neighbouring populations F I G U R E 4 Highest ranking demographic model as inferred from δaδi for 2D site frequency spectrum (SFS) for (a) S. verrucosus (SVER) and S. blouchi (SBLO) and (b) for S. verrucosus (SVER) and S. scrofa vittatus (SSCV). Bottom right panel shows results of goodness-of-fit simulation with empirical estimate for this model shown as vertical line.
Founder event with exponential size growth, followed by a discrete unidirectional admixture event and population size changes "founder_admix_bottleneck_size" Divergence in isolation, continuous asymmetrical secondary contact with instantaneous size change "sec_contact_asym_mig_size"  Frantz et al. (2013), we assume that the violation of this assumption did not result in systematic bias. By presenting divergence times based on two very different Sus speciation estimates (Frantz et al., 2013;Zhang et al., 2022), we accounted for some uncertainty linked to the tree node calibration.

| Support for single species status of S. blouchi from molecular evidence
The single species status for the Bawean warty pig was previously proposed by Groves and Grubb (2011), with reference to Groves (1981 (Ramos et al., 2009). The application of this SNP panel to closely related species may have resulted in ascertainment bias (Helyar et al., 2011). Barbato et al. (2020)  If any ascertainment bias remained after filtering, the effect would be strongest in comparisons with S. s. vittatus (lowest divergence from the species used to ascertain the panel), whereas the relative divergence among other species would not be affected (Goedbloed et al., 2013;Iacolina et al., 2016). Due to incomplete lineage sorting, it is difficult to reconstruct the phylogeny of the Sus genus based on a limited number of markers (e.g. mitochondrial DNA haplotypes), as evidenced by the discrepancies in phylogenies inferred from morphological and mitochondrial DNA markers (Larson et al., 2007;Lucchini et al., 2005). Methods relying on thousands of markers provide a powerful tool to investigate species diversification, especially in cases of speciation with gene flow (Frantz et al., 2013(Frantz et al., , 2020.
Despite clear signals of introgression (Patterson's D) and secondary contact (δaδi), individual genome-wide ancestry levels were estimated to be very low (0%-1.8%). However, this was likely due to the sampling design of this study. Given our main research objectives, individuals with hybrid phenotypes and signs of admixture based on microsatellite markers (Drygala et al., 2020) were not included in this study. The reported level of banded pig ancestry in S. verrucosus is therefore not representative of the whole wild population, but only applies to 'pure'-looking individuals. The results presented here serve to confirm that Javan warty pigs identified as 'pure' using 14 microsatellite loci (Drygala et al., 2020) also lacked a signal of introgression using the 7 K SNP panel set.

| Conservation implications
The long-term persistence of species is strongly dependent on biogeographical processes (Whittaker et al., 2005). The Malay Archipelago is a prime example for species diversification and extinction as a result of periodic climate-driven connections and isolations of the islands (Lohman et al., 2011). Our demographic analyses highlighted that the divergence of S. blouchi can be dated to the Chibanian (Middle Pleistocene) and that S. blouchi has since undergone significant population size changes (i.e. founder effect, population expansion, bottleneck). However, the recent stark decline in the effective population size to ca. 30 is signalling an uncertain future for this species. The observed heterozygosity excess and low genetic diversity may further increase extinction risk due to a limited adaptive potential (Frankham, 2005).
Recent camera-trapping studies put the estimated population size at ca. 300 individuals Rode-Margono et al., 2020). The large difference in effective and census size is not unusual (Frankham, 1995) and likely the result of skewed male mating success if S. blouchi has a mating system similar to the wild boar (Delgado et al., 2008). Mounting human-wildlife conflicts on Bawean are posing a serious threat to the warty pigs as a consequence of increasing human population density. Unlike the protected Bawean deer, the warty pigs are considered a pest and legally hunted in retaliation for crop-raiding . Our results present a strong case that the Bawean warty pig is an evolutionary independent unit from S. verrucosus and requires prompt legal protection and conservation action.
Having survived through large-scale Pleistocene habitat alterations, S. blouchi may be facing its biggest challenge to date that requires a shift in people's perception of the species from a pest to a species worth protecting.
Concerning the Javan warty pigs, remediating population declines and fragmentation, as evidenced by genetic substructuring on Java, should be the main conservation priority. Secondary contact between the Javan warty pig and the banded pig, alongside reports of ongoing hybridisation, warrants further research into the extent of introgression in the wild population and the fitness of hybrid individuals. Both S. verrucosus and S. blouchi would benefit from legal protection to secure the long-term persistence of these endemic species in the face of anthropogenic pressures. In the meantime, well-coordinated ex situ breeding programmes could represent a vital lifeline for these species. Captive Breeding Centre, Indonesia, for co-operation and sample collection. Special thanks to Yudi Irawan from Javan Species Recovery (JaSpeR), Indonesia, for sampling. We also gratefully acknowledge the assistance of Yulianto during DNA extractions.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13689.

DATA AVA I L A B I L I T Y S TAT E M E N T
The 50k