Using recent genetic history to inform conservation options of two Lesser Caymans iguana ( Cyclura nubila caymanensis ) populations

The Sister Islands rock iguana ( Cyclura nubila caymanensis ) is critically endangered and endemic to the Caribbean islands Little Cayman and Cayman Brac. The Cayman Brac population and indeed the entire species is under threat from habitat destruction, invasive species, and anthropogenic impacts on the island. We assessed the genetic diversity, estimated effective population sizes, and tested for differentiation of populations between these two islands to inform potential future translocation should this be needed for the Cayman Brac population. Two mitochondrial DNA markers (cytochrome b and NADH subunit 4) and seven DNA microsatellite markers were used to assess the genetic diversity, genetic structure, demographic history, and effective population size of the two iguana populations. Mitochondrial DNA showed no genetic differentiation between populations; however, we found little to moderate divergence with microsatellites. We compared multiple demographic scenarios and revealed that ongoing gene flow is likely. The demographic history implied a significant genetic bottleneck around 10,000 years ago, coinciding with the sea level rise at the close of the last glacial period, and the start of the Holocene. Estimates of current effective population sizes indicate a small-scale number of breeders on each island of similar magnitude to the census mature population size (between 100 and 800 individuals). The relatively low differentiation between populations supports the possible development of active genetic management plans to manage the declining populations of the Sister Islands rock iguana.


Introduction
Islands are centres of interest in the global biodiversity crisis, with disproportionately higher rates of endemism and extinction when compared with continents, making them key targets for conserving biodiversity (Kier et al. 2009;Spatz et al. 2022).The Caribbean basin is one of 25 worldwide biodiversity hotspots, and despite covering only 0.15% of the Earth's total land area, supports 6.3% (520 species) of the world's known reptiles (Hedges 2006;Habel et al. 2019;Spatz et al. 2022).Within the Class Reptilia, Iguanas are the most endangered, with 77% recorded as near threatened, or threatened (Moss et al. 2018;Perry et al. 2022).West Indian rock iguanas (genus Cyclura), comprising nine extant species and six subspecies, are the most endangered group of iguanas worldwide (Moss et al. 2020;Buckley et al. 2022;Reynolds et al. 2022).Within this genus, Cyclura nubila, comprises two subspecies, the Cuban rock iguana C. nubila nubila, and the Sister Islands rock iguana, C. nubila caymanensis (Buckley et al. 2022).Cyclura nubila caymanensis is listed as Critically Endangered by the International Union for Conservation of Nature (IUCN), with a projected decline of 95-98% over 34 years if conservation efforts are not implemented immediately (Goetz and Burton 2012).Cyclura nubile caymanensis is endemic to Little Cayman and Cayman Brac, two islands located approximately 100 km from Grand Cayman, Cuba and Jamaica (Moss et al.Thea F. Rogers and Ewan H. Stenhouse have contributed equally to this work.2019).Little Cayman and Cayman Brac are approximately similar in size and shape (~ 20 km long and ~ 2 km wide) and separated by 7.5 km.
Cyclura nubila caymanensis populations have been declining primarily as a result of habitat degradation, introduction of invasive species, and direct anthropogenic impacts such as hunting (Moss et al. 2018;Shaney et al. 2020;Perry et al. 2021).Urbanization on Cayman Brac has caused habitat destruction and disruption of sandy beaches used as nesting sites by iguanas (Goetz and Burton 2012).The situation on Little Cayman is less severe, supporting between 1200 and 1500, or 2000 and 4100 iguanas by minimum and maximum estimates, respectively (Goetz and Burton 2012;Rivera-Milán et al. 2015).However, an estimated 6-10% of the C. nubila caymanensis subpopulation are killed by road traffic annually, particularly during the breeding season when gravid females migrate approximately 2 km to access coastal nesting sites (Goetz and Burton 2012;Moss et al. 2020).The invasive green iguana (Iguana iguana), widely thought of as a pest due to the damaging impacts of high population densities on local biodiversity, have been present on Little Cayman and Cayman Brac since 2007 (Moss et al. 2018).In August 2016, three hybrid hatchlings were discovered on Little Cayman Island (Moss et al. 2018;Knapp et al. 2021;Perry et al. 2021).When hybridisation occurs between rare taxa and a widespread sister species, native taxa are at an increased risk of replacement due to competition, but also genetic swamping due to introgression (Moss et al. 2018).Therefore, this discovery may have far reaching implications for biosecurity within the Caribbean (Moss et al. 2018).As with other Cyclura species, C. nubila caymanensis are herbivorous, with repetitive grazing on vegetation by iguanas encouraging additional shoot and foliage growth, promoting plant community diversity (Pasachnik and Martín Vélez 2017).Thus, C. nubila caymanensis may play a key role within Cayman Sister Islands, highlighting the need for the conservation of this species.For small, isolated populations to persist, best practice is to incorporate active genetic management into conservation, twinned with actions that address other threats, including habitat loss, and improved demographic viability (Ralls et al. 2017).
The small population size of the Sister Isles rock iguana makes it especially susceptible to stochastic loss of genetic diversity, as well as inbreeding, due to limited selection and restrictive mating (Moss et al. 2019).This can result in inbreeding depression, reducing fitness of inbred individuals, and can reduce adaptive potential of individuals over time (Moss et al. 2019).To potentially limit the loss of genetic diversity in small populations, active genetic management, including the transfer of individuals from a source population to a translocated population can be implemented, as undertaken with a translocation of Bahamian iguanas (Cyclura cychlura inornana) (Knapp and Malone 2003).
Current C. nubila caymanensis conservation strategies have, until now, primarily involved protecting and restoring habitat (Binns et al. 2011).Given evidence that translocations of Bahamian iguana populations resulted in individuals successfully breeding (Knapp and Malone 2003), it is appropriate to consider active genetic management between Little Cayman and Cayman Brac as a potential conservation strategy.However, there should be caution when exploring active genetic management.Concern still remains regarding outbreeding depression and genetic homogenisation (Edmands 2007;Bell et al. 2019), and migrants can introduce beneficial and deleterious genetic variation (Bell et al. 2019).Beneficial impacts of management can mask deleterious, recessive alleles, and increase genetic variation (Whiteley et al. 2015).However, deleterious impacts of gene flow can be caused by a reduction in local adaptation or genetic incompatibilities between populations (Bell et al. 2019).For these reasons, an assessment of the genetic affinities of the populations needs to be undertaken to ascertain whether active genetic management would be beneficial or detrimental.To implement conservation strategies, it is vital to estimate the effective population size, as species with low effective population sizes are susceptible to losing genetic diversity quickly following genetic drift and inbreeding depression (Weeks et al. 2016).
In this study, we quantify the contemporary genetic diversity of two populations of C. nubila caymanensis and test whether they are differentiated, explore their demographic history, including estimating effective population sizes around the year 2010, and estimate the time of divergence between the two islands.We used explicit demographic simulations under the coalescent framework to evaluate whether the two island populations evolved isolated from each other or are connected by ongoing gene flow, and whether the present-day genetic variation is the outcome of a stable or dynamic demographic history.The results from these analyses reflect a dynamic demographic history in the Lesser Caymans iguanas and are relevant for future conservation, including potentially active genetic management.

