Low‐coverage reduced representation sequencing reveals subtle within‐island genetic structure in Aldabra giant tortoises

Abstract Aldabrachelys gigantea (Aldabra giant tortoise) is one of only two giant tortoise species left in the world and survives as a single wild population of over 100,000 individuals on Aldabra Atoll, Seychelles. Despite this large current population size, the species faces an uncertain future because of its extremely restricted distribution range and high vulnerability to the projected consequences of climate change. Captive‐bred A. gigantea are increasingly used in rewilding programs across the region, where they are introduced to replace extinct giant tortoises in an attempt to functionally resurrect degraded island ecosystems. However, there has been little consideration of the current levels of genetic variation and differentiation within and among the islands on Aldabra. As previous microsatellite studies were inconclusive, we combined low‐coverage and double‐digest restriction‐associated DNA (ddRAD) sequencing to analyze samples from 33 tortoises (11 from each main island). Using 5426 variant sites within the tortoise genome, we detected patterns of within‐island population structure, but no differentiation between the islands. These unexpected results highlight the importance of using genome‐wide genetic markers to capture higher‐resolution genetic structure to inform future management plans, even in a seemingly panmictic population. We show that low‐coverage ddRAD sequencing provides an affordable alternative approach to conservation genomic projects of non‐model species with large genomes.


| INTRODUC TI ON
Many endangered species are restricted to a single or a small number of remnant populations. Management efforts often include introductions from these source populations to other suitable locations to lessen the risk of extinction or because the species in question are ecosystem engineers and can be used to restore degraded habitats elsewhere. However, such interventions have important implications for the genetic future of the newly founded population.
As only a subset of the individuals in the source population can be moved, genetic diversity is at risk to be lost and artificial population structure may be created in the new populations. Genetic diversity is essential for the adaptive potential of a species, particularly in the face of environmental changes and disease outbreaks (Reed, 2005;Reed & Frankham, 2003). Hence, management decisions need to be carefully planned to take the genetic characteristics of the source populations into account to aim at retaining as much genetic diversity as possible (Hoban et al., 2021).
One problem with assessing current genetic characteristics of endangered non-model species is that suitable marker systems, such as simple sets of microsatellites, are often unavailable. Nextgeneration sequencing provides promising tools at decreasing costs (Davey & Blaxter, 2010;Hayden, 2014). However, it can still be financially overwhelming and (if not outsourced) bioinformatically challenging to generate high-quality whole genomes, especially for species with large genomes, and because more than a handful of sequenced individuals are needed for population genomics studies (Corlett, 2017;Shafer et al., 2015). One potential solution is to use reduced representation sequencing, such as restriction-associated DNA (RAD) sequencing, which does not require a reference genome and is generally cost-effective (Andrews et al., 2016;Davey & Blaxter, 2010). Financial and computational costs of whole-genome sequencing of many individuals can be further reduced by adopting a low-depth sequencing strategy (Pasaniuc et al., 2012), where information on the whole genome is obtained, but at low coverage (generally 1-2×). This approach risks loss of genotype accuracy, which can be overcome by inferring genotype likelihoods (Fumagalli et al., 2014;Korneliussen et al., 2014). This genotype-free estimation of allele frequencies has been shown to reduce biases and improve demographic inference from RAD-seq data (Warmuth & Ellegren, 2019). Interestingly, to date, only a small number of studies have combined RAD and genotype-free estimation of allele frequency estimation approaches (Bay et al., 2019;Breusing et al., 2019;Peart et al., 2020;Záveská et al., 2019).
Here, we use low-coverage ddRAD sequencing as a time-and cost-effective approach for the population genetic analysis of Indian Ocean islands (Austin et al., 2003;Palkovacs et al., 2002) and occupies a prominent functional role in shaping and sustaining largescale vegetation dynamics as it is the largest frugivore and herbivore in its island ecosystem (Hansen, 2015;Hnatiuk et al., 1976;Merton et al., 1976). Therefore, A. gigantea are currently used to help restore degraded native ecosystems on several other Western Indian Ocean Islands (Griffiths et al., 2011Hansen et al., 2010).
Aldabrachelys gigantea was on the verge of extinction in the late 19th century due to excessive harvesting, with a population low in around 1870 of somewhere between <1000 and a few thousand tortoises (Bourn et al., 1999;Stoddart & Peake, 1979). Thanks to calls for protection from Charles Darwin and others in 1874, the number of A. gigantea increased quickly to several tens of thousands in the 1960s to today's stable population of well over 100,000 individuals (Turnbull et al., 2015). The Aldabra population is divided into several subpopulations across the different islands that make up the atoll.
Two previously published genetic studies of A. gigantea have involved samples from Aldabra's wild population and were all based on the mitochondrial control region or microsatellite data. The first, by Palkovacs et al. (2002), focused on captive individuals and examined potential genetic differentiation between morphotypes.
Although their sampling included some wild individuals, they did not examine the population structure within the atoll. The second study (Balmer et al., 2011) was based on samples from Malabar, Grande Terre South, and Grande Terre East. They found strong genetic differentiation between the two Grande Terre localities, and between Malabar and Grande Terre. They concluded that movements between different areas and islands are rare. Their sampling did not include samples from Picard and the study relied on eight microsatellite markers originally designed for Chelonoidis niger (split 35-40 mya, Kehlmaier et al., 2019). Using molecular markers developed for other species bears the risk of ascertainment bias and an underestimation of genetic variation (Delport et al., 2006;Ellegren et al., 1995). Similar problems have been encountered in F I G U R E 1 An Aldabra giant tortoise entering the Aldabra Lagoon other microsatellite studies of, for instance, turtles (Çilingir et al., 2019), fish (Carreras et al., 2017), and mammals (Hendricks et al., 2017;Mesnick et al., 2011).
Here, we provide a new sampling scheme for the first time including all the main islands hosting Aldabra giant tortoises and a new analysis that acts as a case study for the conservation genetic analysis of a non-model species using low-coverage sequencing combined with double-digest restriction site-associated DNA sequencing (ddRADseq; Peterson et al., 2012). Our specific aims are as follows: 1. To quantify the overall genetic structure of the endemic A. gigantea population 2. To determine whether there are significant differences in the genetic composition of the species among and within islands.

