Demographic changes and loss of genetic diversity in two insular populations of bobcats (Lynx rufus)

Among felids worldwide, only 6 of 38 species have stable or increasing populations, and most felid species are threatened by anthropogenic influences, especially habitat loss and fragmentation. We documented changes in genetic diversity in an isolated, reintroduced population of bobcats on Cumberland Island (CUIS), Georgia, USA, compared to another bobcat population on Kiawah Island, South Carolina, USA, that was naturally established and experiences limited immigration from the mainland. The CUIS population declined from 32 reintroduced bobcats in 1989 to 10e24 individuals during 2012e2019, and observed heterozygosity declined from 0.742 to 0.634 (SD1⁄4 0.240). Observed heterozygosity of bobcats on Kiawah was 0.699 (SD1⁄4 0.153). We estimated that one bobcat immigrated to Kiawah Island every 5.3 years. We compared the predictions of a novel population viability analysis (PVA) to empirical estimates of abundance and genetic diversity on CUIS and used our PVA to identify management actions that are likely to support long-term viability. Mean heterozygosity from the PVA (0.588, SD1⁄4 0.065) was within 1 standard deviation of the empirical estimate. The PVA estimated the population would decline following population restoration due to loss of genetic diversity and inbreeding depression. Translocations of one female every four years would stabilize allele heterozygosity similar to the Kiawah Island population, but even translocations of two females every two years would not restore heterozygosity to founder levels. The PVA predicted no management action would result in a one in five probability of extinction within 50 years of reintroduction, but all translocation strategies nearly eliminated extinction risk through