Samples
Totals of 21 blood samples from Cayman Brac and 41 blood samples from Little Cayman were collected by the Durrell Wildlife Conservation Trust and the Cayman Islands Department of Environment between October 2007 and September 2011 (Supplementary Material -Table S1).Iguanas were caught by noose or baited live animal traps (Hava-hart®, Woodstream Corp., Lititz, PA 17543) and individually marked with bead tags through the dorsal crest (Rodda et al. 1988) and passive integrated transponders (PIT-tag; ISO FXD; Trovan, Electronic Identification Devices Ltd.) subcutaneously at the dorsal tail base.A quantity of blood was drawn from the tail vein and 0.3 mL stored in 0.8 mL Lysis buffer (800 mM guanidine hydrochloride; 30 mM Tris Cl, pH 8.0; 30 mM EDTA, pH 8.0; 5% Tween 20; 0.5% Triton X-100), and 0.3 ml in 0.8 mL 99% ethanol.
We extracted DNA using a modified salting-out DNA extraction protocol (Supplementary Material, Methods), following (Miller et al. 1988).We vortexed samples for 1 min.Following this, we removed 100 μL of the supernatant, which was added to a 1.5 mL microcentrifuge tube along with 10 μL of proteinase K (20 mg/mL) and 580 μL of TNES buffer (0.273 g Tris, 0.316 g NaCl, 0.263 g EDTA, 0.45 g SDS diluted in 45 mL of polished water).We vortexed the supernatant again for a minimum of 1 min per tube, or until homogenised.We then incubated the samples in an incubating microplate shaker at 55 °C and 250 rpm for ≥ 12 h.Following incubation, we vortexed the samples for 30 s.We added 170 μL of 5 M NaCl (45 ml: 13.15 g NaCl, 45 ml polished water) to each sample which were then vortexed.We centrifuged samples at 13,300 rpm for 6 min, after which we removed the supernatant and split equally between two new microcentrifuge tubes.We added equal volumes of cold 100% ethanol to each new tube, and incubated for ≥ 1 h at −20 °C.Following incubation, we centrifuged the samples at 13,300 rpm for 17 min, after which we discarded the ethanol and washed the DNA pellet with 200 μl 100 mM ethanol sodium acetate (0.347 g NH4 + dissolved in 31.5 ml of 100% ethanol and 13.5 ml of polished water).We centrifuged the samples at 13,300 rpm for 5 min, discarded any remaining ethanol, and air dried the DNA pellets.When completely dry, we added 100 μL of sterile water to each tube to resuspend the DNA and incubated overnight at 4 °C.