| Study system
The endemic distribution of A. gigantea is restricted to Aldabra Atoll, in the southern Seychelles. The atoll consists of four main islands, Grande Terre, Malabar, Polymnie, and Picard (Figure 2a), separated by channels and enclosing a shallow lagoon. On the atoll, giant tortoises are unevenly distributed across the three largest islands (Polymnie, the smallest main island, has no tortoises) due to environmental differences (e.g., terrain, food, freshwater resources, and shade availability) and differences in exploitation history (Bourn & Coe, 1978;Turnbull et al., 2015;Walton et al., 2019). Effective conservation management measures saving the species from extinction in the late 19th century included the reintroduction of tortoises to Picard and atoll-wide invasive species control (Bourn et al., 1999;Bunbury et al., 2018;Stoddart & Peake, 1979;Turnbull et al., 2015). The largest population lives on Grande   Table S1). Absolute ethanol was added to the blood samples in a 1:20 ratio to prevent coagulation (Wietlisbach, 2017). All samples were stored at room temperature until arrival in the lab and then at −80°C until DNA extraction.

| Sample collection and DNA extraction
DNA extraction was performed with 3 µl of blood (in ethanol) per sample, using the sbeadex™ kit (LGC Genomics, Middlesex, UK), following the manufacturer's protocol for DNA extraction from nucleated red blood cells. Genomic DNA concentrations were measured with a dsDNA Broad Range Assay kit (Qubit 2.0 Fluorometer, Invitrogen, Carlsbad).

