Human‐mediated introduction of introgressed deer across Wallace's line: Historical biogeography of Rusa unicolor and R. timorensis

Abstract In this study we compared the phylogeographic patterns of two Rusa species, Rusa unicolor and Rusa timorensis, in order to understand what drove and maintained differentiation between these two geographically and genetically close species and investigated the route of introduction of individuals to the islands outside of the Sunda Shelf. We analyzed full mitogenomes from 56 archival samples from the distribution areas of the two species and 18 microsatellite loci in a subset of 16 individuals to generate the phylogeographic patterns of both species. Bayesian inference with fossil calibration was used to estimate the age of each species and major divergence events. Our results indicated that the split between the two species took place during the Pleistocene, ~1.8 Mya, possibly driven by adaptations of R. timorensis to the drier climate found on Java compared to the other islands of Sundaland. Although both markers identified two well‐differentiated clades, there was a largely discrepant pattern between mitochondrial and nuclear markers. While nDNA separated the individuals into the two species, largely in agreement with their museum label, mtDNA revealed that all R. timorensis sampled to the east of the Sunda shelf carried haplotypes from R. unicolor and one Rusa unicolor from South Sumatra carried a R. timorensis haplotype. Our results show that hybridization occurred between these two sister species in Sundaland during the Late Pleistocene and resulted in human‐mediated introduction of hybrid descendants in all islands outside Sundaland.


| INTRODUCTION
Biogeographic barriers interrupt migration and reproduction among populations and thus are an important force responsible for driving and maintaining genetic differentiation, potentially leading to speciation. Sundaland, a Southeast Asian biodiversity hotspot, is bordered in the East by one of the best known faunal boundaries-the Wallace line (Bacon, Henderson, Mckenna, Milroy, & Simmons, 2013). This barrier is responsible for a sharp break between faunal compositions of Sunda and Wallacea. At the southern border, the Wallace line runs between the islands of Bali and Lombok, and at the northern edge, it divides fauna and flora of Borneo and the Philippines from that of Sulawesi ( Figure 1). Although some species and populations have naturally dispersed across this barrier (squirrel sps.; Mercer and Roth 2003), the presence of the majority of Sundaic mammal species occurring past the Wallace line into the eastern islands of Wallacea is associated with human transportation (Groves, 1983;Heinsohn, 2003;Veron et al., 2014).
Sundaland's dynamic geologic and climatic history, especially during the Plio-Pleistocene epochs, resulted in sea level changes that repeatedly exposed the continental shelf connecting the major islands of this archipelago (Voris, 2000). It is believed that these available land bridges would allow populations previously isolated on single islands to disperse, creating a large panmictic population within this whole system (Latinne et al., 2015;Demos et al., 2016). However, deep geographical barriers like the Wallace line would remain through even low sea level periods, thus creating patterns of genetic divergence between taxa on both sides of these barriers.
Here, we investigated the phylogeographic patterns of two Rusa species: the sambar, Rusa unicolor, and the Javan deer, Rusa timorensis.
While R. unicolor is widespread throughout South and Southeast Asia (from India and Sri Lanka, Southern China and most of Indochina to Borneo and Sumatra, the two largest Sunda Islands), R. timorensis has its native range on Java and Bali only (Figure 1). The presence of Javan deer on islands east of the Wallace line (e.g., Lesser Sunda Islands, Sulawesi, and the Moluccas) is described to be the result of prehistoric to historic human-mediated introductions during the Holocene. These human-mediated introductions were attributed, for example, to the Austronesian-speaking peoples' migrations, ca. 4,000 years ago, responsible as well for the introduction of other deer species (e.g., muntjacs), pigs, macaques, and civets; (Heinsohn, 2003;Groves & Grubb, 2011).
The aim of this study was to compare phylogeographic patterns, genetic diversity, and evolutionary history of these two related species in order to answer the following three questions: (1)

| Sampling and DNA extraction
We sampled 110 individuals labeled as Rusa unicolor (RUN) and Rusa timorensis (RTI) from European museums, aged between 180 and 61 years old. We collected either turbinal bones from the nasal cavity, skin, dry tissue from skeletons, and antler drills exclusively from individuals with known locality. All molecular work, including DNA extraction and sequencing library preparation, was conducted in a laboratory dedicated to work with archival samples to reduce the risk of contamination. DNA extraction followed the DNeasy Tissue and Blood kit protocol (Qiagen, Hilden, Germany), with overnight digestion of samples in Lysis buffer and Proteinase K at 56°C and a preelution incubation for 20 min at 37°C.