Mitochondrial sequencing
We amplified and sequenced two markers: cytochrome b (cyt b) and NADH subunit 4 (ND 4), as these markers are considered highly suitable for assessing mtDNA variation in iguanids (Malone et al. 2000;Welch et al. 2017).Primers were designed on Primer3 v. 4.0 (Koressaar and Remm 2007) based on Cyclura mitochondrial DNA (mtDNA) sequences available in GenBank (Supplementary Material, Table S2).Polymerase Chain Reaction (PCR) was undertaken separately for each fragment in 25 μL volumes consisting of 4 μL 1X Green Flexi Buffer (Promega), 3 μL 25 mM MgCl 2 , 0.4 μL 10 mM dNTPs, 0.3 μL of each of the 20 μM primers, 0.12 μL 5U/μL Taq DNA polymerase (Promega), 1 uL DNA, and 15.18 μL nuclease free water.PCR cycle conditions for both markers were: initial denaturing stage of 95 °C for 15 min, 32 cycles of 94 °C for 1 min, annealing temperature 58 °C for 1 min, extension of 72 °C for 1 min, and a final extension of 72 °C for 10 min.Mitochondrial PCR products were sent for Sanger sequencing to Eurofins Genomics Services, Germany.
We analysed mtDNA sequences using Sequencher v. 5.1 (as used in Steinfartz et al. (2009)) to generate consensus sequences from the forward and reverse sequences for each individual and gene.Cyt b and ND4 sequences were aligned in Mega v. 7.0 and ClustalW multiple sequence algorithm, and further adjusted by eye (Larkin et al. 2007;Tamura et al. 2013).Three separate datasets were generated for downstream analyses: i) cyt b sequences only, ii) ND4 sequences only, and iii) the concatenated dataset of both fragments.Final datasets only included sequences that had a Phred sequence quality of over > 20 and were over 600 base pairs long.
Mitochondrial DNA: phylogenetic analysis and population structure DnaSP v. 5.0 (Librado and Rozas 2009) was used to produce summary statistics for both cyt b and ND4 sequences.We used MrModeltest v. 2.3 (Steinfartz et al. 2009) to determine the most suitable model of sequence evolution for the datasets.A maximum likelihood (ML) tree was generated in Mega for the concatenated dataset with 100 bootstrap replicates to assess nodal support.The outgroup Iguana iguana was used to root the tree, as used in other studies of Cyclura sp.(Shaney et al. 2020).We also generated median joining networks of haplotypes from each dataset using Popart v. 1.7 (Leigh and Bryant 2015).

Mitochondrial DNA demographic history
To explore the demographic history of the two island iguana populations, pairwise mismatch distributions (MD) were calculated with Arlequin v. 3.5.2.2 (Excoffier and Lischer 2010).Mismatch distributions are the number of pairwise differences among all DNA samples within a sample and can be explored to better understand population growth.When plotted on a histogram, in a stationary growth population, distributions from nonrecombinant DNA sequences become irregular, whereas a growing population generates mismatch distributions that are smooth and have a peak (Harpending 1994).The position of the peak reflects the time of the population growth.Bayesian skyline plots (BSPs) were estimated using BEAST v. 1.8 (Drummond and Rambaut 2007).BSPs were retrieved from the BEAST analysis using Tracer v. 1.6 (Rambaut et al. 2018).The BEAST runs had the following parameters for the Markov Chain Monte Carlo (MCMC) algorithm: length = 50,000,000 steps, with the first 10% discarded as burn-in and nucleotide substitution rate 1.13 × 10 -8 per site per year.The MCMC analysis were judged to have achieved when Effective Sample Size (ESS) values exceeded 200 and via examination of the likelihood values in Tracer.This substitution rate was estimated by calculating the average number of pairwise differences between C. nubila caymanensis and the outgroup Iguana iguana using DnaSP.This value was divided by 15 MY, the known divergence time between the two species (MacLeod et al. 2015;Reynolds et al. 2022).This value was then divided by the sequence length to give the substitution rate per site.Both the MDs and the BSPs were made using the concatenated datasets.

Microsatellite genotyping
The ten most polymorphic DNA microsatellite markers in An et al. (2004) were chosen for PCR amplification: 60HDZ13, 60HDZ50, 60HDZ106, 60HDZ148, 60HDZ151, 60HDZ152, 60HDZ181, 60HDZ419, 60HDZ590, and 60HDZ780.PCR products were fluorescently labelled with 6-carboxy-fluorescine (FAM), or hexachloro-6-carboxy-fluorescine (HEX)) using the M13 protocol in Schuelke (2000), enabling products to be genotypically analysed with a capillary sequencing service using the tail sequence AGG GTT TTC CCA GTC ACG ACG ATT on the forward primers.The microsatellite PCRs had final volumes of 10 μL with 5 μL QIAGEN Multiplex PCR Master Mix, 0.2 μL of each 20 μM primer, 0.2 μL 10 mM fluorescent dye (either HEX, FAM, or 6-FAM), 1 μL DNA and 3.4 μL nuclease free water.For the amplification of 60HDZ181 1 μL of QIAGEN Q solution was added to the PCR reaction.The PCR conditions for the microsatellite markers used were as follows: initial denaturing stage of 95 °C for 15 min, 35 cycles of 94 °C for 30 s, annealing temperature specific to the primer for 90 s, an extension of 72 °C for 90 s, and a final extension of 72 °C for 10 min.Microsatellite PCR products labelled with different fluorescent dyes for the same individuals were pooled together and sent for fragment analysis to DNA Sequencing & Services Dundee, Dundee, UK.GeneMapper v. 4.0 (Chatterji and Pachter 2006) was used to determine microsatellite allele size in base pairs.Individuals with > 20% missing data were removed from the final dataset.

