Genetic diversity and demographic history of introduced sika deer on the Delmarva Peninsula

Abstract The introduction of non‐native species can have long‐term effects on native plant and animal communities. Introduced populations are occasionally not well understood and offer opportunities to evaluate changes in genetic structure through time and major population changes such as bottleneck and or founder events. Invasive species can often evolve rapidly in new and novel environments, which could be essential to their long‐term success. Sika deer are native to East Asia, and their introduction and establishment to the Delmarva Peninsula, USA, is poorly documented, but probably involved ≥1 founder and/or bottleneck events. We quantified neutral genetic diversity in the introduced population and compared genetic differentiation and diversity to the presumed source population from Yakushima Island, Japan, and a captive population of sika deer in Harrington, Delaware, USA. Based on the data from 10 microsatellite DNA loci, we observed reduced genetic variation attributable to founder events, support for historic hybridization events, and evidence that the population did originate from Yakushima Island stocks. Estimates of population structure through Bayesian clustering and demographic history derived from approximate Bayesian computation (ABC), were consistent with the hypothesized founder history of the introduced population in both timing and effective population size (approximately five effective breeding individuals, an estimated 36 generations ago). Our ABC results further supported a single introduction into the wild happening before sika deer spread throughout the Delmarva. We conclude that free‐ranging sika deer on Delmarva are descended from ca. five individuals introduced about 100 years ago from captive stocks of deer maintained in the United Kingdom. Free‐ranging sika deer on Delmarva have lost neutral diversity due to founder and bottleneck events, yet populations have expanded in recent decades and show no evidence of abnormalities associated with inbreeding. We suggest management practices including increasing harvest areas and specifically managing sika deer outside of Maryland.