| Mitochondrial genome
All extractions, including negative controls, were built into individual sequencing libraries with single 8-nt indexes (Fortes & Paijmans, 2015), F I G U R E 1 Distribution map of both species and sampling location. Light gray indicates the distribution range of Rusa unicolor (a). Dark gray indicates the native distribution of Rusa timorensis, (b) and dashed dark gray areas indicate introduction range of R. timorensis (C). Filled circles and diamonds indicate R. unicolor and R. timorensis, respectively, and size is proportional to the number of samples. For detailed information about all samples see  (Maricic, Whitten, & Pääbo, 2010). Baits for hybridization were obtained by amplifying three overlapping mitochondrial fragments from one fresh tissue sample of R. unicolor (from the IZW archive) which were consequently prepared into capture baits (Maricic et al., 2010; primers and PCR conditions as described in Martins et al., 2017). After hybridization capture, libraries were amplified for no more than 18 cycles and sequenced again on the Illumina MiSeq platform. Mitogenomes obtained were deposited in Genbank (for accession numbers see Table A1).

| Microsatellite DNA
Microsatellite genotyping was achieved by amplification of 18 loci on all 56 samples for which mitochondrial DNA was obtained as described above (microsatellite loci and references in Table A2). All samples were amplified through PCR with the Type-it Microsatellite PCR kit (Qiagen, Hilden, Germany), with 1 μmol/L of each primer. Annealing temperatures followed a gradient from 63°C to 55°C in 2°C steps and final amplification occurred for 40 cycles at 55°C. Allele sizes were determined on an ABI3130xl Genetic Analyser using GeneScan ™ 500 ROX (both Thermo Fischer Scientific, Darmstadt, Germany) as internal size standard. Alleles were scored with the software GeneMapper v.4.0 (Applied Biosystems, Germany).

| Genetic diversity, phylogeography, and differentiation times
All mitochondrial sequences obtained were aligned using the auto setting as implemented in MAFFT v7.245 (Katoh & Standley, 2013).
The relationship among all haplotypes was reconstructed by a median-joining (MJ) network using the software NETWORK v. 4.6.1.4 (Bandelt, Forster, & Röhl, 1999). Haplotypes were generated by removing noninformative sites and positions with gaps or missing data.
Haplotypic and nucleotide diversities for the full dataset and for each species were assessed with the software DNASP v.5.10 (Librado & Rozas, 2009). We estimated genetic differentiation through F ST as implemented in ARLEQUIN v.3.5.12 (Excoffier, Laval, & Schneider, 2005). For this analysis, we created two datasets: (1) two populations corresponding to species as determined by the museum identification and (2) populations corresponding to major haplotype clades.
The best fitting substitution model for the full mitogenome data-

| Population structure analyses
Analyses of microsatellite data proceeded by removing all individuals with missing data at more than two loci. We estimated the probability for the presence of null alleles on our dataset with the software FREENA (Chapuis & Estoup, 2007

| Mitochondrial genome analyses
The final dataset consisted of 56 individuals for which a mitogenome of 16,064 bp was obtained. Of these, 23 samples were labeled in the museum collections as R. timorensis and 33 individuals were labeled as R. unicolor (Table A1). The origin of three specimens labeled as R. unicolor (RUN37 Moluccas, RUN39 Timor, and RUN61 Java) suggests that these specimens are actually Javan deer, R. timorensis. All other museum labels matched the geographical distribution ranges of the two species. The 56 deer shared 46 haplotypes with an overall haplotype diversity of H = 0.983 (SD 0.010) and a very low nucleotide diversity of π = 0.00941 (SD 0.00118). Within individuals labeled as RTI, we found 18 haplotypes with H = 0.972 (SD 0.026) and π = 0.0085 (SD 0.0024).
The full mtDNA haplotype network separated two major clades by a minimum of 168 mutational steps ( Figure 2). The smaller clade comprised six individuals from Java, Bali, Moluccas, and South Sumatra Sulawesi, and the Moluccas). Both ML and BI tree topologies were concordant with the overall pattern recovered by the haplotypic network, also revealing the existence of two well-differentiated clades ( Figure 3).
Generally, genetic diversity of haplotypes showed geographical structure only in some parts of the tree (subclades A, B, and C). Samples from the Moluccas, Sulawesi, and Timor (outside Sundaland) were present in more than one branch of the phylogenetic tree (subclade D).
The resulting node ages for the estimated age of Cervidae species were similar to those reported in other studies (Table 1)