Microsatellite summary statistics
For each population, Micro-Checker v. 2.3.3 (Van Oosterhout et al. 2004) was used to test for null alleles, and Genepop v. 4.5.1 (Rousset 2008) was used to check for linkage disequilibrium (LD) between pairs of markers and for deviations from Hardy-Weinberg equilibrium (HWE).A strict Bonferroni correction was applied to these results to reduce type 1 error due to multiple comparisons (Armstrong 2014).Genetic diversity for each island population was assessed using inbreeding coefficient (F is ) values, observed (H o ) and expected (H e ) heterozygosity, number of alleles per locus and allelic richness all calculated with Microsatellite Analyser (MSA) v. 4.05 (Dieringer and Schlötterer 2003).Two-tailed heteroscedasticity t-tests were undertaken in R (R Core Team 2020) to check for significant differences in genetic diversity between islands.FSTAT v. 2.9.3.2 (Goudet 1995) was used to test if Fis values were significantly different from zero.

Microsatellite population divergence
MSA v. 4.05 (Dieringer and Schlötterer 2003) was used to estimate the pairwise fixation index (F st ) between the two islands and a factorial component analysis (FCA) was created to explore the proportion of variance between populations using the program GENETIX v. 4.05 (Belkhir et al. 1996).We analysed population structure using Bayesian structure analysis in STRU CTU RE v. 2.3.4 (Hubisz et al. 2009), run with the following parameters: MCMC burn-in length = 20,000, number of repeats after burn-in = 50,000.Values for the number of clusters (K) tested ranged from K = 1 to 10, with each analysis repeated 10 times, as used in previous studies of Iguanids (van den Burg et al. 2018;Zarza et al. 2019).We determined the most likely value of K in Structure Harvester (Earl and Von Holdt 2012) using the Evanno method (Evanno et al. 2005).All runs for most likely K were combined using Clumpp (Jakobsson and Rosenberg 2007).2Mod v. 0.2 (Ciofi and Bruford 1999) was used to determine whether the two island populations evolved under a model of genetic drift (isolation) or a model of ongoing gene flow (migration).2Mod was run three times for 500,000 MCMC iterations, with the first 20% discarded as burn-in.A score of 1 or 0 was assigned to each iteration, with 1 representing migration between populations and 0 representing isolation.The posterior probability of each evolutionary model was estimated by calculating the percentages of 1 and 0 in the output of each run and averaged across runs.The average number of migrants per generation (N m ) was estimated using the private allele method (Slatkin 1985) in Genepop and also estimated from the F st calculated with MSA (Dieringer and Schlötterer 2003).

Microsatellite demographic history
To test whether a population expansion had occurred since the last glacial maximum, macro K and G tests (Reich et al. 1999) were performed in Excel.These tests were complimented by an analysis of the population bottlenecks using the software Bottleneck v. 1.2.2 (Cornuet and Luikart 1996).Simulations were carried out under the two-phase model of microsatellite evolution (TPM) (Di Rienzo et al. 1994) with default parameters, i.e. 70% of the mutations were singlestep and with multistep mutations of up to seven tandem repeats.MSVar v. 1.3 (Beaumont 1999) was used to characterise the effective population size of each population and the timing at which the potential demographic changes may have occurred.MSVar simulations were conditioned with four different demographic scenarios following Gonzalez et al. (2014) to assess lack of dependency of the posterior estimates on the model parameters and starting points.We tested for signatures of a bottleneck, expansion, and stable scenario for each population.Each MSVar run consisted of 200 × 106 simulations of the MCMC algorithm with the first 20% discarded as burn-in.Convergence of the chains was assessed with Gelman and Rubin's diagnostic (GR) (Brooks and Gelman 1998) calculated on the three runs performed for each population with different prior scenarios.GR was conducted using the CODA library (Plummer et al. 2006) in R (Venables et al. 2009).

Effective population size (N e )
To examine population size, we estimated the historical N e based on the expected heterozygosity and calculated theta (Ohta and Kimura 1973).For these analyses three different mutation rates were assumed, i.e., 1 × 10 -3 , 1 × 10 -4 (microsatellite mutation rates used for genetic analysis of Komodo dragons, Varanus komodoensis in Ciofi and Bruford (1999) and 5.37 × 10 -4 (calculated for the microsatellite dataset using MSVar)).

DIYABC analysis
In addition to the above analysis, we undertook approximate Bayesian computation (ABC) analysis.This allows for historical inferences of complex population histories, including divergence events and changes in past population sizes.DIYABC allowed estimated time of divergence to be calculated between Little Cayman and Cayman Brac.The analysis was run with 1,000,000 simulations and a mutation rate of 1 × 10 -4 (Komodo dragon microsatellite mutation rate) with other parameters set as default and adapting the ranges in prior distributions of the effective population size and the divergence time so they would include the complete posterior distribution (Ciofi and Bruford 1999).