| INTRODUC TI ON
The introduction of non-native species may have long-lasting ecological impacts to native communities (Sakai et al., 2001). While non-natives may provide some benefits, many also compete with or harm native species (Kalb, Bowman, & DeYoung, 2018). The factors that influence success of introduced populations are not well understood, and it is difficult to predict which introductions will become invasive (Sakai et al., 2001). The genetic diversity and structure of introduced populations may hold key insights into their ability to expand and adapt to novel environments (Sakai et al., 2001). Although non-natives often display reduced genetic diversity resulting from founder effect and genetic drift, some populations are able to avoid negative consequences of low genetic diversity or inbreeding (Lee, 2002). This observation has been termed the genetic paradox of invasive species (Estoup et al., 2016;Frankham, 2005). The paradox is not always valid, as some non-natives that would otherwise suffer from low diversity, avoid this and retain evolutionary potential, through multiple introductions or admixture of source stocks (Estoup et al., 2016;Kolbe et al., 2004). Furthermore, loss of diversity during bottleneck events may purge deleterious variation in some cases (Hedrick, 2001). Clearly, the success of non-natives depends in part on the introduction and founding history (Lee, 2002;Rius & Darling, 2014). Therefore, documentation of the genetic stocks involved and their demographic history should be useful in management (Cameron, Bayne, & Coltman, 2008;Sakai et al., 2001) as well as understanding the evolutionary response of introduced species to a new environment (Estoup et al., 2016;Lee, 2002).
Native to East Asia, sika deer (Cervus Nippon; Figure 1) have been introduced worldwide for sport hunting or as alternative livestock (e.g., the velvet antler industry). Sika deer are considered invasive in many areas where they have been introduced due to their effects on and interactions with native wildlife (Kalb et al., 2018;Senn & Pemberton, 2009).
In the Japanese Islands, sika deer were geographically isolated from mainland Asia by a series of vicariant events concurrent with cycles of glaciation (Riss-Würm) and changes in sea levels. These pervasive glacial events created small island areas with sika deer populations maintained by few individuals, which has resulted in isolated subspecies. Many of these populations have low neutral genetic diversity and are easily separated into clades based on mtDNA lineages (Tamate, 2009).
The largest free-ranging population of sika deer in North America occurs in the Mid-Atlantic US on the Delmarva Peninsula (a region that encompasses Delaware and the eastern coasts of Maryland and Virginia; hereafter Delmarva; Figure 2) where they have been shown to compete with native white-tailed deer (Odocoileus virginianus; Kalb et al., 2018). Records are sparse, but Japanese sika deer were apparently released to the Delmarva onto James Island in the Chesapeake Bay (2 km off the Delmarva mainland and part of Dorchester County, Maryland, USA) around 1916 as a single introduction of 4 or 5 individuals (Feldhamer & Demarais, 2009;Flyger, 1960;Kalb & Bowman, 2017;Kalb, Bowman, & Eyler, 2013).
The source population for free-roaming sika deer on the Delmarva originated in Japan, but spent several generations at Woburn Abbey, England (Feldhamer & Demarais, 2009;Kalb & Bowman, 2017) prior to their arrival in the United States. While in England, the deer were held in captivity with other stocks of unknown origin and allowed to intermingle (Banwell, 1999). In 1924, roughly 8 years after their introduction to the USA, the population of sika deer in Dorchester County, MD, still primarily on James Island, was divided (unknown quantities), and some deer were moved to Assateague Island in Worchester County, MD, USA (Figure 2; Flyger, 1960). A severe wildfire in 1957 reduced the population of sika deer on James Island by nearly half (Flyger, 1959;Flyger & Bowers, 1958).
Sika deer are now locally abundant on the southern portion of the Delmarva. Despite their founding from so few individuals, sika deer on Delmarva demonstrate remarkable vigor and have experienced near exponential growth over the last century (Davidson & Crow, 1983;Kalb & Bowman, 2017). However, the degree of genetic variation within wild sika deer that may allow future adaptation is not well understood. White-tailed deer have declined throughout portions of the Delmarva Peninsula resulting in areas where sika deer have become the primary cervid species. There is no evidence that it is possible for these two distinctly and evolutionarily different species (Polziehn & Strobeck, 2002) to hybridize with either sex-species configuration and is unlikely to happen in vivo due to differences in breeding season behavior and timing.
The free-ranging sika deer on the Delmarva are believed to be geographically isolated into two populations, one centering in and around Dorchester County MD, and the other around Assateague Island ( Figure 2). In Dorchester MD, Sika deer are a valued, yet controversial game species. Annually there are thousands of harvests throughout several counties in Eastern Maryland with a large hunter interest due to the exotic nature of the animal and unique experience of the sika deer rut behavior. However, neighboring states of Virginia and Delaware, USA (as well as some National Wildlife Areas), do not want established sika deer populations.
F I G U R E 1 Adult male sika deer from Dorchester County Maryland, USA. Image is from early spring just as males begin growing antlers. A large adult male may reach 50 kg dressed weight Introduced populations often experience dramatic changes in effective population size and geographic range, with implications for the maintenance of genetic diversity (Dlugosch & Parker, 2008;Hunter & Gibbs, 2007). Historically, translocations for most species in general were not well monitored or detailed in records.
Additionally, introductions by private citizens are often clandestine.
Although we only have evidence of a single introduction through multiple founder events, an investigation of genetic diversity and demographic history will provide evidence to better understand sika deer on the Delmarva.
The possible effects of low genetic variability could impact long-term management of this species. Low genetic diversity can have implications on the success of isolated populations, such as decreased fitness, and the limited ability of a population to adapt to changes in their environment (Baalsrud, 2011;Reed & Frankham, 2003). The introduced wild sika deer population offers a unique opportunity to investigate the effects of serial founder events on neutral genetic diversity (Sakai et al., 2001 The overall goal of this study was to quantify and compare the genetic diversity in the wild sika deer on the Delmarva and use the data to inform what we currently understand about their introduction to the Delmarva. Therefore, we compared the genetic diversity in free-ranging sika deer populations to the putative source population, as well as a population of captive sika deer on the Delmarva.
We hope that an analysis of the neutral genetic fingerprint of all sika deer on the Delmarva will allow us to better understand how sika deer spread across the Peninsula through identification of unique markers within geographic populations. Because of the different demographic histories we expected to find greater allelic diversity and more private alleles within the captive sika deer and Yakushima, Japanese samples. We also modeled the demographic history of the wild sika deer on the Delmarva to provide additional information about timing of and size of founder events. Finally, we discuss mechanisms that may have led to the success of populations from small founder events. in several locations, we accepted the help of hunters and managers to collect samples and data from their properties. We provided each collector detailed instructions on how to collect tissue and data.

| MATERIAL S AND ME THODS
We collected fecal samples from 12 captive sika deer (6 fecal samples from an all-male pen and 6 from an all-female pen) held in Harrington, DE, USA (Table 1). We collected only fecal samples that were fresh (e.g., pellets were wet, shiny, and covered in a thick mucosal layer) or where we witnessed defecation.
We placed all samples (fecal, muscle, and liver) into sterile 50-mL centrifuge tubes containing 25 ml of 95% ethanol (enough to cover the sample); we changed gloves between samples. For all samples, we recorded the collection date, harvest location, sex of deer, and species of deer on the 50-ml tube as well as on paper with pencil inside the tube. We stored all samples in a dark cabinet until we conducted DNA extraction. We extracted DNA from samples with commercial (Qiagen or Bioline) extraction kits for the appropriate sample type and according to the manufacturer's instructions for maximum quantity yield.
We obtained 14 DNA samples from the Yakushima Island sika deer population, Japan, the hypothesized source of the free-ranging sika deer on the Delmarva (Table 1), courtesy of Dr. H. Tamate and associates from Yamagata University. Samples were shipped as desalted pellets and were resuspended and diluted to a concentration of 100 ng/µl. We stored all DNA at −80°C until analysis.

| Identification of species
On Assateague Island and in Dorchester County, sika deer and white-tailed deer are sympatric and hunting seasons are concurrent.
To confirm that the samples collected were sika deer, we performed a restriction fragment length polymorphism (RFLP) assay for each sample (n = 109).
We identified a set of primers that amplifies a section of the mitochondrial D-loop, resulting in a fragment of ~464 base pairs (bp) in length for both sika deer (Wolf, Rentsch, & Hübner, 1999) and whitetailed deer. We identified differences in base-pair composition of the sequences using the National Center for Biotechnology Information (NCBI: GenBank) for both Yakushima sika deer and white-tailed deer. We slightly modified the light strand primer from Wolf et al. (1999), L14735, at two separate base-pair point changes: transitions at position 9 C to T, and position 16 T to C. Our modifications to Wolf et al. (1999)  The heavy strand primer (H15149) differed in two positions for sika deer and one position for white-tailed deer (Seabury et al., 2011;Wada et al., 2007), all at unique locations, and therefore, we did not modify it. We amplified this fragment using the polymerase chain reaction (PCR), with an initial 3-min denaturation at 96°C, followed by 45 cycles (96°C for 30 s, 60°C for 1 min) and final extension at 72°C for 3 hr (Wolf et al., 1999). The reaction mix included 5 µg of each primer (Invitrogen), 0.05 µg bovine serum albumin (BSA, Promega), 12.5 µl MyTaq Red Mix (Bioline), about 50 ng DNA template, and enough deionized (DI) water to bring the final volume up to 25 µl.
Based on the sequence comparisons (Seabury et al., 2011;Wada et al., 2007), we selected restriction enzyme HinF1, which cut fragments that were large enough to be viewed easily on an agarose gel. HinF1 cut G|AN(R = [A])TC in white-tailed deer, resulting in two fragments (198 and 265 bp), but did not cut the sika deer fragment ( Figure S1). Our restriction cocktail consisted of 5 units of HinF1, 0.025 µg BSA (Promega), 2 µl RE 10X Buffer (Promega: included with enzyme), and 5 µl of the PCR products, and was brought to a final volume of 10 µl with DI water. We digested the samples at 37°C for 3 hr and then held them at 4°C. We electrophoresed the products on a 2% agarose gel containing ethidium bromide (5 µl per 100 µl of agarose mix) and visualized results under UV light.

| Genetic diversity and population structure
We evaluated 16 microsatellite DNA loci developed for cattle (Bos taurus), sheep (Ovis aries), or reindeer (Rangifer tarandus) that successfully amplified in multiple species of cervids (Anderson et al., 2002;Bishop, Kappes, Keele, & Stone, 1994;McDevitt et al., 2009;Tamate et al., 2000;Wilson, Strobeck, Wu, & Coffin, 1997) and displayed high polymorphism and low genotyping error rates (Table S1). Due to low-quality results in preliminary analysis, we omitted six of the loci and used the remaining 10 for all samples.
We carried out amplifications as follows: 4 min denature at 96°C followed by 35 cycles of (30 s denature at 96°C, 60 s annealing (specific temperatures listed in Table S1), and 90 s extension at 72°C) and a final extension of 1 hr at 72°C before being held at 4°C. We created our PCR reaction mix using 10 µl PCR Master Mix (Promega: 400 µM each dNTP; 3 mM MgCl 2 ), 5 µM of both forward and reverse primers, 0.05 µg BSA (Promega), 50 ng template DNA, and enough DI water to bring the final volume to 20 µl per sample. We amplified all loci individually and combined them post-PCR for fragment analysis. We used primers fluorescently labeled at the 5′ end with Note: Counts of samples that were included from each group in each type of analysis are also listed.
TA B L E 1 Wild and captive sika deer sample counts and types from the Delmarva Peninsula, USA, and Yakushima, Japan Applied BioSystems G5 filter set (see Table S1 for dyes) and ran the samples on ABI 3730 DNA analyzer at the Delaware Biotechnical Institute. We scored alleles using GeneMapper Software 3.7 (Applied BioSystems). We only scored peaks in GeneMapper that were 100 relative fluorescence units (RFU's) or greater. Low allele peaks (<300 RFU's) were not common and we only accepted them when they were alleles that were already common in the population (following a relative threshold of calling: Whitlock, Hipperson, Mannarelli, Butlin, & Burke, 2008). We scored heterozygous alleles when the second peak was within 50% of the primary peak's RFU size. All RFU allelelike peaks that fell within the expected range for the loci but were less than 3X the background noise were scored as missing.
We indexed genetic variation among populations based on allelic richness and observed heterozygosity. We compared the observed levels of heterozygosity to expected values by testing for Hardy-Weinberg equilibrium using the program GENEPOP (Raymond & Rousset, 1995;Rousset, 2008). For all loci that had four alleles or fewer per population, we used the Fisher's exact test. For locus OarFCB193, we used a Markov chain method at 100 batches and 1,000 iterations per batch. We determined significance for both tests at the 0.05 level (Goodman et al., 2001). The number of observed alleles in a population is influenced by sample size. Since our sample sizes from each area were different, we quantified genetic diversity and estimated private alleles for each population using a rarefaction procedure in the program HP-Rare based on our smallest sample size (Kalinowski, 2004(Kalinowski, , 2005. We used an analysis of molecular variation (AMOVA) based on Wright's F-statistics (F ST , F IT , and F IS ) implemented in the program GenAlEx (Peakall & Smouse, 2006Wright, 1965).
Measurements of F-statistics were calculated across all 10 loci, and statistical significance was assessed based on 999 permutations among individuals or populations, as appropriate. We also estimated genetic diversity within individuals, within populations, between populations, and between geographical regions with GeneAlEx.
We estimated population structure among sika deer from Delmarva, Japan, and the captive deer using multivariate and Bayesian clustering analyses. We performed a principle coordinate analysis (PCoA) to compare genetic similarity of individuals within and among populations. We performed the PCoA on matrices of genetic distances between individuals based on a converted covariance matrix using the computer program GeneAlEx (Peakall & Smouse, 2006. We also estimated population structure with a Bayesian clustering analysis performed in the computer program STRUCTURE (Prichard, Stephens, & Donnelly, 2000). We used a burn-in of 100,000 Markov Chain Monte Carlo reps followed by 1,000,000 iterations of data collection. We estimated clusters (K) from 1 through 7, with 8 repetitions of each K. We determined the best-fit cluster solution using a modification of the Evanno method (Evanno, Regnaut, & Goudet, 2005) implemented in the computer program STRUCTURE HARVESTER (Earl & von-Holdt, 2012).

| Demographic history
We used approximate Bayesian computation (ABC) to make inferences about wild sika deer on the Delmarva including the potential for serial founder events, using the computer program DIYABC v 2.0.4 (Cornuet et al., 2008). Approximate Bayesian computation analysis can model complex population histories, including bottlenecks or founder events, changes in population sizes, and admixture.
We quantified support for demographic scenarios by generating simulated posterior probability models based on given demographic priors in a coalescent framework (Cornuet et al., 2008;Lawson Handley et al., 2011). We selected uniform probabilities for all scenarios on all parameters.
All priors provided in DIYABC scenarios were based on the best information available regarding historic and current status of wild sika deer on the Delmarva and Yakushima Island, Japan, including a proportional buffer (Table 2). We estimated current effective population sizes based on total population estimates for Dorchester, Assateague, and Yakushima, Japan; we assumed that Yakushima, Japan samples represent the allelic diversity of the founding stocks (Table 2). The DIYABC output includes estimates of all parameters for timing (in generations: T i ), duration of bottlenecks (in generations: db i ) and population sizes after bottlenecks or founder event (Nf i ).
For each ABC analysis we created one million simulated data sets per scenario, using the same number of loci and individuals as the original data set. We ran all mutation rate (µ) models with prior distributions for a generalized stepwise mutation model (GSM) for each locus (uniform mean mutation rate: 1.00E-004 min, 1.00E-3 max; uniform coefficient p: 1.00E-001 min, 3.00E-001 max; gamma individual mutation rate: 1.00E-005 min, 1.00E-002 max; and gamma individual coefficient p: 1.00E-002 min, 9.00E-001 max). We calculated the relative confidence in each set of scenarios via polychotomous logistic regression using the best 0.1 proportion of the data sets simulated. We calculated posterior distributions from this top 0.1 proportion of the data from the best scenario using linear regression of the logit-transformed results (Cornuet et al., 2008).
We considered a range of demographic scenarios in varying com-  Note: Effective population sizes were estimated as between 10% and 25% (Palstra & Fraser, 2012) of total population estimates when provided unless otherwise noted. Resulting posterior estimates shown as means, medians and 95% and 5% bounds. scenarios in groups of 2-4, changing single parameter events (e.g., did bottleneck happen at the same time period as the introduction or were there two separate events) and compared the best scenario from each for a final best case estimation of timing and magnitude of demographic events.
Our simplest scenario was wild Delmarva sika deer were founded directly from the Yakushima, Japan population: a population of sika deer that split into modern Yakushima, Japan samples and a second branch that initiated wild Delmarva sika deer and later split into iso- To verify the performance of the selected model, we used the model check option in DIYABC to estimate the goodness of fit of the simulated data sets with our original data (Cornuet et al., 2014).
Model check was based on 18 summary statistics across all three original populations including mean number of alleles, allele size, F ST , mean genetic diversity, shared allele distance, and (δµ) 2 , a measure of genetic distance between populations. We also estimated our confidence in the selected model using linear discriminant analysis of the summary statistics to provide a confidence interval to differentiate between sets of scenarios (Cornuet et al., 2014).

| RE SULTS
We know that each sample was from a unique individual because all wild Delmarva samples were collected from harvested sika deer, there are roughly 12,000 free roaming throughout the Peninsula (Kalb & Bowman, 2017). Based on microsatellite results, we were able to identify each of the captive sika deer as a unique individual.
All of the sika deer samples from Yakushima, Japan, were collected from unique individuals. Our estimated error rate (program GENEPOP) was 4.9% based on population-locus combinations with more than 3 alleles per locus.
Our genotypic data set was 93% complete across all 10 loci, with 48% of missing genotypes observed at locus BM1225.
Allelic diversity ranged from one allele (OarFCB304) to eight alleles (OarFCB193) between populations, with a mean of four alleles per population (  (Table 4). Mean null allele frequency was 0.13 across 9 loci per population.

| Population structure
The AMOVA showed all populations were differentiated from each other based on F ST , F IS, and F IT (Tables S2 and S3). Most of the genetic variation of these populations was measured (by F ST ) among populations (54%) or within individuals (45%). In our PCoA analysis, the top two axes explained 51% of the cumulative variation between populations ( Figure 3). Wild sika deer (Assateague and Dorchester) samples produced similar, intersecting plots. Sika deer from Yakushima, Japan, both overlapped slightly with wild sika deer samples but were more broadly dispersed across the two axes and only slightly overlapping each other (Figure 3). Samples from captive sika deer clus-

| Demographic history
We had three top model scenarios in our approximate Bayesian computation ( Figure 5)  (Table 2). Bottleneck durations were longer in the US, 5.7 generations (db) and 6.1 generations (db2), compared with the duration in the UK 2.1 generations (db3).
We estimated that (in reverse order from current date) 28.3 (Nf1) effective breeders survived the fire on James Island which oc- Sika deer on many Japanese Islands have lower genetic variation than observed in mainland Asian populations of cervids (Lü, Wei, Li, Yang, & Liu, 2006). The colonization of Japan by sika deer through the rise and fall of sea levels during the Riss-Würm (North American equivalent: Sangamon) interglacial and Würm (North American equivalent: Wisconsinan) glacial periods of the Pleistocene produced small, isolated populations of sika deer (Tamate & Tsuchiya, 1995). These small populations display limited variation within most populations, as well as genetic differentiation between isolated populations (Goodman et al., 2001;Nagata, Masuda, Kaneko, & Yoshida, 1998;Nagata et al., 1999;Tamate, 2009). These geographic separations translated to captive populations on the Delmarva that are both physically (substantially taller and heavier) and genetically very different from the wild populations implying that the captive animals were sourced from several different provinces (Kalb & Bowman, 2017).
We observed more alleles and a greater number of private alleles in the captive sika deer samples compared to both the Yakushima, Japan samples and wild Delmarva sika deer samples (Table 3). The high proportion of private alleles found between Yakushima, Japan, and captive sika deer samples is likely a result of the historic separation between mainland Asian sika deer and Japanese sika deer, and because the captive sika deer are a mixture of different sika deer stocks some of which were known to be from Manchurian (mainland) decent (Olson, Whittaker, & Rhodes, 2013).
We observed a decline in genetic variation in populations congruently with the timing and pattern of establishment. In an increasing order of genetic variation, wild sika deer on Assateague were founded from wild deer in Dorchester, founded from stocks in England and Ireland, which were derived from Japanese sika deer (Goodman et al., 2001;Nagata, Masuda, Kaji, et al., 1998;Senn, 2009). While different loci were used in their evaluation, Senn (2009) and Nagata, Masuda, Kaji, et al. (1998) also observed a pattern of low genetic variation suggesting a long-term, serial loss of diversity in sika deer through bottleneck and founder events.
While most of the alleles that we observed in wild Delmarva sika deer populations were found in the Yakushima Japan samples (12  (Bedford, 1949). Alternatively, the alleles may exist in Yakushima but were not present in this sample. The private alleles may have been lost from Yakushima but retained in Delmarva, which seems unlikely because more variation should be lost from the Delmarva deer during the multiple founding events and bottlenecks. It is also possible that the alleles are novel and derived from mutation after introduction. Finally, the low allelic diversity and proportion of shared alleles observed in both wild Delmarva sika deer populations and the source stock in Yakushima, Japan, is consistent with a single introduction and subsequent founder event.
Genetic variation we observed in Delmarva sika deer is similar to other populations of ungulates that have been founded from few individuals. For instance, introduced populations of elk in Pennsylvania and white-tailed deer in Finland display reduced variation relative to the source stocks (Kekkonen, Wikström, & Brommer, 2012;Williams, Serfass, Cogan, & Rhodes, 2002). Elk in Pennsylvania showed 7 of 10 loci were fixed or had been reduced to two alleles (Williams et al., 2002). White-tailed deer introduced to Finland maintained greater allelic richness and higher heterozygosity (5.36, 0.692) across 14 F I G U R E 5 Top historical models selected by DIYABC. Timelines on the right are not to scale. Change in colors within images represents a population split, or a population bottleneck. Time 0 represents the current (collection) date. Image c, (enlarged) was selected as the best of all models. The right side from time periods t3-t6 are equivalent to sika deer being in the United Kingdom, and time periods 0-t3 represent sika deer on the Delmarva Peninsula (wild deer only). The black color represents some ancestral sika deer lineage (labeled ancestral) that through time evolved into two separate lineages with some level of genetic differentiation, the gray color is an unknown population that was not sampled. The transition on the left wing from black to blue and again from blue to green represents a bottleneck or population reduction in the genetic diversity. Where green Japanese lineages meet the gray unknown lineage, there is genetic transfer creating the purple lineage of sika deer. The yellow lineage represents the midlineage population of sika deer that stocked the Delmarva Peninsula with a bottleneck at transition to the blue and red colors. Color schemes are similar within images a and b, but timings differ and represent different potential introduction scenarios. Image a, was a top model from possibilities that did not have any genetic introgression from other sika deer sources (presuming that there was a straight line introduction from either the Japanese Islands or United Kingdom. Image b, was the top model for only ancestoral introgression prior to sika deer populations leaving the Japanese Islands (it includes a gray "ghost" population) loci than observed in sika deer and elk; however, they were founded from deer that were more highly variable (Kekkonen et al., 2012).
The population in Finland also received additional genetic variation from a secondary introduction 14 years postfounding and continued to grow, providing additional chance for genetic drift (Brommer, Kekkonen, & Wikström, 2015;Kekkonen et al., 2012).

Similar to the translocation of bighorn sheep (Ovis canadensis),
subsequent establishments of sika deer in the wild populations of the Delmarva show progressive decline in neutral genetic variation (Olson et al., 2013). Wild sika deer of the Delmarva were founded from a single introduction of stocks with low diversity, followed by a lag in population growth and at least one bottleneck event.
Additionally, since the population of wild sika deer on the Delmarva remained small for several generations, there was a potential for the loss of neutral genetic variation due to genetic drift.
Reductions in genetic diversity due to inbreeding can lead to reduced sperm count, decreases in birth rates, decreases in juvenile survival, and increased susceptibility to disease (Lawson Handley et al., 2011;O'Brien et al., 1985;Sakai et al., 2001). However, wild Delmarva sika deer have proliferated and in some cases replaced native white-tailed deer. Wild sika deer were observed to have lower susceptibility to parasites (tics and other insects) and disease than native white-tailed deer (Davidson & Crow, 1983) in spite of their lack of genetic diversity. Sika deer on the Delmarva have also been successful reproducing, as emphasized in their near exponential population growth early in the introduction, and extensive, and increasing annual harvest throughout their expanding range (Kalb & Bowman, 2017).
Therefore, despite low neutral diversity, sika populations appear to be robust. In some cases, adaptive diversity can be maintained if the forces of selection outpace genetic drift, or if highly favorable genes are fixed via "allele surfing" (Hedrick, 2004;White, Perkins, Heckel, & Searle, 2013). It is possible that introduced sika have retained sufficient adaptive variation to be successful in their new environment.
Limited genetic diversity precludes fine-scale inferences on population substructure and assignment on the Delmarva sika deer. We calculated high values (Balloux & Lugon-Moulin, 2002;Wright, 1978) for F ST and F IT as a result of the complex nature of the introduced wild sika deer to the Delmarva. The calculated F IS values were not as high (Table   S3), but do suggest some degree of population structure within the sampling area due to social behavior or likely geography; such values may also be a result of allele amplification issues. The assignment into population clusters also confirms a lack of gene flow between captive sika deer and the wild sika deer on the Delmarva. The Yakushima Japanese samples were found to have mixed population assignment with some individuals more similar to captive sika deer samples and some individuals more similar to the wild sika deer samples.
While approximate Bayesian computation is sensitive to low sample sizes, the results strongly supported the recorded history of the introduction of wild sika deer to the Delmarva. The best-fit demographic scenario for our computational analysis involves a single introduction of sika deer forming the wild population. The rate of genetic admixture in our best-fit scenario suggests that there was little introgression of new genes in sika deer that were in the United Kingdom. The duration of bottleneck events were much longer in the US (  (Powerscourt, 1884). In the US, the range and population size of sika deer were most likely restricted from an invasive lag. We define invasive lag as a specific time frame between introduction and establishment of the invader. This could be years or decades before the population begins to experience rapid growth and range expansion consistent with invasive species. In the case of sika deer on the Delmarva, it appears that the lag period was a few generations and was probably increased in length due to the fire event and population collapse.
Generation estimates from ABC follow a similar growth rate estimated in bottleneck durations. The average years per generation were also less in the UK (2.7) than were estimated in the US (4.2) based on known dates (Table 2:  The continued growth and expansion of sika deer range, and increasing population size despite reduced genetic variability are an interesting study example of the genetic paradox of invasive species (Estoup et al., 2016;Frankham, 2005). The timing of introduction relative to sika deer exhibiting invasive behavior, however, was likely restricted by several large population bottlenecks and resulting in a major invasive lag. The population at large, throughout the Delmarva provides a semigeographically isolated and intensely managed study and research opportunity and we recommend its continued utilization for large scale, terrestrial invasive research.
We also encourage the Delaware Department of Natural Resources to design and implement management protocols specifically for sika deer, which should include liberal harvest regulations and ample administration of permits for agricultural damage (ideally to prevent sika deer spread through the state). Most importantly, Delaware officials will need to address the future of captive sika deer on the Delmarva, and if the potential risks associated with an escapee or released animal are warranted.
Future studies regarding sika deer genetics would benefit from a review of known restorations and introductions to help identify available markers that may be more informative at the population level. Studies within the Delmarva could focus on using a wider set of microsatellite loci and evaluating single-nucleotide polymorphisms, which would be informative in addressing if the sika deer of the Delmarva have evolved or maintained any adaptive diversity during their founding. We also suggest a more thorough sampling of sika deer on the Delmarva, especially in Wicomico, MD, and Sussex, DE to evaluate the spread of sika deer as individuals from the two wild populations, Dorchester and Assateague, begin to meet again.

ACK N OWLED G M ENTS
We thank the University of Delaware Department of Entomology

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
This paper is part of the doctoral research of David Kalb at the university of Delaware department of Entomology and Wildlife Ecology. Dr. Bowman, Dr. Delaney, and Dr. DeYoung all served on the dissertation committee and were contributors to all parts of the manuscripts success, from conception to completion.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sampling locations and microsatellite genotypes will be stored in