| ddRAD-seq library preparation and sequencing
To keep sequencing costs as low as possible, we used a reduced representation genome sequencing approach, specifically the doubledigest restriction site-associated DNA sequencing (ddRAD-seq, Peterson et al., 2012). Restriction enzymes were selected based on in silico double-digest runs, using the SimRAD package within R v4.0.3 (Lepais & Weir, 2014;R Core Team, 2020 Quesada et al., 2019). We aimed for approximately 50,000 in silico RAD loci, which was achieved with the selected enzyme combination EcoRI-BfaI with a target size selection window of 300-350 bp (52,000 expected ddRAD loci, Figure S1).
We used 100 ng genomic DNA from each sample (n = 33) for the digestion. A single ddRAD-seq library was prepared by processing the 33 samples following the protocol by Peterson et al. (2012) with slight modifications as described in Çilingir et al. (2021). Briefly, after double digestion, the products were cleaned with a 1.0× ratio of AMPure XP beads. Next, the P1 adapters containing the inline barcodes unique to each sample (Peterson et al., 2012), and with an EcoRI overhang and the P2 adapter with a BfaI overhang, were ligated to the restricted DNA. Then, equal amounts of individually barcoded DNA were pooled. The double size selection was performed with a total of 300 μl pooled aliquot by treatment with 0.5× and 0.12× AMPure XP beads. After the size selection, eight PCR cycles were run using the common PCR1 and the PCR2 primers, which include a standard Illumina index (Peterson et al., 2012). In our case, only one index was used as there was only one sequencing library prepared. A final AMPure XP beads clean-up was followed