Summary statistics of genetic variation
Fragments of 754 and 705 base pairs in length were successfully amplified for the cyt b (47 samples) and ND4 (56 samples) markers, respectively (Table 1).The concatenated dataset contained 46 sequences, 20 from Cayman Brac and 26 from Little Cayman (1459 bp in length).A total of 11 haplotypes were identified in the cyt b dataset and 10 in the ND4.For the concatenated dataset 10 haplotypes were found.Pearson's chi-squared test showed no significant differences in frequency of haplotypes between markers or populations (p = 0.83 and p = 0.28 respectively), nor in the number of polymorphic sites between markers or populations (both p = 1).Haplotype diversity and average number of pairwise nucleotide differences were higher for cyt b than for ND4 marker, with the Little Cayman population showing an apparent higher diversity than the Cayman Brac population, although these differences were not statistically significant.Theta (W) was higher in cyt b than ND4 in each population.Theta values for each marker were similar in Cayman Brac and Little Cayman (Table 1).
Out of ten microsatellites tested, seven generated amplicons that could be genotyped.We thus excluded loci 60HDZ106, 60HDZ590 and 60HDZ780 from all downstream analysis.Micro-checker software suggested a potential null allele in marker 60HDZ50 (Oosterhout statistic = 0.213), and Genepop indicated significant deviations from HWE due to heterozygote deficiency for this marker in both populations (Cayman Brac p = 0.041, Little Cayman p = < 0.001).Deviations from HWE for this locus were therefore most likely due to a null allele rather than the effects of inbreeding, which would cause deviations from Hardy-Weinberg across multiple loci.Marker 60HDZ50 was still used in the analyses of the genotyping results as t-tests indicated that summary statistics (e.g., observed heterozygosity) were not significantly different when this microsatellite was removed from the dataset (Supplementary Material, Table S3).No pairs of microsatellites were found to be in LD.
Values for H e , the average number of alleles per locus, and allelic richness were higher in Little Cayman than Cayman Brac, and H o was higher in Cayman Brac than Little Cayman (Table 2).However, t-tests revealed no significant differences in these values between the two populations.The population in Cayman Brac presented an inbreeding coefficient (F is ) close to zero (−0.0029), while Little Cayman presented a higher inbreeding coefficient (0.0703) indicating inbreeding.Nevertheless, a statistical test with the software FSTAT indicated that these F is values were not statistically different form zero (p > 0.05) (Fig. 1) (Goudet 1995).

Population structure
The pairwise F st estimated with the microsatellite markers was approximately 8% (Table 2) indicating little to moderate population differentiation between the two islands.This pattern of differentiation was reflected in the FCA (Supplementary Material, Fig. S4), where the two populations appeared to form separate groups according to the two main factors explaining the variation in the data (16.97%).The STRU CTU RE analysis supported five clusters as the best clustering solution of the data using Evanno's method (Supplementary Material, Fig. S5).However, a visual examination of the bar plots produced by STRU CTU RE for K = 5 showed no biological meaningful results, as most individuals shared memberships to each cluster (Fig. 2A).When  conditioning the data to K = 2, the Cayman Brac individuals were assigned to cluster one with most of their genotype belonging to that cluster (red in Fig. 2B).In contrast, most of the Little Cayman individuals were assigned to cluster 2, with much of their genotyping belonging to that cluster (green in Fig. 2B).Nevertheless, individuals 39, 42, 45, 67, 68, 70, 79, 82 and 83 from Little Cayman were assigned to cluster 1 with the Cayman Brac samples, as their genotypes more closely resembled those in Little Cayman than the Cayman Brac (Fig. S6 in SM shows the individual labels).
Overall, this analysis indicates minimal difference between the genetic makeup of the iguana populations.
An explicit coalescent simulation of a model of evolution under genetic drift versus a model including gene flow between the two islands using the software 2Mod.This showed that the gene flow model is two times more likely to explain the differences in allelic frequencies between populations (i.e., 64.5% vs 35.5%, respectively).Consistent with these results, estimates of the number of migrants between the two islands (N m ) based on the pairwise F st and on the mean frequency of private alleles in each population (0.076) results in one migrant per generation between the two populations.