| Microsatellite analyses
Of all archival samples for which mitogenomes were obtained, 16 individuals (~29%) could be successfully genotyped at 18 loci. These individuals were distributed relatively well across the clades of the mitogenome tree (Figure 3). Linkage disequilibria were found at 10% of all pairwise loci combinations, yet without any consistency, and percentage of null alleles was 0.18. Therefore, we retained all loci for further analyses. We detected significant deviations from HHWE, which indicated probable population structure within our dataset. Expected and observed heterozygosities at each locus ranged from 0.23 (locus Mu_4D) to 0.92 (locus Mu_1_51) and from 0 (locus Mu_4D) to 0.61 (locus Roe09), respectively (Table A2). Number of alleles varied among loci, with the highest number found at loci Mu_1_51 and NVHRT48 (16 alleles) and the lowest found at locus Mu_4D with only two alleles (Table A2).
According to the ∆K approach, the most likely number of genotypic nDNA clusters was K = 2 ( Figure 5). These two main clusters corresponded well to the two species, sambar (R. unicolor, green cluster) and Javan deer (R. timorensis, red cluster). Population differentiation (F ST ) was always significant between the two species, independent of the grouping method (museum assignments or STRUCTURE analyses, Table 2).

| DISCUSSION
The mitochondrial genomes and the nDNA loci of sambar and Javan deer investigated here revealed an intriguing and surprising pattern of genetic diversity and population differentiation between the two species. Although monophyly of R. timorensis and R. unicolor remain undisputed, our results point to a more complex history of hybridization between species and multiple human-mediated introductions outside the Sunda Shelf.
The presence of two divergent matrilineages clearly indicates molecular differentiation between two groups of Rusa deer, which we interpret as the historical cladogenesis of both R. timorensis and R. unicolor.
Our fossil-calibrated estimates are corroborated by recent studies (Bibi, 2013 (Sathiamurthy & Voris, 2006). This fact raises the obvious question of what then maintained differentiation F I G U R E 2 Haplotypic network for all 46 haplotypes shared among the two species. Circle size is in accordance with frequency and color represents sampling location. Small black circles represent median vectors. All branches represent one mutation step, except when indicated otherwise by numbers on branches. Two major clades were recovered and are indicated by the species names. Haplotypes are described in Table A1 H_39 Borneo and thus allowed differentiation between species based on evolved ecological adaptations (Leonard et al., 2015). The climate on Java is characterized by a West-East gradient, a transition from a slightly seasonal climate in the West to a strongly seasonal one in the East. Central and East Java are characterized by drier, cooler climate (climate-data.org), and the vegetation has more grass areas than on the surrounding islands (Heany, 1991;Mishra, Gaillard, Hertler, Moigne, & Simanjuntak, 2010). Therefore, it is likely that F I G U R E 3 Mitogenome maximum likelihood tree of both species. Colors on tips represent sampling location (as in Figure 2) and stars represent split events with bootstrap values/Bayesian posterior probabilities lower than 90/0.95 (but bigger than 50/0.5  Results (K = 2) 0.14 <.001 Estimations were performed both for the mtDNA and nDNA. Populations were generated by either museum ID (both for mtDNA and nDNA) and by assignment of individuals to one of the two major clades (mtDNA) or to one of the genotypic clusters (nDNA  (Figures 2 and 3).
Genetic distances among RTI hybrid haplotypes and to other RUN haplotypes were very low, indicating a recent, thus human-mediated introduction to these Wallacea islands. Quite recently, it had been suggested to split R. timorensis into seven subspecies according to their occurrence on islands within and outside of the Sunda shelf (Mattioli, 2011;Hedges, Duckworth, Timmins, Semiadi, & Dryden, 2015). However, our data do not support such a suggestion, as all samples from the introduction range shared haplotypes, indicating a lack of differentiation among individuals from these Wallacean islands. Furthermore, the few samples from Java and Timor that could be genotyped showed genetic similarity, again indicating the lack of differentiation.

| Phylogeography and taxonomy of Rusa unicolor
Sambar is currently subdivided into five subspecies: R. u. unicolor R U N 5 8 _ M y a n m a r R U N 5 6 _ S r i L a n k a R U N 4 3 _ I n d o c h in a R U N 3 _ S r i L a n k a R U N 1 2 _ B o r n e o R U N 1 9 _ B o r n e o R U N 2 0 _ M e n t a w a i R U N 6 6 _ S u m a t r a R U N 5 _ B o r n e o R U N 3 2 _ S u m a t r a expansion. The branching order of these three old subclades indicated colonization from northern Indochina southwards to Sri Lanka and to the Sunda Shelf, respectively. The "younger" individuals within subclade D from India, as well as from Sumatra and Borneo, appear then to be descendants from a second/third natural dispersion wave (possibly from Thailand) during glacial periods of the Pleistocene, when low sea level again exposed the shallow Sunda shelf connecting all major islands (Voris, 2000;Bird, Taylor, & Hunt, 2005). During glacial periods, climate was drier and cooler in tropical regions (Gorog, Sinaga, & Engstrom, 2004). However, species that retained a broad ecological niche such as Rusa unicolor would have been able to utilize the newly emerged habitats Such a scenario would likely cause the haplotypic distribution pattern we observed here. During glacial periods, Sundaland was also connected to Southeast Asian mainland, allowing secondary admixture between formerly separated populations, thus generating the patterns we observe between haplotypes from Thailand, India, and Sundaland.
In contrast, the distinct Sri Lankan clade B provides additional evidence for the recognition of Sri Lankan sambar as being distinct. This support comes both from morphological assessments (Groves & Grubb, 2011) and karyotype differences (2n = 56 in Sri Lankan sambar vs. 2n = 58 in Indian and 2n = 62 in Chinese and Malaysian sambar; Leslie, 2011). Sri Lankan populations are often more genetically related to the Western Ghats than to other Indian regions. Very recently, a 40 bp insertion was detected in the control region of the mitochondrial DNA in samples from the Western Ghats (Gupta, Kumar, Gaur, & Hussain, 2015), whose presence we, however, were unable to verify due to method limitations.
Thus, further studies on the mainland populations that could confirm the presence of old splits within R. unicolor are of urgency. Increased sampling and, especially, the inclusion of nuclear markers could confirm the extent of gene flow between these populations or, conversely, true genetic divergence between the subclades recovered here. If confirmed, these populations might represent subspecies of sambar occurring in highly disturbed regions of Mainland Southeast Asia.

| Introductions past the Wallace line
It is generally accepted that the presence of Rusa deer on islands beyond Sundaland (excluding Philippines) was the result of human interference (Long, 2003;Groves & Grubb, 2011;Hedges et al., 2015). However, until now, these individuals were assumed to have been pure R. timorensis, collected, and transported for venison and as game species by humans from the islands of Java and Bali during the Holocene (Heinsohn, 2003). While the nDNA data clearly separated the two species-being highly concordant with their description from the museum collections-the mitochondrial genomes point to a more surprising pattern of past Pleistocene hybridizations. Human-mediated introductions of these deer occurred ca. 3,000 years ago (Heinsohn, 2003), and, to our knowledge, these results could therefore constitute evidence for the earliest transport of introgressed deer.  Hedges et al., 2015) and attained through husbandry before (Leslie, 2011). Hybridizations with fertile offspring have also been reported to occur between other deer species, among others between sambar and red deer Cervus elaphus (Muir et al.,1997) and red deer and Sika Cervus nippon (Smith, Carden, Coad, Birkitt, & Pemberton, 2014).
This fact has potentially important conservation implications for the two Rusa species analyzed in this study. Despite being one of the most widespread deer species in southern Asia, R. unicolor is today no longer abundant throughout most of its native range . Likewise, R. timorensis is currently considered a pest species in areas where it has recently been introduced (e.g., Australia) but has, however, decreased largely in population numbers in native and historical introduction regions (Java, Lesser Sunda Islands, Sulawesi, and the Moluccas). Both Rusa species studied here are now considered vulnerable by the IUCN/Red List of Threatened Species (Hedges et al., 2015;Timmins et al., 2015). Therefore, genetic monitoring of individuals, both at mtDNA and particularly also nuclear genomes, is necessary to assess whether pure RTI and RUN individuals are being introduced (or reproductively assisted in their native ranges) and not introgressed individuals. Moreover, more intensive and extensive sampling of R. timorensis on their native range is necessary to discern whether pure RTI populations still remain in Java and Bali or whether they are composed in their majority by hybrid individuals.

| CONCLUSION
In addition to representing the first comprehensive phylogeographical study on R. unicolor and R. timorensis, this study revealed surprising evolutionary histories of these two sister species.
Answering the questions poised before, we hypothesized that while climate adaptations were likely responsible for maintaining species monophyly, Pleistocene climate changes were responsible for secondary contact and consequent hybridization between sambar and Javan deer. We recovered a pattern of (possibly reciprocal) introgressions between the two species, facilitated by the presence of land corridors during periods of low sea levels in T A B L E A 2 Microsatellite loci details and results