Introduction
As habitats become fragmented, wildlife populations become isolated and small. These populations are at increased risk of extinction from demographic stochasticity, loss of genetic diversity, the accumulation of deleterious mutations, and inbreeding depression (Frankham 2005). Inbreeding depression is a decline in relative fitness resulting from increased inbreeding and reduced genetic diversity, impacting demographic traits such as birth rate and survival (Palomares et al., 2012). Inbreeding and loss of genetic diversity are unavoidable in small, closed populations and have been well documented in predators (Frankham 2009;Ralls et al., 1988). Moreover, population modeling indicates that median times to extinction are accelerated 30%e40% by inbreeding depression (O'Grady et al., 2006), and that loss of neutral genetic diversity occurs faster in smaller populations (Frankham 2009).
The introduction of new individuals into a population through translocation can reduce extinction risk and increase genetic diversity by adding new alleles and changing allele frequencies (Whiteley et al., 2015), reducing inbreeding depression, and decreasing demographic stochasticity. However, translocations also carry risks of potentially spreading infectious disease or parasites (Kock et al., 2010). Also, translocations could be counterproductive if they result in outbreeding depression (reduced fitness of offspring when individuals from genetically different populations are crossed; Edmands 2007). Recent research suggests, however, that outbreeding depression is less of a concern and that genetic rescue can be accomplished even when divergent populations are crossed (Kronenberger et al., 2016).
Genetic management of wild populations is poorly documented and rarely applied (Ralls et al., 2018;Seddon et al., 2007). Because conservation requires urgent action based on available information, there is little opportunity to reduce uncertainty before a conservation action is undertaken. Typically, a population viability analysis (PVA) is conducted to help inform management decisions but there is always uncertainty in any PVA because of the large number of demographic and environmental variables required to assess population fates (Lacy 2018). Inherent uncertainty associated with PVAs means that the accuracy of model predictions, as well as the effect of potential conservation actions, are unknown. Moreover, the lack of genetic rescue attempts means that limited empirical data are available to evaluate the predictive ability of PVAs (Ralls et al., 2018).
Among felids worldwide, only 6 of 38 species have stable or increasing populations, and most, for example, the Iberian lynx (Lynx pardinus) and Eurasian lynx (L. lynx), are threatened by anthropogenic influences, especially the loss and fragmentation of habitat (IUCN 2020)https://www.iucnredlist.org/search?taxonomies&equals;101738&amp;searchType&equals;species. Severe ly threatened, small, fragmented populations may require translocations among subpopulations to maintain the viability of the species. For example, remnant populations of the Florida panther (Puma concolor coryi) declined to approximately 30 individuals by the 1980s (US Fish and Wildlife Service 1987) and exhibited a high frequency of morphological and physiological abnormalities (Roelke et al., 1993). However, the translocation of eight female panthers from Texas into south Florida in the mid-1990s effected the genetic rescue of this population. Florida panthers have increased in abundance and range, genetic diversity and survivorship have increased, and the frequency of physiological and morphological correlates of inbreeding depression have declined (van de Kerk et al., 2019).
The bobcat (L. rufus; Fig. 1) is a mid-sized felid with polygynous breeding that exhibits solitary behavior, and females are the sole provider of parental care (Anderson 1987;Eisenberg 1986). The IUCN lists the bobcat's status as Least Concern because its populations are stable and have expanded to most of its historic range, extending from southern Canada, throughout the United States, to northern Mexico (IUCN 2020). In 1988e1989, 32 bobcats were reintroduced to Cumberland Island National Seashore, a barrier island along the Atlantic Coast of Georgia, USA (Fig. 2). The purpose of the reintroduction was to restore a native predator to the island, and post-reintroduction monitoring indicated that the reintroduction was successful (Diefenbach et al., 1993). However, there is no evidence that bobcats can immigrate naturally to Cumberland Island, so long-term viability of the population is unlikely without human intervention (Diefenbach 1992;Diefenbach et al., 1993).
Almost a decade after reintroduction bobcats continued to persist on the island, where they effected a trophic cascade (Diefenbach et al., 2006). Approximately two decades following reintroduction, the population consisted of about 14 individuals and had an effective population size of 5e8 individuals, with only one individual appearing inbred and overall genetic diversity high compared to other mammalian carnivores (Diefenbach et al., 2015). Demographic information about the reintroduced population was collected immediately post-reintroduction (Diefenbach 1992), but whether a PVA based on this information would provide accurate predictions of longer-term population dynamics was not known.
Uncertainty exists over what type of conservation actions are necessary to ensure the viability of the Cumberland Island bobcat population in the long term. Although PVA models can predict loss of genetic diversity, empirical estimates could validate the accuracy of such models. Also, because isolated populations will always experience genetic loss, the maximum level of genetic diversity that should be expected given population size and population dynamics is unknown. Finally, since supplemental translocations will be required to sustain the population, information is needed to inform decisions about how many bobcats, of what sex, and how often translocations should occur.
Due to its isolation and location on the edge of the species' range, the restored bobcat population on Cumberland Island can be viewed as a proxy for a small population of an endangered species. Our first objective was to use data from the Cumberland Island bobcat population to develop a PVA based on data collected immediately after the reintroduction and to evaluate the predictions of the PVA regarding abundance and genetic diversity after 30 years. We extracted DNA from scats collected during 2012e2019 and compared empirical estimates of genetic diversity and abundance to the PVA. We also estimated the loss in genetic diversity empirically using blood samples collected from bobcats released on the island. Our second objective was to estimate expected genetic diversity for an island population by comparing the genetic diversity of the current Cumberland Island bobcat population to one on Kiawah Island, South Carolina that was naturally established and experiences limited migration from the mainland (Fig. 2). Our third objective was to use the results from these analyses to identify supplemental translocation strategies that are likely to ensure long-term viability of the bobcat population on Cumberland Island and maintain genetic diversity similar to that observed on Kiawah Island.

Study areas
Cumberland Island and Little Cumberland Island, hereafter Cumberland Island (CUIS), are barrier islands located along the Atlantic Coast, 0.5 km north of the GeorgiaeFlorida border (Fig. 2). These two islands are separated from mainland Georgia by 2e4 km of salt marsh and open water. Little Cumberland Island lies to the north, separated from the main island by salt marsh and a tidal creek. The total area of upland habitats on both islands is 6935 ha. Both islands comprise Cumberland Island National Seashore but only the larger island is administered by the National Park Service (NPS). No roads connect the islands to the mainland, and development on CUIS is limited to a few private inholdings, historical structures maintained by the NPS, and infrastructure to support park staff and visitors.
Kiawah Island is a 3200 ha barrier island located along the Atlantic Coast of South Carolina (Fig. 2). As of the 2010 census, there were 3484 housing units and 1626 residents (https://factfinder.census.gov accessed 3 February 2020). Although the island was developed with roads, hotels, houses, and golf courses, there are regulations to minimize removal of trees and other natural vegetation. The island is separated from the mainland by salt marsh and the Kiawah and Stono rivers. Bobcats that immigrate to the island must traverse either a bridge for motor vehicles or at least 50 m of open water at low tide and !550 m of salt marsh.

Sample collection
We collected whole blood samples from all bobcats at the time of reintroduction to CUIS, and all bobcats subsequently captured on the island from 1989 to 1991. Bobcats were captured in accordance with the American Society of Mammologists' guidelines for the use of wild mammals in research (Sikes and The Animal Care and Use Committee of the American Society of Mammalogists, 2016). Unfortunately, some samples were lost to storage equipment failure, but we obtained DNA from 31 bobcats reintroduced to the island and the first generation of kittens born on the island (hereafter referred to as CUIS-Founders). These blood samples were stored at À20 C to À80 C until immediately before DNA extraction.
We collected scat samples (n ¼ 236) from the current population of bobcats on Cumberland and Little Cumberland Islands (hereafter CUIS-Current) in 2012, 2016, 2018, and 2019. We walked roads, hiking trails, and the interdune meadow in late December and early January, except in 2019 in February, to search for scat. We recorded the global positioning system (GPS) coordinates of the path walked and the coordinates of each scat, as described previously (Diefenbach et al., 2015). We handled scats with disposable gloves to prevent cross-contamination of DNA and stored each scat at room temperature in a container with silica gel desiccant. We processed scat samples for DNA extraction within 1 weeke3 months after collection.
We collected 3 mm ear punch skin samples (n ¼ 36; 19 males, 17 females) from bobcats on Kiawah Island (hereafter Kiawah) from 2007 to 2015 during GPS collaring procedures (SCDNR Permit #SC-05-2020, Roberts 2007). Each sample was stored in tris buffer (100 mM Tris, pH 8.0; 100 mM EDTA, pH 8.0; 10 mM NaCl; 1% sodium dodecyl sulfate) at room temperature away from light. The project protocol for bobcat capture and handling was approved by the University of Georgia Institutional Animal Care and Use Committee (IACUC #A2002-10113-M1) and was in accordance with the American Society of Mammologists' guidelines for the use of wild mammals in research (Sikes and The Animal Care and Use Committee of the American Society of Mammalogists, 2016). Ten tissue samples from mainland South Carolina (hereafter SC-Mainland) were donated by fur trappers in 2015 and 2016. These bobcats were trapped in either Marlboro or Newberry counties, SC (Fig. 2). Tissue samples were stored in 70e90% ethanol at room temperature.

DNA extraction
Blood and tissue samples were extracted using the DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA, US), according to the manufacturer's instructions, and were eluted into 200 ml 1Â Tris-EDTA. We extracted DNA from scats in a separate laboratory dedicated to low-quality DNA samples, using either the Qiagen QiaAMP Power Fecal DNA Extraction Kit or the Qiagen DNeasy PowerSoil DNA Extraction Kit (Qiagen, Valencia, CA, USA). Modifications to the manufacturer's protocol are described in Supplementary Information. Negative controls were included to monitor for contamination. Extracted DNA samples were stored at À20 C.

Sexing
We determined the sex associated with each scat sample by amplifying the amelogenin region of the sex chromosomes. This AMELX/Y polymerase chain reaction (PCR) amplifies homologous fragments from the X and Y chromosomes, which differ in size by~20 bp (Pilgrim et al., 2005) and can be visualized by electrophoresis on a 1.5% agarose gel. Cycling conditions are described in Supplementary Information. We recorded sex after a minimum of three independent amplifications produced the same results.
We initially screened low-quality scat DNA samples for suitability for genotyping analysis by PCR amplification in a 4microsatellite panel (MPXB1 ; Table S1). We re-amplified samples that amplified successfully at two or more loci with a second 4-locus panel (MPXB2 , Table S1). We eliminated samples that amplified at four or fewer loci after both rounds of PCR (n ¼ 24). Those that amplified successfully at a minimum of five loci qualified for full analysis at all 11 or 13 loci (n ¼ 236). Multiplexes and cycling conditions are described in Supplementary Information (Table S1). We amplified high quality DNA samples from tissue and blood in a different laboratory and according to a different protocol than scat DNA samples, also described in Supplementary Information (Table S2). For both sample types, we included one negative control (deionized water) on each PCR plate to monitor for contamination, as well as two previously genotyped samples to ensure that genotype calls were reproducible across electrophoresis runs. Amplicons were visualized using an Applied Biosystems genetic analyzer (model 204 3730 XL; Waltham, MA, USA) at the Penn State Genomics Core Facility (University Park, PA, USA). We amplified each DNA sample a minimum of three times at each locus. However, up to 12 attempts were made per locus to amplify a lowquality scat sample before that locus was deemed to have failed for that sample.
Allele sizes were called using GeneMarker software (SoftGenetics, State College, PA, USA), based on a known DNA size standard (GeneScan 500 LIZ Dye Size Standard; ThermoFisher Scientific, Waltham, MA, USA). Genotypes were accepted at any individual locus only if that genotype was observed in at least two independent PCRs. Putative homozygous genotypes were accepted only if a single allele was observed in at least three independent PCRs. Thereafter, a consensus multilocus genotype was generated for each individual.

Population genetics
We identified individuals and potential duplicate genotypes from the scat samples based on multilocus genotypes using GenAlEx v. 6 (Peakall and Smouse 2012). We considered genotypes that were identical at !11 loci as recaptures of the same individual and eliminated duplicate genotypes from further analysis. We evaluated 236 samples at either 11 or 13 microsatellites and retained 90% (n ¼ 212) of the samples, of which 61% were genotyped at all 13 markers. The probability of identity (P (ID) ) and the probability of identity for siblings (P (ID)sib ) were also calculated in GenAlEx. The former is the probability that two unrelated individuals drawn at random from a population will have the same multilocus genotype, and the latter is the probability of observing identical multilocus genotypes in siblings (Waits et al., 2001). Although we collected scat samples at four different time points (2012,2016,2018 and 2019), we considered them as a single group (CUIS-Current) for the population-level statistical analyses because bobcats detected in 2012 could still be alive in subsequent collection seasons.
To examine data quality, we tested for the presence of null alleles, linkage disequilibrium, and deviation from Hardy-Weinberg expectations due to either homozygote or heterozygote excess using Micro Checker (van Oosterhout et al., 2004) and Genepop v. 4.7 (Rousset 2008). We quantified genetic diversity within each population by estimating the average number of alleles per locus adjusted for sample size (AR, allelic richness) in FSTAT (Goudet 2003), as well as mean observed heterozygosity (H O ) and mean expected heterozygosity (H E ) based on Hardy-Weinberg assumptions in Arlequin v. 3.5 (Excoffier and Lischer 2010). We compared these measures of genetic diversity among populations with a Kruskal-Wallis test, implemented in SPSS Statistics v. 26 (IBM Corporation). We also used Arlequin to calculate the modified Garza-Williamson M statistic (Garza and Williamson 2001) to test for evidence of a recent bottleneck. This statistic calculates the ratio of the mean number of alleles to the total range in allele size under the assumption that the number of alleles declines more rapidly than the range in allele sizes after a reduction in population size. Values of M close to 1 are associated with stable populations, whereas M < 0.68 is indicative of a recent bottleneck (Garza and Williamson 2001).

Relatedness and inbreeding
We estimated pairwise relatedness (coefficient of relatedness, r) for each pair of individuals within a population using both the maximum-likelihood method employed by ML-Relate (Kalinowski et al., 2006) and the triadic likelihood estimator in COANCESTRY (Wang 2011), which accounts for inbreeding (Wang 2007). In both cases, we implemented 10,000 randomizations to estimate 95% confidence intervals. Individuals that are unrelated have values of r close to zero. Parents and their offspring, as well as full siblings, share on average 50% of their genes (r ¼ 0.5), therefore values of r > 0.5 indicate higher than expected levels of relatedness, suggestive of inbreeding. We compared relatedness across all four populations with an independent samples Kruskal-Wallis test in SPSS Statistics v. 26 (IBM Corporation). We also used ML-Relate to calculate the likelihood of four possible relationships (parent-offspring, full-siblings, half-siblings, or unrelated) for each pair of individuals. Although both parent-offspring and full siblings have an average pairwise relatedness coefficient r ¼ 0.5, ML-Relate also estimates k-coefficients, which represent the probabilities that two individuals share alleles at a locus that are identical by descent, and which can be used to distinguish between different genealogical relationships (Kalinowski et al., 2006). We performed these relationship analyses only for the insular populations, CUIS-Current and Kiawah, where individuals are likely to be related. The CUIS-Founder and SC-Mainland individuals were collected over sufficiently large geographic distances in mainland Georgia and South Carolina (Fig. 2) that they are all assumed to be unrelated. We estimated population-level inbreeding, F IS , in FSTAT (Goudet 2003). This is the proportion of genetic variance in a subpopulation that is found within an individual. Values of F IS close to zero are expected under random mating, values close to one suggest inbreeding, and negative values indicate an excess of heterozygotes (Excoffier 2007). The latter is expected in an outbred population but can also be indicative of a recent population bottleneck (Cornuet and Luikart 1996).

Effective population size
As a comparison to our abundance estimates (described below), we estimated effective population size (N e ) to provide an indication of the effective number of bobcats that have contributed to the current gene pool. We used NeEstimator 2.01 (Do et al., 2014) to estimate current N e on CUIS and Kiawah using three single-sampler estimators (linkage disequilibrium, heterozygote excess, and molecular coancestry methods) for both islands, and the Nei/Tajima temporal method for CUIS. For the latter, we used the CUIS-Founder population as the first time point (time zero), and CUIS-Current as the second time point, 6 generations later (assuming a generation time of 4 years, as estimated in Vortex (version 10.0.7.0, Lacy and Pollack 2014), described in PVA section below). We were unable to use the temporal method for Kiawah as our samples were collected over a shorter timeframe that represented only one or two generations.
We estimated the historical demography of the CUIS island system using IMa2 v.8.27.12 (Hey and Nielsen 2004). A strict isolation model was implemented for the CUIS system, with the analysis estimating three parameters of interest: q ð¼ 4N e mÞ for CUIS-Current and for CUIS-Founders, and the time of splitting t ð ¼ tmÞ, where m is the mutation rate and t is the number of generations since divergence. Based on preliminary results, priors for this analysis were set to the uniform distributions (0,22) for q founder , (0,0.5) for q current , and (0,0.15) for t. Five independent runs were each allowed to burn-in for approximately 100,000 steps, and then run for an additional~30 million steps. Each run consisted of 50 heated chains with heating parameters h a ¼ 0.9 and h b ¼ 0.8. Substitution rates for most loci were estimated in IMa2 relative to an assumed rate of 10 À3 mutations per locus per year for locus FCA132, which was one of the more variable loci. This assumed mutation rate is of the same order of magnitude as that empirically derived from humans (Brinkmann et al., 1998). Coalescent-scaled parameters (q and t) were converted to natural parameters (N e and t, respectively) as described using the inheritance scalar 2.1797 estimated in IMa2. After verifying that all runs converged on similar posterior distributions, we used IMa2 to calculate the joint posterior densities for each parameter based on the 1,500,000 genealogies resulting from all five separate runs. We also estimated the historical demography of the Kiawah island system using IMa2. While we recognize that there is likely no direct migration between Kiawah and SC-Mainland as the latter samples were collected in the inland piedmont plateau region of SC, previous studies have shown minimal genetic differentiation among bobcats throughout the southeastern United States (Reding et al., 2012). Therefore, it is likely that the SC-Mainland samples are representative of coastal bobcats that do experience migration with Kiawah. An isolation-with-migration model was implemented for the Kiawah island system, with the analysis estimating six parameters of interest: q ð¼ 4N e mÞ for the current Kiawah population, for the current SC-Mainland population, and for their most recent common ancestor, directional migration rates m ¼ Mi m , where M i is the rate of migration into population i) between Kiawah and SC-Mainland, and the time of splitting t ð ¼ tmÞ. Based on preliminary results, priors for this analysis were set to the uniform distributions (0,2) for q Kiawah , (0,25) for q SouthCarolina , (0,50) for q MRCA , (0,50) for m, and (0,1.5) for t. Five independent runs were each allowed to burn-in for approximately 125,000 steps, and then run for an additional~30 million steps. Each run consisted of 50 heated chains with heating parameters h a ¼ 0.9 and h b ¼ 0.8. Substitution rates for most loci were estimated in IMa2 relative to an assumed rate of 10 À3 mutations per locus per year (Brinkmann et al., 1998) for locus FCA132. Coalescent-scaled parameters (q and t) were converted to natural parameters ( b N e and t, respectively) as described using the inheritance scalar 1.4269 estimated in IMa2. Migration was calculated as the effective number of migrants per generation 2N e M ¼ qm 2 After verifying that all runs converged on similar posterior distributions, we used IMa2 to calculate the joint posterior densities for each parameter based on the 1,500,000 genealogies resulting from all five separate runs.

Abundance and survival estimates
We estimated the abundance (N) of the present-day bobcat population on Cumberland Island (2012e2019) using spatially explicit captureerecapture population estimation, secr, version 4.1.0 (Efford et al., 2009) implemented in R (R-Core-Team 2013) using the sampling methods described in Diefenbach et al. (2015). We estimated abundance by sex and year. We analyzed data for males and females separately, and for each sex compared two models: 1) allowed density to vary each year but constrained g0 and s to be the same across years; and 2) allowed density, g0, and s to vary across years. We estimated total abundance for each year by summing abundance by sex and summing their respective variances to estimate variance of total abundance. We used the abundance estimates from the model with the lower Akaike's Information Criterion adjusted for sample size (AIC c ) value (Burnham and Anderson 2002). We estimated the F:M sex ratio for each year and estimated the SE using a first-order Taylor series approximation.
To estimate annual survival, we used a Cormack-Jolly-Seber (CJS) model that estimated the probability of capture and apparent survival (probability of surviving and remaining on the study area; Seber 1982) based on whether we collected a scat from an individual in 2012, 2016, 2018, and 2019. We assumed emigration did not occur, so the estimate of apparent survival was equivalent to true survival. We developed the CJS model in a Bayesian framework using the program JAGS (Plummer 2017) and Package "R2jags" version 1.0e6 (Yu-Sung Su & Masanao Yajima) in program R (R-Core-Team 2013). Because we did not sample the population every year during 2012e2019 we could estimate only a constant survival rate for the 4-year interval 2012e2016 and 2-year interval for 2016e2018. Also, we did not know the age of individuals, so we pooled all individuals for analysis and obtained a single estimate for survival for each year. There is little evidence for sex-specific differences in survival rates, although survival rates for adults are greater than juveniles (Anderson and Lovallo 2003). We assumed capture probability was not different among years because sampling effort was similar. We ran 100,000 iterations of three chains of the Markov Chain Monte Carlo search with no thinning and discarded the first 50,000 iterations. We inspected trace plots for evidence the chains converged and ensured that Brooks-Gelmin-Rubin statistics were <1.10.