| Alignment to a reference genome, estimation of sequencing depth, and downsampling
After quality filtering, the paired reads were aligned to the C.
abingdonii reference genome using BWA-MEM version 0.7.17 ). Calculation of the average per site sequencing depth for each individual was done in three following steps. First, SAMtools ) was used to extract properly paired reads with mapping quality of >20 from the BAM file of individual GrdTr_11 (the individual with the highest number of sequence reads, Table S1). Next, for each individual, all positions with at  (Quinlan & Hall, 2010). Subsequently, per-site sequencing depth per individual was calculated using SAMtools  based on the range given by the bed file (all sites with at least 1× coverage).
Because the average sequencing depth per individual varied considerably, we repeated the major analyses after downsampling the forward and reverse Fastq files of each sample to equalize the number of reads per individual with seqtk v1.3 (https://github.com/ lh3/seqtk) to 154,599 reads (number of reads of individual Picard_2, third-lowest read count, Table S1).

| Estimation of genotype likelihoods
As the mean sequence coverage per sample was low (2.28×; range: 0.2-6.1×, Table S1), the uncertainty of genotypes was accounted for in the subsequent analyses by computing the genotype likelihoods at variant sites instead of calling genotypes. Accordingly, the read alignments of all 33 individuals were processed with ANGSD v0.93 (Korneliussen et al., 2014), a software developed for genomic analyses of low-coverage data. The GATK (Genome Analysis Toolkit) model was used (McKenna et al., 2010), and major and minor alleles were directly inferred from the genotype likelihoods (doMajorMinor 1, doMaf 1). Quality filtering for the subsequent downstream analyses was performed as follows: Only properly paired (only_ proper_pairs 1) and unique reads (uniquieOnly 1) were used, and only biallelic sites were retained (skipTrialleleic 1). Nucleotides with base qualities lower than 20 were discarded. Excess of SNPs around indels and excessive mismatches with the reference were corrected by realignment (C50, baq 1 [Li, 2011]). Reads with a mapping quality lower than 20 were discarded.
Additionally, for the estimation of genotype likelihoods, only SNPs with a p-value <10 −6 (the significance threshold for polymorphism detection) and heterozygosity <0.5 were retained, the latter to exclude potential paralogs (Hardy, 1908;Hohenlohe et al., 2011). Further filters were applied depending on the analysis. For the population genetic structure analyses, sites with read data in fewer than 30 of the 33 samples were excluded (minimum representation among samples >90%, -minInd 30). The minimum depth of sites to be retained was also set to 30, and hence, on average, at least one read per individual was required. The maximum depth per site was set as the sum of the average sequencing depth and two times the standard deviation (373 for the main dataset, 128 for the downsampled dataset). For the estimation of genetic differentiation and diversity, which were calculated per group, at least 50% of the samples in a particular group had to be represented (minInd = 50% of all individuals in a group). The minimum depth for each group was set to the minimum number of individuals allowed (50% of the overall individuals within a group) and the maximum depth was the average plus two times the standard deviation for each group.

| Estimation of kinship
To check for possible familial relationships potentially affecting the population structure analyses, the coefficient of kinship (Jacquard, 1974) was inferred by using NgsRelate v2 (Hanghøj et al., 2019). To achieve this, allele frequencies and genotype likelihoods estimated with the main dataset were used and average coefficients of kinship for all possible individual pairs were calculated.

| Population genetic structure
For a first overview of the population structure, a principal component analysis was carried out with PCAngsd v09.85 (Meisner & Albrechtsen, 2018) with an additional minor allele frequency (MAF) filter of 0.01 or 0.05. As a complementary population structure analysis, we used the clustering tool NGSAdmix (Meisner & Albrechtsen, 2018;Skotte et al., 2013). Similar to the Bayesian clustering method STRUCTURE (Pritchard et al., 2000), NGSAdmix allows the estimation of individual admixture proportions by assigning individuals to different clusters. While a PCA allows the assumption-free visualization of the genetic relatedness among individuals, NGSAdmix tries to minimize the within-group variation to define genetic groups and estimate individual admixture proportions (Meisner & Albrechtsen, 2018;Skotte et al., 2013). To use NGSAdmix, it is recommended to perform LD pruning (i.e., to filter sites based on pairwise linkage disequilibria) as the program assumes the independence of genomic loci (Skotte et al., 2013).
Hence, pairwise linkage disequilibria (LD) were calculated using ngsLD (Fox et al., 2019) and LD pruning was performed by allowing a maximum among SNP distance of 100 kilobases and a minimum weight of 0.5. A total of 100 replicates were performed for each NGSAdmix run and the number of clusters (k) varied between 2 and 10. The results were analyzed and visualized with CLUMPAK (Kopelman et al., 2015), and the log-likelihoods calculated for each run were visualized in R (R Core Team, 2020). Possible differences in genetic diversity among the five groups defined above, Malabar-1, Malabar-2, Grande Terre East, Grande Terre South & West, and Picard (see also Figure 2), were investigated by calculating average number of pairwise differences or nucleotide diversity (π; [Tajima, 1989]) and population mutation rate (Watterson's θ; [Watterson, 1975]). Both measures were based on SFS estimates and performed with the realSFS and ThetaStat modules in ANGSD (Korneliussen et al., 2014). Estimates of Watterson's θ and π were obtained per genome region via a sliding window analysis with a window and step size of 10 kilobases (non-overlapping windows, excluding windows with <10 sites). A Tukey's range test (David & Tukey, 1977) was applied to compare the diversity measures among different groups. Since Tukey's range test is a post hoc test, initially ANOVA was performed on the data.

| Genotype likelihood analysis
The sequencing effort yielded 23,517,270 raw reads for a total of 33 samples. After adapter removal, quality checking, and demultiplexing, an average of 1,188,685 (range: 129,188-2,610,656, see also

| Population genetic structure
We were primarily interested in the overall genetic structure and differentiation among islands. The PCA of the main dataset revealed two distinct clusters (PC1:7.8% and PC2: 4.26%, Figure 2b MAF > 0.05) had no effect on the PCA structure ( Figure S2A,B).
Next, we wanted to investigate if the observed population structure is consistent with two genetic groups as indicated by the PCA and we aimed at estimating admixture proportions. For this analysis, a total of 3781 LD-pruned SNPs with MAF >0.05 were used.
The admixture proportions indicated that when the number of putative clusters was assumed to be 2 (k=2 was the most likely number of k based on log likelihoods, Figure S3A), the M1 individuals were assigned to cluster A and the M2 individuals were assigned to cluster B, except for one individual of M1, which showed mixed ancestry ( Figure 2a). Also, Grande Terre showed high within-island differentiation with most of GE and one GW individual assigned to cluster B and most of GS and the remaining GW individual to cluster A. All P individuals except for two assigned to cluster B showed mixed ancestry. Hence, genetic groups of different islands were assigned to the same clusters (M1 with GS&W, and M2 with GE). This  (Table S1B) (Table S1B). Hence, there is no evidence for potential familial structure within sampling localities explaining the observed population structure.

| Estimation of genetic differentiation and summary statistics
All pairwise F ST estimates calculated among the three Aldabra Islands in the study (Picard, Malabar, and Grande Terre) were 0, suggesting no evidence for among-island differentiation. As expected from the PCA and the admixture proportion analysis, the major differentiation was found between M1 and GE (0.06), followed by GE and S&W (0.041) and M1 and M2 (0.039) (Figures   2f and 3). To account for possible local effects along the genome, we also compared within-and among-island differentiation of Although the absolute differences among groups were small, there was significant variation among the groups, F (4, 596387) = 75.23, p < 2e−16. Mean π values of all the groups differed from each other, except for GE, which did not differ from M2 or GS&W. While both Wattersons θ and π as well the analyses with the downsampled dataset suggested highest diversity in P, differences in diversity among the other groups were small (but significant) and the relative order of groups differed between analyses.

| DISCUSS ION
We used low-coverage RAD sequencing to investigate the population genetic structure and variation in the endemic A. gigantea population. Our data, although relying on a relatively small sample size (5-11 per island/sampling location), not only revealed the subtle genetic structure of previously bottlenecked populations but also suggested a potentially greater role of passive movement between islands via water in a terrestrial species than previously expected.
Our study is one of few to focus on a combination of reduced representation sequencing and the genotype likelihood approach to study the population genomics of a non-commercial and non-model species. Our case study supports the use of low-coverage ddRAD sequencing instead of the low-coverage whole-genome sequencing (Lou et al., 2021), which is still costly for large genomes and/or sample sizes. Sequencing costs depend on the platform, but could be as low as 5.25 USD per sample using our approach (2-3× coverage or 0.3Gb). In contrast, a low-coverage whole-genome resequencing project for a genome of about 2.4 Gb (the approximate genome size of A. gigantea) would result in sequencing costs of about 105 USD per sample (2-3× coverage or 6 Gb). Our reduced representation approach could therefore be particularly useful for species with very large genomes.
We also investigated the effects of unevenly distributed depth of sequencing per individual by repeating all analyses with a downsampled dataset. We showed that the results obtained with both datasets were consistent, but the downsampling led to a loss of resolution, especially for the admixture analysis. The smaller number of loci and among-locus variation in coverage known for RAD (Davey et al., 2013;O'Leary et al., 2018) may mean that there is a minimum acceptable depth of coverage for this technique.

| Unexpected partitioning of genetic structure
We found lower among-island than within-island differentiation.
Specifically, our analysis suggested two main groups of genetic variation ( Figure 2a): M2, and all but one individual from GE represented an eastern group; and M1 and all but one individual from GS&W represented a western group. The P individuals were assigned to both clusters, which was expected, given that the original population of Picard was extirpated in the 1800s, and the current population originates from reintroduced tortoises from Grande Terre and Malabar (Bourn et al., 1999). These findings were supported by the PCA and the pairwise F ST analyses, which showed minor differentiation between Picard and the other islands, but stronger differentiation within Grande Terre and Malabar. The genetic differentiation between GE and GS&W suggests that connectivity along the east-west axis of the island may be limited. This is in agreement with behavioral, ecological, and geographic observations (Bourn & Coe, 1978;Gibson & Hamilton, 1984;Swingland et al., 1989), and a previous study by Balmer et al. (2011). Areas of thick Pemphis scrub and deeply fissured rocks appear to limit the movement of tortoises F I G U R E 3 Heatmap of pairwise genetic differentiation (measured as F ST ), estimated for five different locations using the main dataset (Gibson & Hamilton, 1984). Hence, geographical barriers such as the region around Takamaka (dotted line in Figure 2a) that include deeply fissured limestone and thick Pemphis scrub together with isolation by distance probably explain the observed substructure on Grande Terre already described by Balmer et al. (2011) (see also Bourn & Coe, 1978).
More surprising and different from the previous findings of Balmer et al. (2011) was the low differentiation between M2 and GE. Occasional movement of tortoises carried by tidal currents from the mangrove area in Grande Terre East to the coastal area of M2 may cause inter-island gene flow (Figure 2a). The tortoises often spend days or weeks in the muddy mangroves of Grande Terre East.
Sometimes they move against tidal waters rushing out or in, with a risk of being swept away, and tortoises can even be spotted adrift in the open ocean outside the reef (Hansen et al., 2017). Ocean currents are increasingly acknowledged for their importance in shaping population structure (Arjona et al., 2020;White et al., 2010) and also for terrestrial reptiles (Calsbeek & Smith, 2003;Hawlitschek et al., 2017).
The movement of animals by humans could also explain the low differentiation. Although it is known that animals were transported from Grande Terre and Malabar to Picard for conservation purposes, there is no record of animals being transported from Grande Terre to Malabar or vice versa. It is therefore reasonable to assume that a direct route was taken for the tortoises en route to Picard, given that managing/transporting giant tortoises is a considerable effort. One potential reason for a lack of genetic differentiation between the GS&W and P is the movement of P tortoises to GW via the wide channel and islets between the two islands (Figure 2a), so that admixed P individuals influenced the GW group. However, this genetic similarity is likely to be driven by the founding history of the Picard population, which received individuals from Grande Terre (Bourn et al., 1999).
In summary, our findings suggest a subtle and unexpected signal of east-west population structure in A. gigantea, mainly correlated with landscape features, distance, as well as human-induced reintroductions (primarily on Picard). Seawater may play a less important F I G U R E 4 (a) Per-site estimates of Watterson's θ (b) and nucleotide diversity (π) obtained via a sliding window analysis performed with the main dataset; (c) per-site estimates of Watterson's θ (d) and nucleotide diversity (π) obtained via a sliding window analysis performed with the downsampled dataset. Each group is colored the same as in Figure 2 (orange, Grande Terre East; yellow, Grande Terre South & West; blue, Malabar Group 1; and green Malabar Group 2; pink, Picard) and the average value per each group is indicated with a black dot role as a barrier than has been previously assumed (Balmer et al., 2011;Grubb, 1971), instead water currents may support movements. Balmer et al. (2011) found no variation at the mitochondrial control region and there is currently no evidence for an ancient split into genetic groups. Given the very long generation time of giant tortoises, the substructure could still be several hundred years old and predate the species bottleneck.

| Limitations of the study
Our study provides the first genomic insight on the wild Aldabra giant tortoise population and the number of tortoises per each main island of the atoll included in this study was limited to 11. Population structure analyses have shown to be robust to extremely low (e.g., 0.125×; Lou et al., 2021) and highly uneven per-sample coverage (e.g., 0.5× to 6×; Skotte et al., 2013). But the sequencing effort (i.e., the combination of the number of samples and the sequencing depth per sample) affects the population genetic inferences obtained with genotype likelihood-based allele frequency estimations (Buerkle & Gompert, 2013). In a recent review including experimental design recommendations for different types of population genomic analyses using low-coverage whole-genome sequencing data, it was suggested to prioritize the total number of samples (≥10 samples per population) over per-sample coverage and to aim for ≥10× coverage per population both for population structure analyses (i.e., PCA and admixture analyses) and relative estimation of rare allele-dependent metrics (i.e., pairwise F ST and genetic diversity estimates). Our study design includes 11 samples (average coverage ~25×) per island (and originally expected population), and 5-6 samples (average coverage 11.4-13.7×) per sampling locality/genetic group. The unexpected outcome that there is more population structure than foreseen led to a rather low sample size per genetic group, while the coverage is still within the recommended range. Theoretically one would expect that these recommendations for whole-genome sequencing would apply for the reduced representation sequencing approach as well, given that the latter is a representative of the former. Our study does indeed show that it is possible to find subtle population structure with this kind of data. Nevertheless, we believe future population genomics studies of Aldabra giant tortoises would highly benefit from more samples per locality and per population as well as more extensive geographical sampling, for example, further western parts of Malabar and small islands in the Aldabra Lagoon.

| Conservation and research implications
Our study has several conservation implications. Specifically for the study species, our findings suggest that if giant tortoises from Aldabra are to be used for translocations, translocated individuals should ideally represent and potentially retain the overall genetic variation in the wild population. The local population on Picard was extinct and the current population is based on several bouts of turtle translocations since the early 1900s (Bourn et al., 1999) with the last occurring in the 1980s. Interestingly, individuals from Picard showed a mixture of the genetic assignments found on the other islands ( Figure 2), suggesting that the translocations likely involved more than one source. This is in accordance with the relatively high genetic diversity. Since the island also hosts the research station, individuals taken from Picard could be a valuable alternative and logistically more feasible than trying to capture individuals across the entire atoll.
However, our study does not yet give insights into the potential benefit of using captive or already rewilded populations for future translocations. Because A. gigantea were heavily exploited and exported to the outside of Aldabra Atoll in the 19th century (potentially before and during the species bottleneck) (Stoddart & Peake, 1979), it is not impossible that some of the original diversity now lost in the wild can still be found elsewhere. In any case, it is advised to translocate as many individuals as possible to minimize founder effects (Frankham et al., 2007). This should facilitate genetic management and monitoring of ongoing and future rewilding projects, including spatially larger projects in Madagascar (Pedrono et al., 2013), to maximize the evolutionary potential and survival of rewilded populations. As previously suggested by Balmer et al. (2011), we found no evidence for large differences in genetic diversity among the main islands and we currently do not see the need for translocations between islands.
The similar diversity among islands also suggests that the observed population structure is unlikely to be explained by the species bottleneck, for instance, by much stronger reduction and then isolation of the eastern part of Grande Terre. However, we caution that our method of low-coverage ddRAD has not been tested sufficiently for its reliability on the estimation of exact diversity measures (see recommendations from Lou et al., 2021).
Our study underlines the importance of genetically informed management decisions by showing unexpected population structure as previously discovered in Iberian wolves, Peruvian diving petrels, and Atlantic puffins, among others (Cristofari et al., 2019;Kersten et al., 2021;Silva et al., 2018). Considering that both reduced representation sequencing and low-coverage approaches aim to decrease costs, our approach could be used in the population genomics of other vertebrates to address similar research questions. Our approach could be particularly suitable for systems with large genomes and no or only little genetic knowledge is available when aiming at a first overall look at population structure. Finally, our study shows that land-and seascape genetics should go hand in hand because terrestrial organisms living close to the sea could be influenced by both.

ACK N OWLED G M ENTS
We would like to thank the Seychelles Islands Foundation and Aldabra research staff for their invaluable help in the field and their support in exporting the blood samples. The former Zurich-Aldabra Research Platform was involved in the blood sampling, and we especially thank Gabriela Schaepman-Strub for her support and contributions. We thank Xenia Wietlisbach for her help with sample management. We would also like to thank Claudia Michel for her contributions to the optimization of the ddRAD-seq protocol.
All the wet-lab work and the sequencing efforts required for this study were conducted in the Genetic Diversity Center, ETH Zürich.
During the course of this study, F. Gözde Çilingir was funded by the Georges and Antoine Claraz Foundation contributed to the sequencing costs of this project.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The sequencing data that support the findings of this study are avail-