Demographic history and effective population size
Mitochondrial DNA MDs (Supplementary Material, Fig. S3) showed no population expansion on either island, as both distributions were multimodal and not unimodal as expected under a demographic expansion.Despite the observed multimodality in both islands, only the Little Cayman population showed a significant deviation from the expected model when tested with the sum of square deviations statistics (p = 0.04).The remaining measures of fit to an expansion model were not significant (p > 0.05).In the Cayman Brac population, the main peak in the observed distribution was slightly steeper than the expected line, suggesting a bottleneck (Steinfartz et al. 2009).The BSPs are similar for both populations, showing a mild decline in effective population size in both islands between 500,000 and 100,000 years ago, and a stark decline starting approximately 100,000 years ago resulting in an effective population size decreased by more than one order of magnitude (Fig. 3A).
The microsatellite demographic analysis using the K and G tests detected no expansion in either population (Cayman Brac K test p value = 0.75, G test = 0.46 and for Little Cayman K test p value = 0.75, G test = 0.64).In contrast, we detected evidence of a genetic bottleneck in the Little Cayman population using the heterozygosity excess method (Wilcoxon test p = 0.0008).The Cayman Brac population also showed an excess of heterozygosity but was borderline non-significant (Wilcoxon test p = 0.06).
Consistent with trends implied from results of the Bottleneck analysis, MSVar detected a reduction in N e that took place approximately 10,000 years ago (Fig. 3B), and which reduced N e by 2 and 3 orders of magnitude reaching as few  S4).The ratio between the estimated N t (population size at time t) and N 0 (current population size) was consistently larger than 1 reflecting a larger ancestral population size compared to the current population (Supplementary Material, Fig. S7).The MSVar runs using different prior distributions for the parameters measured were convergent, as indicated by the GR statistics, which were lower than 1.2 in each population (1.19 for Cayman Brac and 1.11 for Little Cayman).
The historical N e was estimated from the expected heterozygosity using the mutation rate estimated by MSVar (5.37 × 10 -4 ) as well as two additional mutation rates from the literature (10 -3 and 10 -4 ).These represent rates that were half an order of magnitude faster or slower (Ciofi and Bruford 1999).Based on the population specific inferred rate by MSVar, the N e values estimated were 1552 for Cayman Brac and 2102 for Little Cayman.Using the higher mutation rate results in N e values as low as 833 and 1,129 for Cayman Brac and Little Cayman, respectively, while the corresponding estimates with the lower mutation rate were 8332 and 11,286 (Supplementary Material, Table S5).Estimating N e from the neutral evolution rate parameter theta generated similar N e values to estimates based on expected heterozygosity (Supplementary Material, Fig. S5), consistent with MSVar's estimates of large Ne values prior to the bottleneck.
The results were consistent with Bottleneck and MSVar, detecting a reduction in N e approximately 13,000 years ago.This reduction in population size also implies a Holocene bottleneck.The distribution of effective population sizes (N e ) shows that the ratio of ancestral to modern population sizes supports a bottleneck with overlapping ranges of the estimated occurrence of the bottleneck.Estimates of expected heterozygosity are similar across all methods -both MSVar and DIYABC analysis results overlap (Supplementary Material, Table S6).

Phylogenetic reconstruction
MrModelTest identified the General Time Reversible model (GTR) as the most suitable model of sequence evolution for the concatenated mtDNA dataset.Using this substitution model, the ML tree inferred for this dataset resulted in haplotypes shared by individuals of both islands with no island specific clade (Supplementary Material, Fig. S1).The tree topology was dominated by small branch lengths reflecting the accumulation of few nucleotide substitutions between sequences.This shared distribution of genetic variation between islands was also observed with the haplotype network (Fig. 1).The same result was obtained when each mitochondrial DNA dataset was analysed separately (Fig. S2 SM), preventing identifying individuals' geographic provenance based on these molecular markers.

Discussion
Here, for the first time, we explore the genetics of the only two populations of Cyclura nubila caymanensis using mitochondrial and nuclear DNA data.Overall, we found that the two populations of Cayman iguanas present relatively high levels of genetic variation and few differences between them with both the mtDNA and nuclear markers, although the latter seemed to find more differences between the two islands than the materially inherited markers.
Island populations tend to have higher extinction rates than mainland ones, partially due to limited population sizes (Perry et al. 2021).Small census population sizes are often associated with small effective population sizes which exacerbates the effects of genetic drift, but can lead to increased inbreeding (Bell et al. 2019).Additionally, small populations can be impacted by stochastic environmental factors that dramatically modify the entire species range and lead to an "extinction vortex" (Blomqvist et al. 2010).
Cyclones occurring in the Caribbean region can have a radius extending well over 50 km with varying wind speeds reaching over 300 km/h, and as a result can be extremely destructive (Mumby et al. 2011).The two populations studied here present relatively high levels of genetic variation (observed heterozygosity was higher than 60% in both populations) in comparison to other iguanids (Lanterbecq et al. 2010;MacLeod et al. 2015;Moss et al. 2019), and do not show significant inbreeding (F is ).Overall, Little Cayman shows higher levels of genetic variation, although this may simply reflect larger sample size (41) in comparison to Cayman Brac (21).Haplotypes of mtDNA were shared between islands (Fig. 1 and Supplementary Material, Fig. S1).This observation implies either incomplete lineage sorting or gene flow between the two islands.
The microsatellite markers showed low differentiation between the two islands as measured with the pairwise F st (~ 8%; Table 2) and visualised in the FCA analyses.The Evanno analysis of the STRU CTU RE results indicated that K = 5 was the best supported STRU CTU RE result; however, this clustering solution showed that most individuals shared genetic variation across the inferred clusters, thus suggesting that over-splitting the data into clusters has no biological meaningful value (Fig. 2A).The STRU CTU RE result for K = 2 showed a separation between the two islands (Fig. 2B) although individuals in both islands seemed admixed.As microsatellites have a higher mutation rate and represent (putatively) independent loci when compared to the mtDNA, they generate a higher statistical power to uncover genetic differences between individuals, for example microsatellites can resolve genetic differences between siblings, which mtDNA cannot do as it is uniparentally inherited (Lanterbecq et al. 2010).Further, the shared mtDNA haplotypes may also reflect sex-biased gene flow with higher female dispersal than by males (Moss et al. 2019).
To obtain a better understanding of the causes of contemporary patterns of genetic variation in these two populations, we analysed a scenario of evolution under genetic drift against another that included gene flow between the two islands (Ciofi and Bruford 1999).The outcome supported a model where the two islands have experienced historical and, possibly intermittent, gene flow, a result supported by the estimate of one migrant per generation between islands based on the number of private alleles and the infinite island model (Slatkin 1985).Overall, these results (limited nuclear DNA divergence and lack of mtDNA divergence) imply a recent connection between the two islands, which may be explained by the relatively shallow depth of the sea basin between them.Currently Little Cayman and Cayman Brac are separated by ~ 7.5 km of sea with depths varying between a few meters to over 200 m.However, much of the water body between these islands is shallow with a section of approximately 1-2 km of length deeper than 200 m (Kobara and Heyman 2008).Therefore, it is likely that during the last glacial period, while sea levels were on average 121 (± 5) meters lower than today, the inhabitable area on each island was larger and the width of sea separating them was much narrower (Fairbanks 1989).Although the current space between these two islands could be overcome via rafting, closer proximity of the islands during the period of lower sea level is likely to have facilitated this contact, thereby limiting the divergence between the two populations (van den Burg et al. 2023).Additionally, anthropogenic facilitated migration of individual iguanas may have occurred (F.J. Burton, pers.comm).