Population viability analysis
We created six models of the viability of bobcats on CUIS using program Vortex version 10.0.7.0 (Lacy and Pollack 2014) that addressed demographic stochasticity, effects of inbreeding depression, and four different supplementation strategies. The baseline model (1) was a stochastic model solely using demographic parameters estimated from monitoring the population following the initial reintroduction (Diefenbach 1992). Post-reintroduction monitoring provided data on initial population size, adult survival, density-dependent reproductive rates, juvenile survival, age at first reproduction, and lifespan (complete parameterization of the Vortex model available in online data repository). We randomly selected from the allelic diversity of 31 individuals from CUIS-Founders to initialize the genetic diversity of the population in simulations.
The second model (2) was the same as the baseline demographic model except that inbreeding depression occurred with lethal equivalents ¼ 3.9 and 50% due to recessive alleles. We used 3.9 lethal equivalents because this is the number of diploid lethal equivalents estimated to affect fecundity (O'Grady et al., 2006). The remaining simulations were the same as the second model except that we implemented different strategies for supplementation: 3) add 1 female every 4 years beginning in year 30 (i.e., beginning in 2020); 4) add 1 adult female every 4 years beginning in year 30 and if the population is > 15 then also remove a juvenile female from CUIS; 5) add 2 adult females every 2 years and if the population is > 15 then also remove 2 juvenile females from CUIS; and 6) add 1 adult male every 4 years beginning in year 30. For each scenario we repeated the simulation 999 times for 100 time steps (years) and monitored abundance, genetic diversity, and extinction probability for each scenario.

Genetic diversity
The probability of identity (P (ID) ) and the probability of misidentifying siblings as the same individual (P (ID)sib ), based on 13 loci, were extremely low, both for mainland (CUIS-Founders and SC-Mainland) and island populations (CUIS-Current and Kiawah; Table 1). One locus (FCA391) showed evidence of null alleles and possible deviation from HWE (p ¼ 0.007) due to homozygote excess, but only in the CUIS-Founders population. A second locus (FCA096) deviated from HWE (p ¼ 0.001) due to heterozygote excess in the Kiawah population. However, neither locus showed evidence of null alleles or deviation from HWE in any other population, and therefore they were not excluded from the statistical analyses.
We identified 37 distinct bobcats from the 236 scats that comprised the original CUIS-Current sample. Each genotype was observed in an average of 6.378 ± 5.736 scats. Nine unique individuals were identified in 2012 (4 males, 5 females), 16 in 2016 (9 males, 7 females), 21 in 2018 (7 males, 14 females), and 22 in 2019 (9 males, 13 females), with multiple individuals recaptured across years. As expected, allelic richness and expected heterozygosity were higher (Kruskal-Wallis p < 0.019) in mainland populations (CUIS-Founders and SC-Mainland) than island populations (CUIS-Current and Kiawah; Table 1). Although the pattern was consistent in that observed heterozygosity was also lower in insular populations than mainland populations, statistically we did not detect a difference in H O among the four populations (p ¼ 0.538). CUIS-Current had the lowest genetic diversity of the four populations (AR ¼ 3.814 ± 1.264, H O ¼ 0.634 ± 0.240, H E ¼ 0.591 ± 0.213; Table 1). SC-Mainland had the highest allelic richness (AR ¼ 5.587 ± 1.327), H O (0.774 ± 0.152), and H E (0.748 ± 0.110). Genetic diversity on Kiawah was higher than CUIS-Current, but lower than the mainland populations (Table 1). The inbreeding coefficient, F IS , was weakly negative in all four populations, but there was an excess of heterozygotes only in CUIS-Current (F IS ¼ À0.073, p ¼ 0.007) and Kiawah (F IS ¼ À0.086, p ¼ 0.003; Table 1). The Garza-Williamson modified index, M, was high in both mainland populations (M ! 0.746) but low in both island populations (M 0.600), indicating the latter experienced recent bottlenecks (Table 1).

Relatedness and inbreeding
Similar relatedness (r) values were obtained from ML-Relate and COANCESTRY. The overall Pearson's correlation coefficient between estimators was 0.958, and it ranged from 0.955 to 0.974 for individual populations (Table S3). Therefore, only ML-Relate results are reported here. Overall relatedness was very low in both mainland populations (Fig. 3A), and they did not differ from each other (Kruskal-Wallis p ¼ 0.230). Relatedness was highest among CUIS-Current bobcats (mean r ¼ 0.106, Fig. 3A), but did not differ from Kiawah (Kruskal-Wallis p ¼ 0.803). However, relatedness in both insular populations was higher than in mainland populations (Kruskal-Wallis p < 0.001). In CUIS-Founders and SC-Mainland, over 86% and 93%, respectively, of pairwise comparisons between individuals produced r-values close to or equal to zero, indicating the individuals were unrelated (Fig. 3B). In contrast, for CUIS-Current, 76.1% of 703 possible pairwise comparisons were unrelated, 12.7% were most likely to be half-siblings, 2.4% were likely full-siblings, and 8.8% were parent-offspring. Ten of these pairs of individuals had r ! 0.60, with three pairs of full-siblings and one parent-offspring pair having high r-values of 0.670e0.803, suggesting they are inbred. On Kiawah, of 630 possible pairwise comparisons, ML-Relate estimated that 77.5% of pairs were unrelated, 12.9% were half-siblings, 3.5% were likely full-siblings, and 6.2% were parent-offspring. Only two pairs of fullsiblings had high r-values of 0.652 and 0.698 (Fig. 3B).
The coalescent-scaled time of splitting was estimated as 0.006 (95% CrI ¼ 0.002e0.032). These coalescent-scaled parameters were converted to N e and t, respectively (Fig. 4a), assuming a mutation rate of 10 À3 mutations per locus per year and a generation time of 4 years, and using the inheritance scalar 2.1797 estimated from the data in IMa2. This gave us effective population sizes of 691 (95% CrI ¼ 181e2921) for CUIS-Founders and 14.13 (95% CrI ¼ 3.51e61.47) for CUIS-Current. The time of splitting, which we interpret as the time of founding for the Cumberland Island population, was estimated at 13.24 years (95% CrI ¼ 3.43e70.45).
Population parameters were also estimated for the naturally-maintained Kiawah system using IMa2 (Fig. 4b).

Abundance and survival estimates
For both male and female analyses, abundance estimates were best modeled (DAIC c ¼ 0) when density differed each year but g0 and s were constrained to be constant across years. At the conclusion of the reintroductions on Cumberland Island, Diefenbach et al. (1994) estimated 29 bobcats on the island but abundance could have ranged from 27 to 40 depending on fates of offspring born in 1989. In 2012 we estimated 9.9 bobcats on the island, but abundance increased over time to 23.9 by 2019 (Table 2). Abundance estimates for 2016e2019 were within 1 standard deviation of the average predicted population size for the second scenario PVA (inbreeding depression), and the abundance estimate for 2012 was 1.3 SD below average.

Population viability analysis model results
By 22 years post-reintroduction (2012) the population averaged about 25 bobcats under the baseline demographic scenario (Fig. 5). With inbreeding added to the demographic scenario, the average abundance steadily declined, and by year 50 there was almost a 20% chance of extinction ( Fig. 5; Fig. S1).
By year 30, average genetic diversity declined from about 0.72 to below 0.6, in which inbreeding depression slightly reduced genetic diversity compared to the demographic scenario (Fig. 5). Adding 1 female bobcat every 4 years stabilized average genetic diversity, regardless of whether an island bobcat was removed (Fig. 5). Adding 2 female bobcats every 2 years and removing an island-born female increased average genetic diversity to about 0.65 but it took 20e30 years. Adding 1 male every 4 years did little to improve population viability and was the least effective translocation strategy for increasing genetic diversity and population persistence (Fig. 5). Scenarios where 1e2 female bobcats were translocated to the island every 2e4 years beginning in year 30 ensured population persistence was >95% to year 100.

Discussion
Although PVAs have been incorporated into the management plans for large and endangered felids such as the Florida panther (van de Kerk et al., 2019), we are not aware of any study of a mid-sized felid where the population was not already endangered and that documented changes in demographics and genetics over time such that it could be used to test PVA predictions. Population estimates 23e30 years after 32 bobcats were reintroduced to Cumberland Island in 1988e1989    (Diefenbach et al., 1993) were consistent with modeled predictions at the time of reintroduction (Diefenbach 1992) and with our PVA. Abundance declined to approximately 10 individuals in 2012 but doubled to 24 bobcats by 2019, and these estimates were within 1 SD of the average of simulated scenarios, except for 2012 (1.3 SD below average). Also, observed heterozygosity after 30 years (0.634, Table 1) was slightly greater than the average predicted heterozygosity but was within one standard deviation (Fig. 5). Parameter values used in the PVA were based on data collected on the bobcat population immediately after reintroduction on the island (Diefenbach 1992;Diefenbach et al., 1993). The original CUIS-Founders from mainland Georgia had high observed heterozygosity (H O ¼ 0.718), and moderately high allelic richness (AR ¼ 5.371), similar to those of present-day mainland South Carolina bobcats (Table 1) and to other continental populations (Anderson et al., 2015;Croteau et al., 2012;Jane cka et al., 2007;Lee et al., 2012;Millions and Swanson 2007;Reding et al., 2012;Smith et al., 2020). Previous studies of bobcats from the southeastern United States, for example, have reported H O ¼ 0.736 and AR ¼ 6.49 (Reding et al., 2012). Palomares et al. (2012) monitored population parameters of the Iberian lynx in Doñana, Spain for 25 years for evidence of inbreeding depression, and found no change in overall mortality rates, but did document reduced litter sizes, increasing female-biased sex ratios, and reduced heterozygosity. On Cumberland Island, survival estimates from our study were consistent with survival estimates reported for nonhunted populations of bobcats (Anderson and Lovallo 2003). The founder bobcats had annual survival rates of 0.93, which was higher than we estimated for 2012e2019, probably because reintroduced bobcats were vaccinated for feline diseases (Diefenbach et al., 1993). The sex ratio of founder bobcats was nearly equal, but during 2018e2019 we found the sex ratio was female-biased and differed from 1.0.
The insular bobcat populations we studied, CUIS-Current and Kiawah, had lower allelic richness and expected heterozygosity than mainland populations, and all measures of genetic diversity were lower in CUIS-Current than Kiawah (Table 1). Furthermore, the observed heterozygosity of CUIS-Current bobcats appears to have declined since we examined this population in 2012 (H O ¼ 0.742 ± 0.074 SE; Diefenbach et al., 2015). In the absence of prompt management intervention, such as additional translocations, this small, isolated population of bobcats on CUIS could be adversely affected by genetic drift, and genetic diversity will continue to decline.
Despite their small population size and isolation from the mainland, however, the overall levels of genetic diversity of the CUIS-Current and Kiawah bobcat populations were higher than has been estimated for populations of the Canada lynx on the islands of Newfoundland (AR ¼ 2.92, H O ¼ 0.40) and Cape Breton (AR ¼ 2.83, H O ¼ 0.41; Prentice et al., 2017), and for isolated populations of other mid-size felids that have experienced significant population bottlenecks, such as the Iberian lynx (Casas-Marce et al., 2013) and the Eurasian lynx (Sindi ci c et al., 2013). Two remnant Iberian lynx populations in Spain have very low allelic richness (2.10 and 3.28), observed heterozygosity (0.31 and 0.46), and effective population size, consistent with demographic contraction and isolation (Casas-Marce et al., 2013). Translocations between these two Iberian lynx populations began in 2007 (Casas- Marce et al., 2013), and second-generation reproduction has been documented (Sim on et al., 2012). However, concerns about genetic erosion in the Iberian lynx remain (Palomares et al., 2012), and long-term conservation efforts may also require supplementation with genetically-healthy captive-bred individuals (Casas-Marce et al., 2013). Similarly, the Eurasian lynx was extirpated in western Europe by the mid-20th century and only four remnant populations remained (Linnell et al., 2009). However, improved ecological conditions, legal protection, and successful reintroduction programs halted the population decline and facilitated recovery (Bagrade et al., 2016;Linnell et al., 2009). Although observed heterozygosity remains low (H O ¼ 0.47) in reintroduced populations of Eurasian lynx in the Dinaric Mountains in Slovenia, Croatia, and Bosnia and Herzegovina (Sindi ci c et al., 2013), moderate levels (H O ¼ 0.52e0.57) have been reported in autochthonous populations in Latvia (Bagrade et al., 2016). It should be noted, however, that these remnant Iberian and Eurasian lynx populations have been purging genetic diversity for significantly longer than the CUIS and Kiawah bobcats, which could explain why their genetic diversity has remained lower than the insular populations we studied.
We believe the PVA we developed could be a useful model for assisting with conservation decisions to maintain the viability of the bobcats on CUIS because the PVA matched both the population dynamics and genetic characteristics of empirical estimates of abundance and genetic diversity. The CUIS population will always be small and isolated and will require continued human intervention to be sustained. However, effective management actions could establish goals and objectives that are likely to result in acceptable population viability. We believe the genetic characteristics of Kiawah bobcats could provide insights for setting goals and objectives for the CUIS bobcat population.
We discovered similar genetic characteristics between the CUIS-Current and Kiawah populations, even though the Kiawah population is not completely isolated from mainland bobcats. The coalescent analysis indicated Kiawah and CUIS populations had similar time-to-founding estimates of~13 years based on CUIS-Current samples and~14 years for the Kiawah population. Although the time-to-founding was underestimated for CUIS by a decade, we assumed errors in the mutation rate and generation time were similar between populations. If these errors affected both data sets similarly, then the true divergence time between the Kiawah and mainland South Carolina populations is likely older than our estimate. Average overall relatedness between CUIS-Current bobcats was moderately low, approaching that associated with first cousins (r ¼ 0.125), had not increased since 2012 (Diefenbach et al., 2015), and was similar to the Kiawah population (Fig. 3). Although the majority (>76%) of individuals on both islands were unrelated (Fig. 3), several individuals had high relatedness (r >> 0.5). Of note, three pairs of full-siblings and one parent-offspring pair on CUIS had r-values of 0.670e0.803, suggesting these individuals were inbred (Fig. 3). On Kiawah, we identified two pairs of full-siblings that had r-values of 0.652 and 0.698 (Fig. 3). The Garza-Williamson modified index was low for both CUIS-Current (M ¼ 0.591) and Kiawah (M ¼ 0.600), in which values of M < 0.68 are suggestive of a population bottleneck (Garza and Williamson 2001). Furthermore, the inbreeding coefficient, F IS , was negative for both CUIS-Current (À0.073) and Kiawah (À0.086), indicating a significant excess of heterozygotes. Although inbred populations are typically expected to have positive values of F IS close to 1, bottlenecked populations may develop a transient excess of heterozygotes at selectively neutral loci following a population bottleneck (Cornuet and Luikart 1996). This leads to negative F IS values that last for a few generations until a new equilibrium is established between the loss of alleles through genetic drift and the generation of new alleles by mutation (Cornuet and Luikart 1996). The F IS for CUIS-Current has increased (i.e., become less negative) since we examined the population in 2012 (F IS ¼ À0.255) and is now similar to that of Kiawah.
We were able to evaluate the natural rate of effective migration in the Kiawah Island system using coalescent-based analyses (IMa2). Using a model with bidirectional migration, we estimated low rates of genetic migration, with immigration onto Kiawah estimated at an average of 0.76 individuals every 4 years, and emigration from the island particularly low (Fig. 4b). Likelihood ratio tests of nested isolation-with-migration models indicated that we cannot reject a model with zero emigration from Kiawah (p ¼ 1), but there was not clear support for a model with no immigration from the mainland into the island (p ¼ 0.092). The persistence of a moderately diverse bobcat population on Kiawah in the absence of managed translocations is supporting evidence of a low but nonzero rate of immigration into the island. The results from Kiawah suggest that bobcats will naturally disperse into accessible coastal islands, but that these islands are likely to act as demographic sinks.
For the CUIS-Current population, estimates of effective population size were generally consistent among methods (11e16.4 individuals) and similar to those reported for highly fragmented populations of the Iberian lynx, which ranged from 8.5 to 13.2 in Doñana, and 15.2e23.1 in Andújar (Casas-Marce et al., 2013). Our temporal estimate of 33 effective breeders on CUIS was somewhat larger than estimates from single-sample methods, but confidence intervals associated with all estimates overlapped in many cases. We acknowledge that the estimate using the heterozygote excess method was not reliable for CUIS-Current because the confidence interval included infinity, but we note that the point estimate was consistent with the other methods. In addition to the different genetic-based estimators being consistent with each other, they also were consistent with our demographic estimates of abundance ( Table 2). The discrepancy between b N e and b N on CUIS was likely due to the effect of both inbreeding and genetic drift acting in a small and isolated population to decrease key measures of diversity, such as allelic richness and heterozygosity, and thereby decrease b N e .
Genetic estimates of b N e for Kiawah were similar among estimators, 5.6e9.6 individuals, or had large confidence intervals (linkage disequilibrium estimator, b N e ¼ 21.7). Our estimates of b N e for the self-introduced population on Kiawah were similar to those for Canada lynx on much larger islands (Newfoundland N e ¼ 6.5, Cape Breton Island N e ¼ 7.6; Prentice et al., 2017).

Point estimates of b
N e were generally lower for the Kiawah population compared to the CUIS-Current, but not statistically different. Nevertheless, Kiawah (3200 ha) is less than half the size of CUIS (6935 ha) and overlap of female home ranges, which may limit successful reproduction (Diefenbach et al., 2006), would occur at lower abundance and result in a smaller effective population size.
Without any intervention, the PVA predicted further decline in abundance and genetic diversity (Fig. 5) and increased likelihood of extinction ( Supplementary Fig. S1). On Kiawah, we estimated one bobcat immigrated to the island every 5.3 years, which would be similar to our scenario of translocating one bobcat to CUIS every 4 years. Under this model, we predicted that genetic diversity on CUIS would be maintained at a level similar to Kiawah and that the probability of extinction would be near zero ( Supplementary Fig. S1, Fig. 5).
Higher rates of translocations (e.g., every 2 years) resulted in greater levels of heterozygosity ( Fig. 5) but still lower than that of the founding population. Furthermore, translocations of males did increase genetic diversity but were less effective at reducing the risk of population extinction than translocating females ( Supplementary Fig. S1). Given that the current CUIS population has limited evidence of inbreeding, it is possible that translocations of one female bobcat every 4e6 years would be sufficient to sustain the viability of the population.
If translocations of bobcats were to occur on CUIS, as part of a properly designed research and monitoring program, we believe there are several areas of uncertainty that could be reduced regarding the effectiveness of translocations, which could also be applied to endangered felids. First, GPS technology is available to monitor the movements of bobcat-sized animals. A greater understanding of space use could be obtained by monitoring movements of CUIS resident bobcats and translocated individuals, especially immediately after translocations occurred. Diefenbach et al. (2006) found significant overlap in space use when there were >30 bobcats on the island, which may affect reproduction if females cannot find areas of exclusive use. Understanding where translocated individuals establish home ranges and how resident bobcats respond to new individuals could help identify whether simply adding individuals when populations are low is more effective than periodic translocations with replacement (Fig. 5). Second, monitoring the survival, spatial ecology, and breeding success of introduced and resident females may identify mechanisms by which different translocation strategies are successful. Third, continued monitoring of genetic structure of the population could document the effectiveness of translocation strategies. In particular, do translocated females breed and contribute their genetics to the population, and how quickly does that occur?
The PVA indicated that, regardless of translocation strategy, genetic diversity is unlikely to be restored to levels of the founding population on CUIS. Therefore, it seems likely that some action will be required to ensure the persistence of the population. If no action is taken, the PVA predicted that the probability of extinction increases to 1 in 5 within the next 20 years. Bobcats have secure populations throughout their entire range and so research on bobcats on CUIS does not carry the risks of conservation efforts on globally endangered mid-sized felids, such as the Iberian lynx and Eurasian lynx. Conservation actions taken to sustain the CUIS bobcat population provide an opportunity to learn, and could provide insights into population improvement for endangered mid-sized felids in the future.
Funding Support for this project was provided by The National Park Service, The Pennsylvania Game Commission (Cooperative Agreement # 1434-03HQRU1548), and by Penn State Beaver.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.