Evidence of fine-scale genetic structure for reef manta rays Mobula alfredi in New Caledonia

: Our understanding of the genetic connectivity of manta ray populations and the drivers that shape genetic population structure is still limited. This information is crucial to identify the spatial boundaries of discrete populations and guide decisions on units to conserve. In this study, we used genome-wide single nucleotide polymorphisms (SNPs) to assess the genetic structure and diversity of reef manta rays Mobula alfredi at a local scale within New Caledonia and regionally in the western Pacific Ocean. We provide the first evidence of fine-scale genetic differentiation in M. alfredi , found between the 3 cleaning station aggregation sites in New Caledonia (n = 65 samples, N = 2676 SNPs, F ST = 0.01, p < 0.0001). Furthermore, population structure was evident at the regional scale between individuals from New Caledonia (n = 73) and East Australia (n = 19) on the basis of genetic differentiation statistics (3619 SNPs, F ST = 0.096, p < 0.0001) and clustering algorithms, with unidirectional gene flow detected from east (New Caledonia) to west (Australia). These results reveal that reef manta rays can form genetically distinct groups within a relatively small geographic range and highlights the need to consider genetic structure when designating management units for conservation action and planning.


INTRODUCTION
The metapopulation paradigm provides a useful framework to understand connectivity and population dynamics. In this context, local populations are viewed as discrete spatial entities that interact through gene flow and migration (Hanski & Gaggiotti 2004). Across species, the factors that influence the rate of gene flow between local populations are diverse and include dispersal or mobility, geographic distance between populations as well as life history traits (Allendorf et al. 2013). Environmental heterogeneity, and the resulting variation in habitat quality, is also a significant factor (Allendorf et al. 2013) as demonstrated by the impact of biogeographical barriers across different taxa, including insects (Rabasa et al. 2008), amphibians (Smith & Green 2005), birds (Esler 2000), fishes (Schtickzelle & Quinn 2007), and mammals (Matthiopoulos et al. 2005). Such biogeographical barriers are frequent in the terrestrial environment and significantly partition the landscape, creating a mosaic of ecosystems, such as forests, deserts, and grasslands (Briggs 1974).
One of the most challenging aspects of ecological studies is to identify population structure in apparently continuous ecosystems with no obvious physical barriers to dispersal, such as in the marine environment, and particularly for highly mobile species (Palumbi 1994). Over the past decade, studies have demonstrated that even when dispersal potential seems high, population structure can still occur in marine groups such as sharks (Geraghty et al. 2014, Vignaud et al. 2014), mammals (Andrews et al. 2010, and sea turtles (Roden et al. 2013, Read et al. 2015. In such cases, genetic population structure in the marine environment is shaped by an array of processes including oceanic conditions, currents, bathymetry, and geographical distances (Allendorf et al. 2013). Habitat choice, in particular site fidelity, is another key process that can drive population structure by limiting migration rates and encouraging assortative mating (Bowen et al. 2016). The main driver for site fidelity is often resource availability, as has been suggested for Hawaiian spinner dolphins (Andrews et al. 2010) and reef manta rays (Deakos et al. 2011, Perryman et al. 2019. To attain the key goal of conservation, i.e. protecting biodiversity by supporting the long-term persistence of viable, natural populations of wild species, requires the identification of population units and the geographic boundaries between them. This enables management and conservation programmes to appropriately focus and prioritise efforts to maximize evolutionary potential and minimize extinction risks (Allendorf et al. 2013). Whilst the presence of population structure may increase genetic differentiation through local adaptation, it may also lead to reproductive isolation and a reduction of genetic diversity (Futuyma 2019). Genetic isolation could diminish the resilience of a population or species and its capacity to adapt to a rapidly changing environment (Allendorf et al. 2013, Futuyma 2019. In addition, certain species are more vulnerable to extinction due to a combination of characteristics such as small population size (Pimm et al. 1988), limited geographic range (Gaston 1994), specialized habitat requirements (Brown 1995), and conservative life history traits (e.g. small litters, slow individual growth rate, late maturation, long inter-birth intervals, high survival rates; MacArthur & Wilson 1967). In the marine environment, these traits are found primarily, but not exclusively, in marine mammals and elasmobranchs.
The reef manta ray Mobula alfredi (Krefft, 1868) is a K-selected species and displays several traits that make it vulnerable in a changing environment (Mac -Arthur & Wilson 1967). These characteristics include large body size, long life expectancy, low fecundity , an often fragmented distribution (Couturier et al. 2011, Marshall et al. 2011, Rohner et al. 2013, and small local population sizes as estimated in East Australia ), Hawai'i (Deakos et al. 2011), Japan (Kashiwagi 2014), and Mozambique (Marshall et al. 2011). As planktivores, these animals depend on specific environmental conditions and processes that shape the abundance and distribution of their prey (Rohner et al. 2013. Food resource availability has been hypothesised to be the main driver of manta ray movements, often resulting in aggregations from a few to hundreds of individuals at feeding grounds and cleaning stations (Anderson et al. 2011, Couturier et al. 2011).
Patterns of seasonal or year-round aggregations have been found to differ with locations. On the one hand, seasonal migrations were correlated with monsoonal shifts and productivity in Indonesia (Dewar et al. 2008), the Maldives (Anderson et al. 2011, Kitchen-Wheeler et al. 2012, Harris et al. 2020, and East Australia  while large-scale movements occur along continuous continental coastlines (up to 1150 km, Armstrong et al. 2019) and between island chains (Germanov & Marshall 2014). On the other hand, physical barriers such as open expanses of sea and deep-water channels are thought to be a factor that may reduce the chance of reef manta rays transiting even between geographically close aggregation sites (e.g. 150 km, Deakos et al. 2011). Long-term residence patterns (spanning years to decades) have been recorded for reef manta rays at sites across the globe, including West Australia (Armstrong et al. 2020), East Australia (Couturier et al. 2011), Hawai'i (Deakos et al. 2011, Mozambique (Marshall et al. 2011), the Red Sea (Braun et al. 2015), British Indian Ocean Territory (Andrzejaczek et al. 2020), and Indonesia (Setyawan et al. 2018(Setyawan et al. , 2020. Such site fidelity could leave reef manta rays even more exposed to anthropogenic pressure if activities such as targeted fishing and bycatch, coastal development, or unmanaged tourism occur in critical habitats (Anderson et al. 2011, Couturier et al. 2011, Stewart et al. 2018. Listed as Vulnerable to Extinction on the IUCN Red List (Marshall et al. 2018), reef manta rays have been increasingly targeted by fisheries due to the high value of their gill-plates on the Asian market , O'Malley et al. 2017. Globally, the last IUCN assessment reports a suspected population reduction of 30−49% over the past 3 generations, with further reduction predicted over future generations (Marshall et al. 2018). More information is needed on population structure and the nature of the drivers that shape the movements and connectivity of manta rays to face these conservation challenges. Such information is essential to define management units and prioritise management efforts, for example the establishment of protected areas for genetically isolated groups.
Molecular tools provide opportunities for identifying species boundaries, patterns of gene flow, genetic diversity, and spatial structure in elasmobranchs (Ovenden et al. 2018). To date, studies on manta rays using molecular tools have focussed primarily on taxonomy , Hinojosa-Alvarez et al. 2016, Hosegood et al. 2020. In the first application of genomic tools to examine population structure in manta rays, Stewart et al. (2016) applied double-digest restriction siteassociated DNA (ddRAD) sequencing methods and found significant population structure in oceanic manta rays M. birostris at large geographic scales within the Indo-Pacific region. However, a more recent examination of a wider subset of M. birostris samples by Hosegood (2020) failed to find any population structure and suggested global panmixia for the species. This species is a close relative of M. alfredi and occupies mostly pelagic environments, demonstrating greater migration distances than reef manta rays (up to 1400 km, Hearn et al. 2014). In reef manta rays, genomic studies demonstrated population structure at large scales between the Indian and Pacific Oceans (Hosegood et al. 2020) and between regions separated by the Indian Ocean basin (Venables et al. 2021). To date, only 1 study has used genomic methods to investigate fine-scale population structure in reef manta rays and did not find significant genetic structure between reef manta rays sampled at aggregation sites separated by less than 100 km (within Mozambique) or less than 1000 km (between Mozambique and South Africa) (Venables et al. 2021). Molecular tools have also been used to assess structure in ecologically similar species, namely large filter-feeding elasmobranchs such as whale sharks (Vignaud et al. 2014) and basking sharks (Hoelzel et al. 2006), and revealed low levels of genetic differentiation among ocean basins.
Here we applied genomic methods to infer population structure and connectivity of reef manta rays at a local scale across New Caledonia, and compare populations regionally within the western Pacific Ocean to East Australia. The coastal waters of New Caledonia are listed as a UNESCO World Heritage Site, with 6 marine clusters representing a total of almost 16 000 km 2 of reefs, lagoons, and mangroves. Along with the Natural Park of the Coral Sea (1 292 967 km 2 ) these protected areas are part of an essential management and conservation process that integrates threatened species (GNC 2018).
In New Caledonia, reef manta rays are observed at different aggregation sites, used as cleaning stations or feeding grounds, along the outer slopes of the barrier reef of the archipelago (Fig. 1). These sites are either connected with continuous habitats, such as along coastlines and the barrier reef, or on isolated reefs and islands separated from other sites by deep open oceanic waters (> 2000 m). Preliminary observations using photo-identification (photo-ID) suggest that site fidelity is high but differs among sites, with some movement of individuals between sites (up to 330 km apart, H. Lassauce unpubl. data). In this study, we used single nucleotide polymorphisms (SNPs) from genotyping by sequencing (Kilian et al. 2012) to assess the genetic structure and diversity of the population of reef manta rays in New Caledonia. Specifically, we compared the genetic differentiation at a regional scale, with the inclusion of samples from East Australia, and at a local scale (levels of heterozygosity among sites in New Caledonia). In doing so, we provide the first estimates of genetic structure and diversity for this species in New Caledonia.

Study sites
This study was conducted at 4 sites around New Caledonia, an archipelago principally consisting of a larger Main Island and 3 smaller islands approximatively 100 km off the east coast known as the Loyalty Islands (Fig. 1). The Main Island is surrounded by a 1600 km barrier reef that shelters the shallow waters of a 16 874 km 2 lagoon (Andréfouët et al. 2009). The continental shelf ends at the barrier reefs, where it drops to depths greater than 2000 m. Bathymetry around the Loyalty Islands is similar, with a relatively narrow continental shelf and a deep channel (> 2000 m) between the Loyalty Islands and the mainland (Fig. 1).
Reef manta rays were sampled at 3 sites off the Main Island: Noumea in the south and Pouembout in the north are located on the west coast, while Touho is in the northern part of the east coast. One site, Ouvea, is located on the northern island of the Loyalty Islands (Fig. 1). At all sites, except Pouembout, manta rays are observed primarily on the reef slope (between 10 and 15 m deep) near reef channels, at cleaning stations. In Pouembout, samples were collected opportunistically from an aggregation of manta rays that were found at the surface in shallow waters (<12 m), feeding in nutrient-rich waters of the lagoon.
Additional specimens from 2 sites in Queensland, Australia (hereafter referred to as 'East Australia', see Fig. 1), were included in the study to serve as a reference group when determining the genetic population structure within New Caledonia.

Sampling methods
Tissue samples were collected from the pectoral fin of reef manta rays using a biopsy tip attached to a 2 m spear pole (Pneudart; Stewart et al. 2016). Sampling was conducted via either SCUBA diving or free diving, and samples were stored in 95−100% ethanol at −20°C. All individuals sampled were photo-identified and designated as either adult males (based on the presence of fully developed claspers), adult females (based on the presence of mating scars or pregnancy), or juveniles (defined by the absence of fully developed claspers or mating scars). The exception to this were manta rays sampled from Pouembout, where the high turbidity of the water did not allow for photo-identification or confident observations of sex. Sampling information is summarized in Table S1 in the Supplement at www.int-res.com/ articles/suppl/n047p249_supp.pdf.

Laboratory procedures
Genomic DNA was extracted from samples using a Qiagen DNAeasy Blood & Tissue Kit following the manufacturer's instructions. DNA was quantified using the Qubit dsDNA Broad range assay (Thermo Fisher Scientific) and then diluted with Qiagen AE buffer to 50 ng μl −1 for use in genotyping by sequencing.
Development and genotyping of SNPs for all samples was undertaken following the standard DArTseq protocol (Kilian et al. 2012). DArT (Di versity Arrays Technology) is a genotyping-by-sequencing platform that uses an enzymatic double-digest genome complexity reduction following San saloni et al. (2011) and next-generation sequencing on the Illumina HiSeq 2500 platform as described by . In this study, the restriction enzymes PstI and SphI were used for the double digest of genomic DNA.

Data quality check and cleaning
Sequencing reads were processed using proprietary DArT analytical pipelines (Jaccoud et al. 2001, Kilian et al. 2012). In the primary pipeline, poor-quality sequences are filtered, applying more stringent selection criteria to the barcode region in comparison to the rest of the sequence. Identical sequences were then compiled into .fastq files and cleaned through the secondary pipeline using the DArT algorithm (DArTSoft14, Kilian et al. 2012). This pipeline workflow first calls sequence clusters for all pooled samples and then for each individual. The DArT pipeline retained SNPs on the basis of the balance of average count for each SNP allele, the reproducibility values (> 95%), the average count for each se quence or row sum (sequencing depth: ≥ 5), and the call rate (proportion of samples for which the marker is scored: > 95%). Potential contamination was identified by comparing all reads to bacterial and viral sequences from GenBank and a custom DArT database. Lastly, we applied the following filters to generate the final dataset: < 5% of individuals with missing data (call rate of > 95%), average reproducibility of alleles at a locus (> 95%), across SNPs that share a sequence tag, only 1 random SNP was retained per locus, and monomorphic loci were removed. Individuals with a call rate < 80% were also removed. After this filtering process, average heterozygosity per locus per individual was plotted and visually inspected for potential outliers indicative of sample contamination. The fine-scale analysis of structure focussed on the cleaning station sites, to the exclusion of Pouembout. This is because the relationship between cleaning and feeding aggregations is unclear, and feeding aggregations could contain manta rays from multiple cleaning stations and vice versa, decreasing any population genetic signals. Two datasets were generated, one including the outgroup (East Australia), referred to as the regional dataset, and the other excluding this outgroup, called the local dataset (New Caledonia only). For the analysis of genetic diversity only, loci with minor allele frequencies > 0.05 were filtered out of these 2 datasets. All filters were applied using the 'gl.filter' functions in the R package 'dartR' ).

Genetic diversity
We calculated genetic diversity for each sampling site in New Caledonia, for New Caledonia as a whole, and for East Australia using the statistical programming language R version 3.6.2 (R Core Team 2019). Specifically, we estimated observed (H o ) and expected (H e ) heterozygosity (Nei 1987) with the function 'gl.basic.stat' in the R package 'dartR' ). The inbreeding coefficient (F IS ) (Weir & Cockerham 1984) for each sample partition was also calculated using the 'gl.basic.stat' function. The rarefacted allelic richness (A R ) was estimated for each sample partition using the function 'allelic.richness' in the 'adegenet' package in R (Jombart 2008). Since sample partitions have different sample sizes, the function used the minimum number of individuals (n = 19 for the regional dataset and n = 18 for the local dataset). Differences in the diversity between sample partitions were tested using ANOVA and ttests (equal variance) or Welch's t-tests (unequal variance). Levene's tests were used to test the homogeneity of variances, and a Shapiro-Wilk test assessed the normality of the data.

Population differentiation and structure
We inferred genetic differentiation and population structure at a regional scale between East Australia and New Caledonia and then at a local scale among our study sites. We also compared the genetic diversity at both regional and local scales.
We conducted genetic clustering using 2 complementary methods: discriminant analysis of principal components (DAPC) (Jombart et al. 2010) and the program TESS3 (Caye et al. 2016). The outcome of the DAPC depends on the number of principal components (PCs) used in the analysis and the number of clusters in the analysis. To assess the number of PCs to be used for each dataset, we examined the mean alpha score (difference in reassignment probabilities between observed and permuted clusters), using the function 'a.score' and 'optim.a.score' of the R package 'adegenet' (Jombart 2008). Alpha scores were calculated using the highest recommended number of PCs (N/3, were N is the total number of individuals) as recommended by the developer of 'adegenet' (Jombart 2008). The DAPC was then performed with the 'dapc' function in the R package 'adegenet' (Jombart 2008) without assigning individuals to clusters a priori. This method derives a number of potential populations from the analysis and produces membership probabilities for each sample based on the data from all samples. Classification of a given sample to a cluster is based on the highest assignment probability. The expected classification rate for each population (prior) was compared to the number of correct classifications (posterior) generated using the DAPC.
Population genetic structure was further assessed using Bayesian clustering analysis of individual geno types with the 'TESS3R' package (Caye et al. 2016) in R to estimate individual ancestry coefficients assuming admixture of K ancestral genetic clusters in our data set. TESS3 implements a spatially explicit Bay esian clustering admixture model. By including information on individual geographic coordinates, TESS tends to perform better than non spatial models implemented in programs using a similar algorithm, such as STRUCTURE (Pritchard et al. 2000), when the ancestral level of differentiation is low (Durand et al. 2009, François & Durand 2010. It also performs better when the dis-tribution of the studied species tends to be continuous, with individuals geographically close and more prone to share ancestral genotypes (Durand et al. 2009). In the latter case, individuals are assigned to the most likely cluster on the basis of multi-locus genotypes sampled at distinct geographical locations without assuming predefined populations (Du rand et al. 2009). Here, the specified spatial prior is weak compared to the amount of information contained in the molecular dataset, and therefore, deviations from spatial smoothness are still al lowed in the posterior inference (Corander et al. 2008). We added UTM coordinates for each individual and ran 20 replicates of the admixture model for each value of the maximal number of clusters (K). The best value of K was examined from 1 to 10 with a maximum number of 10 000 iterations per run (20 repetitions) and a tolerance value of 10 −7 . A cross-validation procedure with 10% of masked data was used to select the best value of K according to the asymptote in the plot of cross-validation scores (Caye et al. 2016). Given that fine-scale population structure has not previously been described for this species, we did not have an a priori expectation of multiple genetic clusters within New Caledonia, and performed the analysis at both the local and regional scale. The analysis was processed for K-values from 2 to 5, although only the results for K = 2 and K= 4 are presented, as only these clustering results warrant a biological interpretation.

Gene flow analysis
Relative directional migration rates were analysed using the 'divMigrate' function in the 'diveRsity' R package (Sundqvist et al. 2016). This function gives relative migration rates proportionally to the highest migration rate between 2 sites, which always equals 1. We compared 3 differentiation metrics (see Section 3.3.4): 'd' (Jost's D), 'gst' (Nei's G ST ), and 'Nm' (Alcala et al. 2014). Migration rates were estimated with 10 000 bootstraps, and directional gene flow estimates with a rate < 0.2 were filtered out of the visual outcome. We used samples from New Caledonia (Noumea, Touho, and Ouvea) and from East Australia for this analysis.

SNP dataset
The regional dataset initially comprised a total of 5124 bi-allelic SNP loci in 73 reef manta rays from 4 locations in New Caledonia and 20 individuals from East Australia (Table 1). One individual (MA138) from East Australia was removed due to a call rate of less than 80%. The local dataset had a total of 3973 bi-allelic SNP loci in 65 reef manta rays from 3 locations in New Caledonia. One of the sites (Pouembout) was removed from the finescale analysis due to the small sample size (n = 8) and difference in habitat type compared with the rest of the sites. After the filtering, the total number of polymorphic SNP loci was 3619 for 92 individuals in the regional dataset and 2676 for 65 individuals in the local dataset. All filtering steps are sum marized in Table 1. The mean and distribution of heterozygosity was similar among individuals within sample sites and regions with no outliers (Fig. S1)

Genetic diversity
Genetic diversity statistics for all sample partitions are shown in Table 2. H e was significantly higher for the East Australian outgroup (t 3978 = 6.35, p < 0.001) compared with the overall New Caledonian dataset while no significant difference was observed for A R (t 3978 = 0.76, p > 0.05) ( Table 2). When considering the local dataset, there were no significant differences in any genetic diversity statistics between sample partitions within New Caledonia (F = 0.05, p > 0.5, Table 2).
Inbreeding coefficients (F IS ) were relatively high at all levels of population structure, and all were significantly different from zero (p < 0.001) ( Table 2). At the sample site level, the average proportions of heterozygosity in individuals were significantly different (t 3816 = 2.18, p < 0.05) between samples from East Australia (mean ± SD, F IS = 0.0662 ± 0.03020) and individuals from all groups of New Caledonia in the regional dataset (F IS = 0.0844 ± 0.02093). Among groups from New Caledonia, values ranged from 0.0699 ± 0.02762 in Noumea to 0.0812 ± 0.02918 in Ouvea, with no significant differences among sites in the local dataset (F = 0.93, p > 0.05).

Population differentiation and structure
3.3.1. Pairwise F ST and Nei's genetic distance D Statistically significant genetic differentiation was observed between all sampling sites ( Table 3). The highest estimates of genetic differentiation were found between East Australia and New Caledonia (F ST = 0.0958 ± 0.006 SD, p < 0.001, Table 3). The degree of differentiation between New Caledonian sampling sites was lower but still statistically significant (pairwise F ST = 0.0116−0.0165, Table 3). Overall F ST values for the regional and local datasets are provided in Table 3. Nei's genetic distances followed the same pattern as the fixation indexes.  Table 3. Pairwise F ST values (above diagonal, with 95% CI in parentheses) and Nei's genetic distance (below diagonal) calculated for reef manta rays Mobula alfredi for regional and local datasets using 3619 and 2676 SNPs, respectively. n: sample size. All F ST values are significant at p < 0.001

DAPC
Examination of the mean α-score values revealed that the appropriate number of PCs was 30 for the regional dataset and 21 for the local dataset (Tables S2 & S3).
The DAPC showed individuals clustered into their country and sampling site of origin with high confidence (Tables S4 & S5). Considering the regional data set first, DAPC revealed a clear distinction across all samples from East Australia and New Caledonia, with PC1 providing clear discrimination (Fig. 2). Using the classification algorithm, all individuals had a 100% probability of being assigned back to their country of origin (Table S5).
Considering the local dataset, the first of the two DA eigenvalues represented most of the variation within this dataset (Fig. 3). The DAPC discriminated sites along both the north−south (PC1, discriminating Noumea) and east−west axes (PC2, differentiating Ouvea and Touho). Within New Ca le donia, the highest assignment probabilities correspond to the population of origin for each pre-defined group (Table S5, Fig. S2). The average correct classification probability was highest for a sample from Noumea (mean ± SD, 94.12 ± 13.84), then Touho (85.43 ± 27.23), followed by Ouvea (79.79 ± 30.71) with only 3 (12%) individuals per site likely to be incorrectly classified (Table S5, Fig. S2).

Estimation of spatial structure: TESS3 analysis
For the regional dataset, the TESS3 analysis revealed the most likely number of ancestral populations/genetic clusters or best value of K as 4 (the cross-validation test exhibited a minimum value or a significant plateau at K = 4) (Fig. S3), supporting subtle levels of population structure (Caye et al. 2016). At K = 2, a strong distinction is revealed between Australian individuals and all individuals from New Caledonia, with distinct admixture compositions with 2−5 and 0−9% ancestry proportions of each other's dominant cluster, respectively. At K = 4, each site has its own distinct admixture composition, with a dominant ancestral population (Fig. 4). On average, New Caledonia individuals have less than 1.4 ± 0.003% ancestry proportions of the dominant East Australia cluster. The dominant ancestral populations comprise, on average, between 76 ± 0.006% (Ouvea) and 84 ± 0.008% (East Australia) of ancestry proportion per individual.

Gene flow analysis
Relative directional migration rates were estimated using 3 differentiation metrics, but, as the patterns were consistent, only Jost's D results are presented (Fig. 5). The values displayed in the figures do not represent the magnitude of each mi gration (as a 256 Fig. 2. Inference of population structure of reef manta ray Mobula alfredi from East Australia and New Caledonia (regional dataset) using 3619 single nucleotide polymorphisms: plot of individuals based on the first 2 principal components (PCs) of the discriminant analysis of principal components (DAPC). The inset shows the eigenvalues of the analysis, representing the proportion of genetic structure captured by each of the PCs. Black PCs were used in the DAPC whereas greyed out ones were not measurement of the proportion of migrants per generation would) but compare migration rates proportionally to the highest, which always equals one.
The relative migration rates between East Australia and all of the New Caledonian sites are low but show statistically significantly higher migration rates from New Caledonia to East Australia than vice versa. Considering the local dataset, the highest bi-directional migration rate was found between Ouvea and Touho. Comparatively, the relative migration rates between Noumea and both Ouvea and Touho were slightly lower but also bi-directional.

DISCUSSION
Previous findings have indicated ocean-wide genetic differentiation and limited regional genetic connectivity in mobulid rays (Stewart et al. 2016, Hosegood et al. 2020, Venables et al. 2021). Here we build on this work to provide evidence of fine-scale population genetic structure of reef manta rays, highlighting new considerations for the management of the species. Despite many calls for studies such as ours to provide information on units to conserve and to identify drivers of population structure in mobulids, only a few studies have investigated the population genetics of this taxon (Deakos et al. 2011, Stewart et al. 2016, 2018, Setyawan et al. 2018, Perryman et al. 2019, Hosegood et al. 2020). The only comparable study that investigated genetic structure of reef manta rays at a similar spatial scale (within hundreds of kilometres) and type of DNA marker (3057 SNPs) was conducted within southern Mozambique and found no evidence of genetic differentiation (Venables et al. 2021). Our findings provide evidence that fine-scale population genetic structure can exist for the species and may depend on its habitat use and the geographical context, highlighting the need to be assessed on a case-by-case basis.

Regional scale population genetic structure
We found strong genetic differentiation at the regional scale, which may indicate that large expanses of water lead to increased genetic differentiation in reef manta rays. This statistically significant population structure at a regional level is consistent with previous studies on the reef manta ray and its close relatives. For example, 2 studies reported clear genetic distinction at an ocean basin level (> 7000 km) for populations of Mobula alfredi between the Indian and Pacific Oceans (F ST = 0.16, Hosegood et al. 2020, F ST = 0.38, Venables et al. 2021. In comparison to our findings, the larger magnitude of F ST in these papers may be due to the greater geographic distance be tween examined populations, which can be an important driver of genetic differentiation as in other taxa (e.g. sharks, Geraghty et al. 2014, Vignaud et al. 2014; humpback and southern right whales, Rosenbaum et al. 2017, Carroll et al. 2020b). However, further in vestigations would be required to validate this statement. While Venables et al. (2021) used a similar SNP discovery process to our study, a comparison with the study by Hosegood et al. (2020) might be biased due to the use of a different procedure. The use of common methods and sets of genetic markers would allow for comparison at a broader geographic extent and enable a better understanding of the genetic structure and diversity of the species (Domin gues et al. 2018). Stewart et al. (2016) discovered substantially lower levels of genetic population structure for the oceanic manta ray (F ST < 0.004) compared with the reef manta ray (F ST = 0.02) across the Indian and Pacific Oceans, using 3108 SNPs from ddRAD sequencing, but employing different discovery and filtering processes to those em ployed here. We propose that the interspecific difference in population structure found between oceanic and reef manta rays is largely due to different biology and ecology, with some variation also due to the SNP genotyping procedure. Reef manta rays under take smaller migrations, have a smaller home range, and demonstrate higher site fidelity than oceanic manta rays . Oceanic manta rays have been observed to undergo long-distance travels be tween sighting locations (up to 1500 km, Hearn et al. 2014), and vagrants are occasionally re corded outside their known distribution (Couturier et al. 2015). Strong genetic differentiation has also been found in populations of other reef-associated elasmobranchs separated by vast expanses of water (Vignaud et al. 2014, Momigliano et al. 2017. Physical barriers such as 258 Fig. 5. Population network and relative migration rates of reef manta rays Mobula alfredi based on Jost's D estimates of genetic differentiation for samples from East Australia and New Caledonia for the combined regional (East Australia and New Caledonia) and the local (only New Caledonia) datasets. The thickness of connecting lines is proportional to the relative rate of migration. Abbreviations are as follows; A: East Australia; N: Noumea; O: Ouvea; T: Touho. Statistically significant asymmetrical rates are marked with an asterisk, and migrations with a relative rate < 0.2 are not displayed open expanses of sea and/or deep-water channels are thought to re duce the chance of reef manta rays (Deakos et al. 2011) and other reef-associated species (Heupel et al. 2019) to transit even between geographically close aggregation sites (e.g. 150 km, Deakos et al. 2011, Setyawan et al. 2018). Our analysis found that the admixture composition of individuals from East Australia showed a greater proportion of ancestral genetic clusters most common in New Caledonia than the reverse, and migration rates were higher from New Caledonia toward East Australia. This difference seems to indicate that East Australia is a recipient of gene flow from New Caledonia. Previous studies on reef manta rays along the east coast of Australia have reported relatively high connectivity between sites separated by 380 km, with a high proportion of individuals re-sighted at more than one site using photo-identification (Couturier et al. 2011 and detected with satellite ) and acoustic tagging (Couturier et al. 2018). In New Caledonia, satellite tagging does not indicate any offshore movements extending beyond New Caledonian waters (H. Lassauce unpubl. data), supporting the finding of genetic isolation by suggesting that migration is a rare event. Further research is needed on the genetic structure of other populations of reef manta rays in adjacent countries such as Vanuatu, Fiji, or the Solomon Islands to extend our understanding of the species' regional connectivity.

Local-scale population genetic structure
The present study is the first to report fine-scale population structure for M. alfredi. In New Caledonia, genetic differentiation was detected in all of our analyses between the 3 study sites: pairwise F ST values were low but significant between all sites; DAPC distinguished 3 separate clusters and revealed high assignment probability rates for most individuals to their sampling location; and the TESS3 analysis revealed distinctive admixture compositions between samples from each site. Population structure within such a relatively small area (maximum oceanic distance between sites is 335 km) has not been recorded for this species, and was not an a priori hypothesis going into this study. However, the low level of differentiation indicates gene flow between the study sites. Therefore, we suggest that the New Caledonian reef manta ray exists as a metapopulation with at least 3 distinct local populations connected by gene flow.
Our understanding of these results is limited by the uncertainty in the spatial and reproductive ecology of reef manta rays, which appears to vary between regions of the world (Stewart et al. 2018). On the one hand, large-scale movements reported from multiple locations suggest that populations of reef manta rays occupy large areas that include several key aggregation sites (Anderson et al. 2011, Marshall et al. 2011, Germanov & Marshall 2014, Armstrong et al. 2019. On the other hand, reef manta rays have also been recorded to exhibit high site fidelity with sometimes clear segregation between geographically close aggregation sites around the world (Dewar et al. 2008, Deakos et al. 2011, Marshall et al. 2011, Setyawan et al. 2018, Perryman et al. 2019. For instance, in Hawai'i, Clark (2010) and Deakos et al. (2011) reported strong evidence through photo-ID of longterm, high site fidelity and no connection between 2 known aggregation sites less than 150 km apart. Similarly, in Indonesia, Setyawan et al. (2018) found high levels of site fidelity with acoustic telemetry and the absence of connection between study sites located only 150 km apart. High site fidelity has been documented in other regions of the world, such as in Mozambique (Marshall et al. 2011) and Indonesia (Setyawan et al. 2020). The habitat types or quality that could drive such differences in behaviour are not well understood. In New Caledonia, ongoing studies involving photo-identification and satellite telemetry seem to indicate high site fidelity and only a few connections between all 3 study sites (H. Lassauce unpubl. data). Only a few migrants per generation are needed to reduce genetic differentiation (Wang 2004). Therefore, it could be that long-term residency with limited movement is the norm, with rare migration events linking local regions with gene flow.
Our results did not suggest that the deep ocean was a barrier to the dispersal of manta rays within New Caledonia as hypothesised by Deakos et al. (2011) in Hawai'i for M. alfredi, where no connection was discovered between 2 islands (150 km apart) separated by a 2000 m deep channel. In the context of our study, the continuous coastal environment, specifically the relatively shallow water between Noumea and Touho, could favour genetic connectivity (Couturier et al. 2011 while the deep-water channel separating these sites from Ouvea could act as a barrier discouraging gene flow (Deakos et al. 2011). In New Caledonia, deep diving behaviour of reef manta rays has been recorded using satellite telemetry (Lassauce et al. 2020). One individual trav-elled between Touho and Ouvea via a 2000 m deep channel, and several animals made offshore foraging excursions (H. Lassauce unpubl. data), behaviours consistent with research on other manta ray populations (Germanov & Marshall 2014). These observations seem to indicate that deep waters may not be a barrier across a small spatial scale and that, in fact, other factors, possibly associated with the ecology and behaviour of the species, may drive fine-scale genetic differentiation.
Explanations for population structure include social behaviour and/or distribution of food resources. In New Caledonia, all of our study sites (except Pouembout) are described as cleaning stations. Cleaning stations are critical habitats for manta rays. These sites are gathering points that enable social interactions between individuals and have been identified as essential areas for reproduction (Stevens et al. 2018) and other social behaviour (Perryman et al. 2019). The reproductive ecology of reef manta rays is still unclear, but some evidence suggests that females reside longer in areas with high mating potential and sufficient food resources while males move between aggregation sites (Marshall & Bennett 2010, Deakos et al. 2011, Stevens 2016, Stevens et al. 2018. For instance, the seasonal migration of M. alfredi in the Maldives was influenced by monsoon currents that promote phytoplankton blooms (Anderson et al. 2011). Migrations associated with food availability have also been documented in Indonesia (Dewar et al. 2008) and East Australia . It is possible that high productivity attracts manta rays and variations in conditions provoke migrations (Deakos et al. 2011, Setyawan et al. 2018, and so further work is needed to investigate correlations be tween migration events and prey availability.

Genetic diversity
Genetic diversity was significantly higher in East Australia (statistically significantly higher H e and lower F IS , Table 2) compared with New Caledonia, consistent with the larger population size of the former  or its role as a recipient of gene flow. Indices of genetic diversity were similar among sites in New Caledonia (Table 2).
Inbreeding coefficients (F IS ) revealed significantly higher values in New Caledonia (F IS = 0.084) than in East Australia (F IS = 0.066), both of which differ from the F IS values not significantly different from zero presented by Venables et al. (2021) for the species using the same standard DArTseq protocol. Another study using a similar SNP-calling procedure found similar F IS values (ranging from 0.065 to 0.070) for a highly migratory species of shark (school shark Galeorhinus galeus) with no genetic structure at a regional level (Devloo-Delva et al. 2019). In contrast, Glaus et al. (2020) revealed lower inbreeding estimates for bull shark Carcharhinus leucas populations from widespread locations across the Indian and Pacific Oceans. Additional work that employs next-generation sequencing is necessary to assess genetic diversity in other populations of reef manta ray. This would allow for a broader comparison and a better understanding of the genetic diversity of the species and its drivers at a global scale to achieve effective conservation and management.

Management and conservation
The challenge for conservation and management of mobile species is to identify the relevant management units. Information on the extent to which populations are genetically subdivided is crucial to establish effective conservation measures (Palsbøll et al. 2007). The findings presented here suggest that (1) the New Caledonian manta ray is a distinct population, and potentially a separate evolutionary significant unit (Waples & Gaggiotti 2006) from East Australia, and (2) New Caledonian manta rays exist in a metapopulation with sites or sub-populations within the reef system linked by gene flow. This work contributes to the use of genetic tools to identify ecological units for the creation of appropriate conservation measures rather than geographically or politically based legislation (Deakos et al. 2011, Stewart et al. 2016, 2018, Setyawan et al. 2018, Perryman et al. 2019, Hosegood et al. 2020. Future genetic monitoring should be used to detect changes in connectivity and diversity (Leroy et al. 2018), which can presage declines across the metapopulation (e.g. Carroll et al. 2020a). Manta rays currently have no legal protection in New Caledonian waters despite the fact that populations of long-lived, slow-reproducing species are thought to have poor resilience to even low amounts of mortality (Stevens et al. 2018). Therefore, a cautionary approach would be to combine genetic sampling as part of a broader long-term monitoring project that would provide estimates of important population parameters over time (Clutton-Brock & Sheldon 2010). Such a monitoring regime should also include the monitoring of any bycatch and population viability analyses to assess the impact of possible threats.
Acknowledgements. We thank Conservation International, the Aquarium des Lagons, and The Manta Trust for their overall support. We especially thank Franck Bouilleret, Hugues Gossuin, and Thomas Auger for their invaluable assistance in field operations, as well as Rochelle Constantine, Olivier Chateau, and François Tron for their help and support. We also thank the Hô-üt Association and all of its members, and all those involved in the field operations for their valuable help (Ludovic Mazens, Marie-Renée Pabouty, Amaury Durbano, Jean-Baptiste Badiou, Taiki Kido, and the team of divers from the Aquarium des Lagons). We finally gratefully acknowledge the spokespersons for customary authorities of Ouvéa and Touho, and the 3 Provinces for allowing such operations. H.L. was supported by a scholarship from the Southern Province of New Caledonia, E.L.C. was supported by a Rutherford Discovery Fellowship from the Royal Society of New Zealand Te Apārangi. Laboratory and field work were funded by the Keidanren Nature Conservation Fund (KNCF) through the SATO YAMA UMI Project and Conservation International. Australian sample collections were funded through the Australian Research Council Linkage Grant: LP150100669. Sample collection was conducted with authorizations from the Southern Province (permit no: 34584) and the Northern Province (permit no: 609011-33) of New Caledonia. In the Loyalty Islands Province, no permit was required by the local authorities, although permission was granted by the local spokespersons for customary authorities. Transportation of tissue samples from New Caledonia to New Zealand was conducted with the authorization of the CITES authorities under permit no.