Demographic history
Prior to the sea level rise towards the end of the last glacial period, Little Cayman and Cayman Brac presented a larger area above sea level, approximately 30% to 50% larger than today (Fairbanks 1989).Circa 11.5 KYA, North American and Eurasian glaciers began to retreat, and sea levels rose until reaching present levels.The reduction of inhabitable island surface concomitantly likely resulted in a decline of the iguana population size, as evident for now extinct vertebrates in the Bahamas (Perry et al. 2021).While we have no certainty that this change in sea level caused the bottleneck evident in the genetic data of the Sister Islands rock iguanas, it is likely that it was the driver as it has been shown to have had equivalent effects on other taxa including iguanas of the genus Cyclura (Aplasca et al. 2016;Reynolds et al. 2022).This bottleneck dramatically reduced the iguana populations, which in turn resulted in a stronger effect of genetic drift and the appearance of the small levels of divergence (F st ) between the Little Cayman and Cayman Brac populations.
A further factor contributing to this bottleneck could be human arrival on Caribbean islands, dated to around 8000 to 5000 years ago, based on evidence from sites in Cuba, Hispañola and Puerto Rico (Moreno-Estrada et al. 2013).Although the mode of the time of the bottleneck estimated from MSVar was ~ 10,000 KYA, the 95% HPD of the time estimate ranges between 3667 and 73,930 years ago in Cayman Brac and 2175 and 61,292 years ago in Little Cayman, resulting in an overlap between the potential times of the bottleneck and the human colonisation of these islands, indicating that more recent bottlenecks could be due to anthropogenic influence.Human colonisation resulting in habitat alteration and the introduction of invasive species has been known to cause dramatic declines in other island populations of reptiles (Moss et al. 2018;Perry et al. 2021).Gonzalez et al. (2014) found a significant genetic bottleneck in populations of the La Gomera giant lizard (Gallotia bravoana) in the Canary Islands that occurred approximately 2500 years ago, coinciding with first human arrival, their modification of the habitat, and the introduction of nonnative species.

Effective population size (N e )
Across analyses we observed that the N e values for Cayman Brac were smaller than those for Little Cayman.This is in line with the higher genetic diversity in Little Cayman (e.g., expected heterozygosity, average number of alleles per locus, and allelic richness).The MSVar estimates of current N e (between 100 and 500 animals; Fig. 3D) are consistent with the estimated current numbers of sexually mature iguanas on these islands, i.e., approximately 100 individuals for Cayman Brac and 800 individuals for Little Cayman (Goetz and Burton 2012).The combination of small N e , an N e /N c possibly close to one, and the small islands these reptiles inhabit puts them in a precarious situation where genetic drift may quickly remove genetic variation (Aplasca et al. 2016), and catastrophic or demographically stochastic changes may result in population collapses (Hudson et al. 2016;Valbuena-Ureña et al. 2017).Our values of N e differed slightly from the published estimate of 71.8 published in (Moss et al. 2019).As our N e estimates were much higher, this implies that populations may have gone through historical bottlenecks, and show ongoing genetic erosion (Moss et al. 2019).This may increase the risk of genetic drift and inbreeding for future Cyclura generations (Moss et al. 2019).

Conservation implications
The Species Management Plan for Cyclura nubila caymanensis (Binns et al. 2011) outlines the most urgent conservation measures to address the current threats to C. nubila caymanensis: habitat and nest-site protection, feral predator control, decreasing road mortalities, reducing the incursion of invasive species, and mitigating the impact of non-native predators.Additionally, methods to augment the populations, such as head-starting of hatchling iguanas, have been trialled on Little Cayman in recent years (J.Haakonson and T. Laaser, pers.comm.), employing some of the very successful methods of captive breeding and head-starting used for decades for Cyclura lewisi on Grand Cayman (Grant and Hudson 2015).A further potential conservation option, the translocation of iguanas between Cayman Brac and Little Cayman, has not been evaluated yet, as data on genetic diversity on and the genetic structure between the islands were missing.We are addressing this knowledge gap with this study.
Descriptive statistics, the mtDNA network, analyses of genetic diversity and genetic structure, as well as assessments of demographic history using both mtDNA and microsatellites all imply little to moderate genetic divergence between the Cayman Brac and Little Cayman populations of Cyclura nubila caymanensis.These genetic patterns support a model of some gene flow between the two island populations which would principally support translocations as a further management tool for the Lesser Caymans iguana, should these become needed in the future.However, moving animals from Little Cayman to Cayman Brac may result in the decline of the total Lesser Caymans iguana population, due to the larger anthropogenic pressure on the environment of Cayman Brac making it a population sink unsuitable for this option (Moss et al. 2019).This problem is exacerbated by the relatively low reproductive success rate of the species, given that the life history of Cyclura predicts relatively high fecundity but high juvenile mortality, which might limit the possibility of a quick population increase (Goetz and Burton 2012).Conversely, all Cayman Brac individuals could be relocated to Little Cayman, as this island has not yet experienced the same level of anthropogenic disturbance as Cayman Brac.However, this would leave the ecological niche that C. nubila caymanensis individuals currently fill on Cayman Brac unoccupied, with possibly detrimental consequences the biota of this island.
A potential augmentation of either population in terms of genetic variation might, if needed, most effectively be undertaken by the translocation of head-started juveniles taken from a proportion of hatching clutches.The potential of disease transfer must be considered in any translocation and any of the above scenarios; but especially in the face of a newly emerging bacterial disease identified in Cyclura lewisi on Grand Cayman (Conley et al. 2021).This causal bacterium -a novel Helicobacter species -is possibly transmitted through invasive common green iguanas (Iguana iguana), a species which has also formed non-native populations on Cayman Brac and continues to reach Little Cayman through supply transport links from Grand Cayman (Perry et al. 2021).Further, thorough population dynamics and viability analyses Lacy and Pollak (2022) should be undertaken in the planning stages to explore the consequences of various augmentation and translocation scenarios.All decisions should be guided quantitatively, and would ideally include mean kinship at a population or individual level (Ralls et al. 2017).Minimising mean kinship maximises the genetic diversity retention while conversely minimising inbreeding depression in subsequent generations (Ralls et al. 2017).
This study has highlighted the urgency to protect Cyclura nubila caymanensis through exploration of the genetic structure and diversity between two Cayman Island populations.While this study showed that translocation may be a possible option, it is important to note that the results and subsequent conservation implications are based on mtDNA sequences and only seven nuclear markers.Future work should include genomic analysis such as whole genome sequencing, which has successfully been undertaken on individuals of the critically endangered Lesser Antillean iguana (Iguana delicatissima) (Miller et al. 2019).This would provide a useful genomic resource, informing applied conservation and management for the critically endangered Sister Islands Rock Iguana.

Fig. 1
Fig. 1 Median joining haplotype network for the concatenated mtDNA dataset.Each circle represents a haplotype, and the size of the circle is proportional to the number of individuals belonging to that haplotype.Island populations are represented by different col-

Fig. 2
Fig. 2 Bar plots from Structure analysis for (A) K = 5 and (B) K = 2.The first 18 bars in the plots represent the genotypes of Cayman Brac individuals and the last 41 bars represent those of Little Cayman individuals.The proportion of each colour in the sample's genotype represents the proportion of the genotype belonging to a given genetic cluster

Fig. 3
Fig. 3 Plots showing different aspects of the demographic history of Cayman Islands iguanas.(A) Bayesian Skyline Plots based on mtDNA data show a very similar bottleneck scenario in both populations within the last 100,000 years.Middle lines represent median values, outer dashed lines represent upper and lower limits.Different colours represent the two different island populations.(B)

Table 1
Descriptive statistics of the mitochondrial DNA Summary statistics generated by DnaSP for cyt b and ND4 sequences for each iguana population.n sample size, PL DNA fragment length, NPS number of polymorphic sites, NH number of haplotypes, HD haplotype diversity, TW Watterson theta, π Tajima's nucleotide diversity.Chi-squared tests showed no significant differences in number of haplotypes between markers or populations, or number of polymorphic sites between markers or populations

Table 2
Summary statistics for the microsatellite dataF st pairwise fixation index, Fis inbreeding coefficient, Ho observed heterozygosity, He expected heterozygosity, ANAPL average number of alleles per locus, AR allelic richness, Nc census size (Goetz and Burton 2012), and Ne effective population size.T-tests showed no significant differences in Fis, H o , H e , ANAPL and AR between islands