Genomic evidence indicates small island-resident populations and sex-biased behaviors of Hawaiian reef Manta Rays
BMC Ecology and Evolution volume 23, Article number: 31 (2023)
Reef manta rays (Mobula alfredi) are globally distributed in tropical and subtropical seas. Their life history traits (slow growth, late maturity, low reproductive output) make them vulnerable to perturbations and therefore require informed management strategies. Previous studies have reported wide-spread genetic connectivity along continental shelves suggesting high gene flow along continuous habitats spanning hundreds of kilometers. However, in the Hawaiian Islands, tagging and photo-identification evidence suggest island populations are isolated despite proximity, a hypothesis that has not yet been evaluated with genetic data.
This island-resident hypothesis was tested by analyzing whole mitogenome haplotypes and 2048 nuclear single nucleotide polymorphisms (SNPs) between M. alfredi (n = 38) on Hawaiʻi Island and Maui Nui (the 4-island complex of Maui, Molokaʻi, Lānaʻi and Kahoʻolawe). Strong divergence in the mitogenome (ΦST = 0.488) relative to nuclear genome-wide SNPs (neutral FST = 0.003; outlier FST = 0.186), and clustering of mitochondrial haplotypes among islands provides robust evidence that female reef manta rays are strongly philopatric and do not migrate between these two island groups. Combined with restricted male-mediated migration, equivalent to a single male moving between islands every 2.2 generations (~ 64 years), we provide evidence these populations are significantly demographically isolated. Estimates of contemporary effective population size (Ne) are 104 (95% CI: 99–110) in Hawaiʻi Island and 129 (95% CI: 122–136) in Maui Nui.
Concordant with evidence from photo identification and tagging studies, these genetic results indicate reef manta rays in Hawaiʻi have small, genetically-isolated resident island populations. We hypothesize that due to the Island Mass Effect, large islands provide sufficient resources to support resident populations, thereby making crossing deep channels separating island groups unnecessary. Small effective population size, low genetic diversity, and k-selected life history traits make these isolated populations vulnerable to region-specific anthropogenic threats, which include entanglement, boat strikes, and habitat degradation. The long-term persistence of reef manta rays in the Hawaiian Islands will require island-specific management strategies.
Reef manta rays (Mobula alfredi) are planktivorous elasmobranchs that inhabit tropical and subtropical oceans between north and south latitudes of about 35 degrees [1,2,3,4]. Reef manta rays are known to have a strong affinity to specific coastal reef habitats [4,5,6,7,8,9,10,11,12,13,14,15,16,17], and spend most of their time at depths less than 50 m  feeding and visiting cleaning stations, while making occasional visits to nearshore pelagic habitats for foraging [6, 18]. Several discrete populations worldwide are reported to be in decline in large part due to anthropogenic threats including fishing and loss of coral reef habitat [3, 9, 19,20,21,22]. In addition, their conservative life history traits of slow growth, late maturity, and low fecundity [3, 23,24,25,26] can hinder recovery [21, 27] and contribute to their Vulnerable to Extinction status on the IUCN Red List of Threatened Species . Island populations may be particularly vulnerable as remote archipelagos are more likely to be demographically and genetically isolated. However, in most island regions we have limited knowledge of genetic connectivity, effective population size, and patterns of migration, which will be critical to effective conservation and management of manta ray populations.
In the past decade, genetic studies have revealed phylogenetic and geographic partitioning both among and within species of manta rays. Phylogenetic studies have provided robust support for the species level discrimination of M. alfredi from oceanic manta rays (M. birostris) and their relationship in the family Mobulidae [29,30,31,32]. For reef manta rays, strong genetic differentiation across ocean basins suggests that large areas of open ocean are effective barriers to long distance dispersal [31, 33,34,35]. Until recently, the few studies that have examined population genetic structure on smaller regional scales detected no reductions in gene flow within archipelagos [33, 36] or along continuous continental coastlines spanning up to ~ 400 km [33, 34]. However, Lassauce et al.  provided the first evidence of fine-scale genetic differentiation between reef manta ray aggregation sites on two islands of New Caledonia. While this demonstrates the potential for restrictions in gene flow between island groups, the scale at which gene flow occurs among reef manta ray aggregation sites remains unclear for most populations.
In Hawaiʻi, reef manta rays occur throughout the archipelago , but are frequently observed and well-studied at one aggregation site on Maui Nui and two on Hawaiʻi Island (Fig. 1). Regular photo-identification cataloging of the manta ray population along the western “Kona” coast of Hawaiʻi Island (hereafter referred to as “Hawaiʻi Island”) by the Manta Pacific Research Foundation  has logged 318 unique individuals with photos dating back to 1979. The Maui Nui reef manta ray, which includes individuals from the islands of Maui, Molokaʻi, Lanaʻi and Kahoʻolawe, has been well studied since 2005 with some photo-identifications dating back to 1990. The Maui Nui manta ray photo-identification catalog consists of 600 unique individuals to date , which have demonstrated regular movements between the 4-island Maui Nui complex  and are therefore presumed to be a single population. No evidence exists from either tagged individuals [9, 14] or from photo-identification matches  that reef manta rays cross the 48 km ‘Alenuihāhā Channel that separates Maui Nui and Hawaiʻi Island, however no formal genetic evaluation has been conducted.
Here, we test the hypothesis suggested by Deakos et al.  that each island group supports a demographically independent population of reef manta ray. To test the island-resident hypothesis, nuclear genome-wide single nucleotide polymorphisms (SNPs) and whole mitochondrial genomes were genotyped from 38 individuals of both populations (Hawaiʻi Island = 18, Maui Nui = 20). Using this population genomics approach, we aim to determine the magnitude and direction of gene flow among populations and explore the role of sex-biased dispersal and estimates of effective population size (Ne) to gauge population size and resilience. Overall, this study aims to provide insight into genetic stock structure and demographic parameters that can be used to inform management of reef manta rays in Hawaiʻi.
Mitogenome assembly and diversity
Mitogenomes were successfully assembled and haplotyped in 34 individuals, with an average sequencing coverage of 82.6x (min = 16, max = 366) (Fig. S2) and high-quality base calls across 86.2% of the mitogenome (min = 62.4%, max = 99.4%)(see Additional File 4 for summaries). Across individuals, an average of 10,969 reads (range: 2,802 − 45,143) were mapped to the reference mitogenome from M03 (OP562409.1; described in ). The 34 individuals were haplotyped at nine variant sites (≥ 4x coverage) with an average of 88% of variable allele calls (mean calls per individual was 7.8 of 9 variable sites). Mean read depth (across individuals) per variant site was at least 33x (grand mean = 95.7x). All nine variable sites were biallelic and segregating (all transitions, no transversions) with three synonymous changes and six replacement sites. No insertions or deletions were included in the final variant dataset; however, the control region is rich in AT-repeats that were difficult to align and could contain INDELs that are not represented in this dataset.
Mitochondrial molecular diversity indices are summarized in Table 1. Thirteen haplotypes were observed among 34 reef manta rays from the two island groups (Fig. 2). Overall nucleotide diversity was 0.00018 and haplotype diversity was 0.877. Within each population, Hawaiʻi Island had only two variable sites and Maui Nui had nine variable sites. There were no fixed differences, two shared mutations, seven unique mutations (all unique to Maui Nui), and no unique mutations to Hawaiʻi Island. Maui Nui had roughly 3-times higher genetic diversity (π = 0.00017), number of haplotypes (H = 11) and average number of nucleotide differences (K = 3.02) compared to Hawaiʻi Island (π = 0.00005, H = 4, K = 0.98). The control region was the most diverse single region with four variable sites, and the remaining five variable sites were spread among four genes: NADH5 (2), NADH4 (1), 16S ribosomal (1), and cytochrome b (1) (Table 2). The remainder of mitochondrial genes did not exhibit any variation across individuals.
Neutrality tests calculated from the mitogenome were not significant (Table 1), except for Fu’s FS in Maui Nui (FS = -4.44, P = 0.009), which suggests a recent population expansion among the Maui Nui population. This result combined with the pattern that all unique variable sites were restricted to the Maui Nui population provides evidence that Hawaiʻi Island is an ancestral population, which has since expanded northward to Maui Nui. Among island populations, we found the mean number of nucleotide differences (Kxy = 4.44), mean number of nucleotide substitutions per site (Dxy = 0.00024), mean number of net substitutions per site (Da = 0.00013), and nucleotide divergence (K = 0.00024).
Differentiation across the mitogenome
The median-joining mitogenome haplotype network revealed thirteen haplotypes among two haplogroups coinciding with island populations (Fig. 2). Based on mtDNA haplotype frequencies, we found significant genetic structure among the two island groups (AMOVA ΦST = 0.488, P < 0.0001, Table 3). Genetic differentiation in the five mitochondrial regions exhibiting variation (i.e., with variable sites) among haplotypes was consistently high and significant (ΦST = 0.544–0.731, P < 0.001) for 4 of 5 genes including 16S, NADH4, cytochrome b, and the control region (Table 2). Only NADH5 showed less structure and was not significant (ΦST = 0.169, P = 0.054). This degree of differentiation suggests the existence of strong barriers to gene flow between Maui Nui and Hawaiʻi Island populations. Furthermore, the presence of a Hawaiʻi Island haplotype from a single individual in Maui Nui (Fig. 2) supports the direction of exportation of diversity from Hawaiʻi Island to Maui Nui.
Coalescence-based migration estimates determined by MIGRATE-N showed an overall pattern of net migration from Hawaiʻi Island to Maui Nui (Fig. 3). The mitogenome showed an effective number of migrants per generation (NeM) of 0.023 (95% CI: 0.011–0.043) and 0.021 (95% CI: 0.008–0.047) moving northwesterly from Hawaiʻi Island to Maui Nui and southeasterly from Maui Nui to Hawaiʻi Island, respectively (Table 4). Migration estimates overall gene flow (mean NeM = 0.022; 95% CI: 0.01–0.045) to be equivalent to 1 female migrant moving between these island groups every 45 generations, or approximately 1305 years, based on the estimated 29-year generation time [10, 25, 28]. For every one migrant from Hawaiʻi Island to Maui Nui, the relative migration network analysis estimated that there were 0.20 migrants from Maui Nui to Hawaiʻi Island (Fig. 3), further supporting net migration from Hawaiʻi Island northwesterly to Maui Nui.
Nuclear genome scans
A total of 2048 filtered and informative nuclear genome-wide SNPs were successfully genotyped in 38 individuals. Employing the OutFLANK approach across island groups, we detected 10 SNPs putatively under divergent selection (Fig. S3) and 2038 neutral nuclear SNPs. For all subsequent analyses using nuclear loci, we analyzed the neutral loci (2038 SNPs), and outlier loci (10 SNPs) separately. While we also conducted all analyses using these two datasets combined (i.e., all 2048 nuclear loci) we do not present those results in the main text as they were very similar to neutral loci datasets and presented no change in inference. Nuclear molecular diversity indices are summarized in Table 5. For neutral nuclear loci, the mean number of alleles per locus, effective number of alleles, and the observed heterozygosity were all higher in Maui Nui (Na = 2.022, Neff = 1.594; HO = 0.508) compared to Hawaiʻi Island (Na = 2.018, Neff = 1.555; HO = 0.471). The same pattern was true and more pronounced for the outlier loci with Maui Nui (Na = 1.800, Neff = 1.508; HO = 0.446) presenting higher diversity indices than Hawaiʻi Island (Na = 1.700, Neff = 1.256; HO = 0.221). Inbreeding coefficients revealed that the influence of inbreeding is negligible across all populations for both neutral and oultier loci datasets (Table 5).
Population structure between Maui Nui and Hawaiʻi Island was detected in both neutral and outlier nuclear loci (Table 6; Fig. 4). Genetic differentiation of the neutral nuclear loci (FST = 0.003, P = 0.033) was significant but notably weaker relative to the mitochondrial genome (ΦST = 0.488). The STRUCTURE analysis based on outlier loci recovered segregating populations between island groups (K = 2, Table S5), which coincided well with the DAPC analyses (Fig. 4). The STRUCTURE analysis based on neutral loci also resolved two populations (K = 2, Fig. S4, Table S5), but the overall pattern showed a high degree of admixture and the clustering signal between islands was much lower than was observed in the oulier loci (Fig. 4). The DAPC analyses detected a stronger clustering pattern among island groups compared with STRUCTURE, but overall results coincide well from both tests. The AMOVA and STRUCTURE output using all 2048 nuclear loci (combining neutral and outlier) produced similar patterns to the neutral loci (Table S6; STRUCTURE results not included here). In contrast to the relatively low genetic differentiation at neutral loci (FST = 0.003), differentiation at outlier loci was more than 60x higher (FST = 0.186, P < 0.001), and showed a clear pattern of clustering among islands, yet with still some degree of admixture (Fig. 4).
We attempted to annotate RAD contigs containing outlier loci to explore potential ecologically relevant functions. Of 10 contigs containing outlier SNPs identified with OutFLANK, three blasted without hits (i.e., no sequence counterpart in GenBank), three had blast hits to known sequences but were unmapped, and the remaining four were mapped and annotated with gene ontology terms (Table S2). Sequence similarity scores to annotated genes were relatively low (< 75%), and there were no obvious patterns of enrichment or ecologically relevant adaptive functions of the four genes associated with outlier loci.
Coalescence-based migration estimates determined by MIGRATE-N varied among datasets but showed an overall pattern of net migration from Hawaiʻi Island to Maui Nui (Table 4; Fig. 3). Migration estimates from neutral loci showed relatively low migration rates, that are only slightly higher in the northwesterly direction NeM = 0.45 (95% CI: 0.06–0.88) compared to NeM = 0.44 (95% CI: 0.18–0.87) in the southeasterly direction. The magnitude of migration varied across nuclear datasets, however the overall pattern of net migration from Hawaiʻi Island to Maui Nui was concordant across all datasets (Table 4). Based on nuclear loci, for every one migrant from Hawaiʻi Island to Maui Nui, the relative migration network analysis estimated that there were 0.87 migrants from Maui Nui to Hawaiʻi Island (Fig. 3). Migration analysis on neutral nuclear loci estimates overall gene flow (mean NeM = 0.45; 95% CI: 0.12–0.88) equivalent to 1 effective migrant moving between these island groups every 2.2 generations, or approximately every 64 years, based on the estimated 29-year generation time. Migration estimates were 20 times greater for nuclear loci than the mitogenome. The patterns of net migration in both the nuclear and mitogenome, coupled with the results of the haplotype network, which revealed the presence of Hawaiʻi Island haplotypes in Maui Nui (but not the inverse), supports a Hawaiʻi Island to Maui Nui net direction of migration.
Estimates of contemporary genetic effective population size (Ne) based on neutral nuclear loci resulted in Ne(raw) estimates of 104 (95% CI: 99–110) for Hawaiʻi Island and Ne(raw) = 129 (95% CI: 122–136) for Maui Nui populations. Ne(adj) estimates were found to be slightly less conservative with Hawaiʻi Island and Maui Nui, estimated to be 122 (95% CI: 110–127) and 155 (95% CI: 141–157), respectively.
We investigated the fine-scale genetic structure of reef manta rays (Mobula alfredi) in the Hawaiian Islands using a combination of genome-wide nuclear loci and whole mitogenome sequences. Our assessment of genetic connectivity revealed that patterns of gene flow are associated with the geographic separation of island groups, with significant differentiation observed in the mitochondrial genome (ΦST = 0.488, P < 0.0001) as well as neutral (FST = 0.003, P = 0.033) and outlier loci (FST = 0.186, P = 0.001) across the nuclear genome (Figs. 2, 4 and 3). This pattern aligns with the expectations for island-associated populations of reef manta rays, which are known to have limited home ranges and high site fidelity [6, 7, 11, 13, 41]. Our results of fine-scale genetic structure among island groups in Hawai’i are consistent with that of Lassauce et al.  from New Caledonia and reinforce evidence that insular populations can be genetically isolated at small spatial scales (< 100 km). Furthermore, we describe new insights regarding patterns of sex-biased migration and genetic connectivity and discuss potential mechanisms that may be limiting the dispersal of reef manta rays between islands.
Patterns of dispersal and connectivity
Sex-biased dispersal among islands
Patterns of genetic differentiation of mitochondrial genomes and nuclear outlier loci both show a clear clustering among island populations (Figs. 2 and 4B). Outlier loci can reflect genomic regions associated with local adaptive differences [42, 43]. Strong differentiation at these regions putatively under selection indicates that localized selection could be overpowering the homogenizing force of male dispersal. These patterns of male-mediated dispersal and signatures of local selection are consistent with the hypothesis advanced by Portnoy et al.  that the combination of philopatric females and dispersing males may favor local adaptation by simultaneously allowing dispersal and the localized sorting of adaptive alleles.
Strong divergence in the maternally-inherited mitogenome (ΦST = 0.488) relative to the biparentally-inherited nuclear genome (FST = 0.003) is robust evidence that female reef manta rays are strongly reproductively philopatric. This signal of differentiation and geographic clustering of mitochondrial haplotypes (Fig. 2) indicates the existence of strong barriers to gene flow between Maui Nui and Hawaiʻi Island, with female reef manta rays remaining resident to, and reproducing predominantly around, the island where they were born. Comparatively, the low, but significant divergence in neutral loci across the nuclear genome indicates reduced gene flow mediated by weak male-biased dispersal between island groups. Coalescence-based migration analysis estimates movement between islands equivalent to 1 individual migrant every 64 years for males (nuclear loci) and 1305 years for females (mitogenome) (Table 4). The combination of restricted female- and male-mediated migration provides evidence these populations are demographically isolated. This asymmetry in spatial genetic patterns between nuclear loci and mitochondrial haplotypes among reef manta rays is consistent with female philopatry and weak male-biased dispersal [45, 46], which has been documented in batoids [47,48,49,50] and sharks [51,52,53,54]. Broadly, this molecular evidence for male-biased dispersal in reef manta rays adds further support to sex-biased dispersal as a recurrent pattern in viviparous elasmobranchs (reviewed in Phillips et al. ).
Female reproductive philopatry is common among elasmobranchs and widespread across batoids including several species of stingrays , skates [56, 57], and sawfish [47, 50, 55]. Previous studies (reviewed in Flowers et al. ) have provided compelling evidence of site fidelity and residency in M. alfredi that are consistent with philopatry [6, 8,9,10, 12, 15, 59,60,61,62]. However, this study is the first to confirm reproductive philopatry in reef manta rays using genetic evidence.
Evidence supporting the island-resident hypothesis
Here, we review evidence evaluating hypotheses for the drivers of population breaks between island groups. We provide support for the hypothesis suggested by Deakos et al. , that large islands provide sufficient coastal resources to support resident populations of reef manta rays, thereby making inter-island dispersal unnecessary.
The genetic patterns presented here are in concordance with evidence from photo-identification and tagging studies and thus confirm that Maui Nui and Hawaiʻi Island reef manta ray populations are distinct population stocks with restricted movement between them. First, no photo-identification matches have been made between catalogs on Hawaiʻi Island, which to date contains 318 unique individuals (1979–2023; ), and within Maui Nui, which contains 600 unique individuals (2005–2023; ). The number of individuals in the catalog that have perished since their last sighting is unknown. Second, active tracking and satellite tagging of 53 unique individuals (40 in Clark , 13 in Deakos et al.  and Deakos, unpublished) demonstrated that reef manta rays readily moved between the islands of Maui Nui but did not migrate to Hawaiʻi Island (and vice-versa). Therefore, photo-identification and tagging studies all provide no evidence of animal movements between Maui Nui and Hawaiʻi Island [9, 14].
Availability of resources may be driving this limited movement of reef manta rays between islands. Reef manta rays forage for zooplankton in generally nutrient-poor oligotrophic waters in tropical and subtropical oceans. To meet energetic needs of this large planktivore, reef manta rays require high densities of planktonic prey . Islands provide localized hotspots of productivity compared to surrounding pelagic waters. This fertilizing effect of islands, known as the Island Mass Effect , is driven by several mechanisms including island-induced mixing, nutrient flux from freshwater runoff, mesoscale eddies, and increased internal wave activity, all of which combine to enhance phytoplankton productivity near islands . Biological production scales with total reef area, thus larger islands with greater reef area exhibit increased productivity enhancements that translate up the food web .
At small scales, reef manta ray feeding events coincide with localized high biomass of zooplankton driven by fine-scale oceanographic processes, such as strong tidal currents interacting with island topography [5, 63], cold-water bores created by breaking internal waves , and surface slicks generated by Langmuir cells . On Hawaiʻi Island, surface slicks are ubiquitous and prevalent along the western coastline and have been found to accumulate dense concentrations of zooplankton . These features, which are driven by a variety of mechanisms including internal waves and Langmuir cells , have been correlated with reef manta ray feeding events in other regions  and could provide enhanced foraging opportunities near islands.
Collectively, these fine-scale and meso-scale oceanographic features, all enhanced by the Island Mass Effect, provide more biomass and greater predictability of planktonic prey near large islands . Reliable coastal resources likely eliminate the need to travel to other islands to forage and thus could be driving the strong population differentiation among islands for insular reef manta rays.
Within the Hawaiian Archipelago, similar patterns of population and genetic breaks between Hawaiʻi Island and Maui Nui have been observed in several wide-ranging insular species including common bottlenose dolphins [69, 70], spinner dolphins , rough-toothed dolphins [72, 73], and pantropical spotted dolphins . These studies similarly cite the increased productivity around the Hawaiian Islands due to the Island Mass Effect as the leading explanation driving patterns of high site fidelity and population differentiation among neighboring islands [69, 70, 72, 73].
We reject the alternative hypothesis that the lack of exchange between island groups results from isolation by distance. The shortest distance between Maui Nui and Hawaiʻi Island is only 49 km, and linear movements up to 91 km have been documented elsewhere in Hawaiʻi  and more than 200 km between nearby islands in other archipelagos (Table S3) [11, 12, 16, 60, 75, 76]. Along continental shelves with continuous coastlines, linear movements over 500 km are common [18, 60, 77, 78].
There is evidence that deep-channel crossings create habitat breaks that can be barriers to dispersal for reef manta rays, particularly when separating large islands. First, the deepest area transited over following an acoustically tagged reef manta ray in Hawaiʻi was 360 m in Maui Nui  and 300 m off Hawaiʻi Island . Tracked individuals generally followed the bottom contour depth, remaining in relatively shallow water with maximum recorded dive depths of 218 m  and 308 m (Deakos, unpublished). Hawaiʻi Island and the Maui Nui island complex are large, high volcanic islands separated by the ʻAlenuihāhā channel, which has a minimum crossing depth of 1900 m, and plummets to depths > 4700 m on both eastern and western sides (Fig. 1). Second, the few inter-island movements over deep-water that have been recorded globally all occur between small atolls with island areas less than 46 km2 (Table S3) [11, 12, 41, 75]. For example, in French Polynesia, Carpentier et al.  reported crossings between Bora Bora and Maupiti, a channel similar in span (50 km) and depth (> 3000 m) to the ʻAlenuihāhā Channel. However, these islands are only 0.2% of the area of Hawaiʻi Island (Table S3). Third, in contrast, regular inter-island movements up to 450 km have been recorded between relatively large islands in Indonesia that are connected by shallow shelves typically < 300 m [15, 16, 60, 62, 79]. These patterns reinforce the hypothesis that shallow channels and shelves create more continuous habitat that is crossed regularly, regardless of island size, but deep channels (> 300 m) are crossed infrequently and only when separating small islands or atolls. Altogether this suggests that the interaction of island size (as a proxy for coastline generated resources), channel depth, and philopatric behavior together play an important role in determining movement patterns and ultimately gene flow in reef manta ray insular populations. It remains unclear why deep-water inhibits movement, but could be due to increased exposure to predators, less foraging opportunities, and reduced ability to navigate via bathymetry.
Oceanic and archipelagic patterns of genetic connectivity
At the ocean basin scale, evidence of genetic structure among populations of reef manta rays demonstrates little potential for long-distance dispersal and migration from distant populations [33, 80]. Ocean basins are common barriers to dispersal in other elasmobranchs , particularly reef-associated sharks [54, 82, 83]. With that said, an opportunistic sighting of a pregnant reef manta ray at Cocos Island in the Eastern Tropical Pacific , nearly 6000 km from the nearest aggregation site in the Marquesas (and ~ 7500 km to Hawai’i), reminds us of their potential for oceanic dispersal, even if extremely rare.
Perhaps not surprisingly, the pattern of regional genetic structure observed in reef manta rays is in direct contrast to genetic patterns in the oceanic manta ray (Mobula birostris), which shows relative panmixia across their circumtropical distribution  albeit with indications of structure in the Eastern Tropical Pacific . The contrast in reef manta rays with high site fidelity versus oceanic manta rays with wide-ranging behavior, can explain these observed differences in genetic connectivity.
Our results of genetic structure among neighboring islands within the Hawaiian Islands provides evidence of finer-scale restrictions in gene flow than observed in surveys on continental coastlines  and within the Maldives archipelago [33, 36]. In Mozambique, Venables et al.  report high genetic connectivity along ~ 400 km of continuous coastline. Similarly, Hosegood  detected no genetic sub-structuring among islands in the Maldives spanning up to 350 km across the archipelago. In both regions, animals move more regularly between aggregation sites [59, 80, 86], which likely explains higher gene flow than that observed among island groups in Hawaiʻi (current study) and New Caledonia  where inter-island movement is restricted. Our results are comparable to the fine-scale structure observed between Grande Terre (New Caledonia) and Ouvea, which are approximately ~ 120 km apart across a deep-channel (> 2000 m) . Patterns from genetics and photo-identification in both regions [9, 35] indicate high site fidelity and few connections between aggregation sites, suggesting movement and gene flow across islands are restricted.
Patterns and implications of small effective population size estimates
Our estimates of contemporary effective population size (Ne) are 104 for Hawaiʻi Island and 129 for Maui Nui. When sampled from a mixed-age group and overlapping generations, as in this study, Ne is an approximate estimate of the harmonic mean of the number of breeders (Nb) in the population over several generations [87, 88]. There are not yet robust abundance estimates of reef manta ray populations in the Hawaiian Islands (see [9, 89]), however minimum population size can be approximated using the photo-identification catalog sizes of 318 and 600 unique individuals for Hawaiʻi Island and Maui Nui, respectively. These minimum population sizes do not consider individuals that have died. A decade or more can sometimes pass between sightings of certain individuals making it difficult to determine when an individual is no longer part of the population. While no constant relationships exist between Ne and the minimum estimates for census size (Nc) across taxa [90, 91], our results are intermediate to the only two reef manta ray Ne estimates published to date [34, 92]. The Yaeyaema, Japan population has an Ne estimated at 89 and a catalog size of 305 unique individuals . The population along the southern Mozambique coast has an Ne estimated at 375  and a catalog size of 1209 [80, 93]. Using catalog size as minimum estimates for census size (Nc), the Ne/Nc ratios in these four populations are similarly between ~ 0.2–0.3 (Table S4), suggesting that this ratio is relatively consistent among reef manta ray populations and could be useful for estimating one metric in the absence of the other.
The relatively high Ne/Nc ratios observed in reef manta rays are consistent with species with Type I survivorship curves, which well characterize manta rays and other viviparous elasmobranchs . Several studies have demonstrated the variation in Ne/Nc among taxa is driven primarily by age-at-maturity, adult lifespan, and variation in reproductive success [94,95,96]. Reef manta rays have a delayed age-at-maturity (8–17 years), long lifespans up to 45 years [10, 25], and low and variable reproductive output of a single pup every 1 to 7 years for mature females [23,24,25, 79, 92].
We acknowledge that the relatively small sample sizes used in this study may reduce precision in our estimates of contemporary Ne. Precision and accuracy of Ne estimates derived from low numbers of samples (and loci) can decline resulting in infinite parameters or upper confidence limits . However, our study utilized high numbers of loci (> 2000 SNPs), which improve precision in Ne for small population sizes . Furthermore, the Ne confidence limits we report are finite and relatively narrow for both Hawai’i Island (104; 95% CI: 99–110) and Maui Nui (129; 95% CI: 122–136) populations, suggesting these data have sufficient power and precision for estimating Ne. Additionally, our relative number of samples to Ne estimates (13–15% of Ne) are above the 10% of Ne guideline proposed by Palstra & Ruzzante  and similar to estimates of other elasmobranchs [100, 101]. Efforts to increase sample size and geographic coverage are ongoing and are expected to improve precision in future genetic surveys of reef manta rays in the Hawaiian Islands.
Most regions support relatively small populations of reef manta rays, typically less than 1,000 individuals. Long-term photo-identification studies have produced minimum estimates (i.e., catalog sizes) as low as 54 animals in Yap, Micronesia  and up to 4,411 individuals in the Maldives . Larger population sizes (over 600) tend to be associated with continental shelves such as the East Coast of Africa , Australia , as well as archipelagos consisting of many islands connected by shallow water like the Maldives [59, 86], and Indonesia [16, 76]. Smaller populations tend to be associated with remote archipelagos with islands separated by deep-water including Hawai’i [9, 14], French Polynesia , New Caledonia , and Seychelles [11, 103] and generally have small home ranges. These smaller home ranges can likely be attributed to having access to sufficient coastal resources, cleaning stations, mates, and protection from predation [5, 104].
Patterns and consequences of low genetic diversity
We conclude that the low genetic diversity observed in reef manta rays is due to low mutation rates combined with inherently small, localized populations. The low levels of genome-wide diversity observed in Hawaiian reef manta rays (Tables 1 and 5) are generally consistent across elasmobranchs (e.g., [105, 106]). For example, a comparison of whole mitogenomic diversity of the endangered speartooth shark (Glyphis glyphis) revealed strikingly similar patterns of genetic variation (i.e., π = 0.00019, h = 0.76; ) to those of M. alfredi (π = 0.00018, h = 0.88; Table 1). In G. glyphis, the low genetic diversity was attributed to low mutation rates and a low effective population size, which follow patterns in other elasmobranchs [106, 108]. Although empirical estimates of mutation rates do not yet exist for reef manta rays, the rates of mitochondrial mutations in elasmobranchs are slow relative to other taxa . With concordant patterns within a variety of evolutionary distinct elasmobranchs, it is reasonable to conclude low mutation rates are naturally occurring phenomena in mobulid rays and, along with small effective population size, contributes to the low diversity observed in reef manta rays.
The concern regarding the consequences for populations with low genetic diversity is compounded when population sizes are small, as in reef manta rays. When low genetic diversity is observed in a population, it is often interpreted as an indication of inbreeding depression  and is thought to compromise reproductive fitness and capacity for population growth . However, recent research has challenged the assumption that high levels of genetic diversity are necessary for increased fitness and survival and reframed the importance of genetic diversity when considering variation that is relevant to future environmental and climatic changes . Based on patterns in other elasmobranchs and other populations of reef manta rays, low population size may be the natural state of this species and low genetic diversity may not be a suitable barometer for evaluating population health and extinction risk. There is evidence of a population expansion for the Maui Nui population based on a single metric (Table 1; Fu’s Fs = -4.44. P < 0.01), however, other neutrality tests (i.e., Tajima’s D or Fu and Li’s D) did not show support of a bottleneck or recent expansion, suggesting these Hawaiian populations have not fluctuated dramatically in size. Thus, small population sizes and low genetic diversity may be the natural biological state for this species at least in this region.
Management implications of island-resident populations
Combined with evidence from photo-identification and tagging studies, these genetic results indicate reef manta rays in the Hawaiian Islands have small resident island populations that are significantly genetically isolated and should be managed as discrete stocks. The lack of female migration among islands means extremely little potential for replacement females to enter and establish from other islands. The relatively low levels of male-mediated migration still indicate that replacement males cannot be counted on to replenish island populations should they decline. On the ocean basin scale, evidence of genetic differentiation between Hawai’i and the South Pacific and Indian Oceans  provides further evidence that Hawaiian reef manta rays are demographically isolated from other regions. The high degree of residency, low genetic connectivity, and geographic isolation of the Hawaiian archipelago all suggest there is little potential for replenishment from distant aggregation sites. Together with small population size, restricted gene flow, low genetic diversity, and conservative life history traits, this leaves reef manta rays at extremely high risk to human-induced perturbations. Reef manta rays have low intrinsic growth rates due to their delayed age-at-maturity (8–17 years for females) and low fecundity [9, 28, 89]. These extremely conservative life history traits are expected to severely restrict the potential for recovery from any potential population reductions in the future.
Globally, several reef manta ray populations have been reported to be in decline  and has resulted in their listing on Appendix II of the Convention for the International Trade in Endangered Species of Wild Fauna and Flora (CITES) in 2013. Much of the decline has been attributed to direct manta ray hunts to provide for the global demand of gill plates [4, 8, 21, 22]. In Hawaiʻi, manta rays are neither fished intentionally nor known to be caught as bycatch in any fishery. In 2009 a state bill was passed that protects manta rays from being killed or captured. Despite these much-needed protections, sighting rates at a reliable cleaning station on the island of Maui have declined by over 90% in the past decade (Deakos, unpublished). Whether this is due to a reduction in population size or manta rays relocating elsewhere due to degradation of reef habitat remains unknown and requires further investigation. In addition, more than 10% of the Maui Nui population has evidence of entanglement in fishing line, primarily from shore-casting fishing gear used to target giant trevally (Caranx ignobilis) . Injury to the cephalic fins from fishing line could negatively impact feeding efficiency and mate attraction and thus may have long-term consequences even if animals survive.
Monitoring long-term population trends and breeding stock will be important for evaluating the population level impacts of local anthropogenic threats, which include fishing line entanglement [9, 20], pollutants and contaminants , plastic ingestion , boat strikes [20, 41], pressures from commercial manta ray dive tours [14, 22], habitat degradation as a result of coastal development [86, 114], as well as projected declines in zooplankton  and climate change.
Further research to assist with effective management strategies should include the following. First, identifying essential habitat areas that are used for cleaning, feeding, mating, and pupping  could be achieved by expanding telemetry studies to develop core use density maps . Second, with critical habitats defined, a focus should be made on quantifying and reducing regional anthropogenic threats, especially habitat degradation around nursery and aggregation sites, both important for reproduction. In addition, an assessment of the impact of commercial tourism activities should be conducted, including evaluating mitigations such as establishing codes of conduct , setting carrying capacity for boats/divers and implementing rest-periods in heavy use areas. Third, robust mark-recapture population models should be built by expanding photo-ID and acoustic tagging efforts to improve survival and abundance estimates (e.g., ) or adopting a modification of the standard POPAN model that incorporates per capita recruitment and transience parameters to estimate annual population sizes . Fourth, genetic sampling should be expanded across the archipelago and years to evaluate gene flow, monitor changes in contemporary Ne , and assess population trends. Fifth, building population projection models that incorporate target prey availability (informed by empirical study of prey density thresholds) and examine extirpation risk can be used to predict population trends under different management and climate change scenarios. Collectively, these efforts coupled with engagement with the public and stakeholders will be critical to ensuring the long-term persistence of healthy populations of island-resident reef manta rays in the Hawaiian Islands.
Long-term stability of aggregation sites is beneficial for the stability of social structure, particularly mating behavior. In addition to benefits of cleaning, coral reef cleaning stations also serve as hubs for reproductive activity [9, 23, 25]. Degradation of nearshore coral reef habitat serving as mating aggregation sites, like in Maui , is expected to negatively influence reproductive success and could have long-term demographic consequences. The loss or degradation of embayments that can serve as pupping/nursery grounds [62, 113] can shift patterns of natural selection and negatively influence demography . Degradation of coastal coral reefs and embayments from coastal development and climate change represents a serious threat to reef manta rays in Hawaiʻi and elsewhere.
We provide evidence of strong female reproductive philopatry and weak male-mediated dispersal that indicate genetic isolation of small resident island populations of reef manta rays in Hawaiʻi. Despite the proximity of these island-associating populations, they represent two genetically distinct stocks with varying geographic characteristics. The Hawaiʻi Island population demonstrates a much tighter home range but has immediate access to deep-waters whereas the Maui Nui population shares a relatively shallow bathymetry between the 4-island Maui Nui complex. The threats affecting Maui Nui’s population may be more dependent on entanglement in coastal fishing line and habitat degradation, whereas the Hawaiʻi Island population face a different set of challenges from direct interaction with boats and divers. These distinct differences between neighboring island populations, combined with their reproductive isolation and vulnerable life history characteristics, highlight the importance of local, island-specific management strategies.
Most data were collected opportunistically while freediving or with open-circuit SCUBA, either from a boat or from a shoreline entry at known manta ray aggregation areas (Table 7; Fig. 1). Biopsy samples were collected from March 2012 to October 2015. For each manta ray encountered, attempts were made to get a photo-identification of the ventral side, a gender and age-class determination (juvenile or adult) based on clasper development in males [9, 23, 119] or mating scars and visible pregnancy in females , and body size measurements using paired-laser photogrammetry as described in Deakos . Biopsy samples were obtained using a modified Hawaiian sling containing a stainless-steel cylindrical biopsy tip (13 mm in length x 5 mm in diameter) that extracts a sample of skin and muscle. Biopsies were taken from the caudal end of the manta ray’s disc, to avoid sampling close to the main body trunk. Samples were preserved in either 20% salt-saturated dimethyl sulphoxide (DMSO) or 95% ethanol and stored at -20 °C. Biopsy tips were washed and sterilized in a 10% bleach solution for 5 min, rinsed with fresh water, soaked for 5 min in 95% ethanol, air dried, and placed individually into small bags for reuse.
Surveys on Maui took place during the daytime, and primarily focused on a known manta ray cleaning station located at the south end of the Olowalu Reef (20.7913°N, -156.5880°W) and within a kilometer of the West Maui shoreline  (Table 7; Fig. 1). Hawaiʻi Island daytime surveys targeted manta rays opportunistically visiting cleaning stations or surface feeding in current lines within a kilometer off the West Hawaiʻi Island shoreline (Table 7; Fig. 1). Nighttime surveys were conducted at two popular commercial manta ray snorkel and dive locations at Keahou Bay (19.5585°N, -155.9669°W) and Makako Bay (19.7347°N, -156.0567°W) where manta rays feed on zooplankton attracted to underwater lights .
Whole genomic DNA was extracted from each sample using an Omega E-Z 96 Tissue DNA Kit (Omega), following the manufacturer’s protocol. DNA extracts were quantified using the AccuBlue High-Sensitivity dsDNA kit (Biotium, USA) on a Spectramax M3 fluorescent plate reader (Molecular Devices, USA) and visualized using gel electrophoresis. Samples were normalized (to 40ul) and 1-3ug of DNA per sample were digested overnight with DpnII restriction enzyme (NEB). We used the ezRAD approach  to construct restriction-associated digest (RAD) reduced representation libraries with DpnII (GATC cut site) following the ezRAD protocol  modified to a with-bead protocol (Additional file 3) using the Kapa HT TruSeq library preparation Kit (Roche Sequencing). In summary, fragmented DNA is end-repaired, 3’ ends adenylated, and ligated with Illumina TruSeq HT dual-indexed adapter sequences (IDT). DNA fragments from 300 to 425 bp (target insert sizes 200-300 bp) were isolated using a PippenPrep automated electrophoresis system (Sage Science). Adapter-ligated, size-selected fragments were then amplified using PCR (see Additional File 3 for conditions). Following each step, samples were cleaned using AMPureXP paramagnetic beads (Beckman-Coulter), which were left in the cleaned samples and reactivated by adding 2.5 M NaCl 20% PEG (Polyethylene glycol) to the solution at various steps . DNA concentration was quantified following each step using Accublue Quantitation (Biotium). A Bioanalyzer (Applied Biosystems) was used to check size-distribution and quality of final amplified libraries. Individually barcoded libraries were normalized in equimolar concentrations (150ng) and combined in equal proportions into a single library per island. The two pooled libraries were cleaned with a final 1:1 bead cleanup and sequenced on one lane of an Illumina HiSeq 3000 (PE 150) at the UCLA Technology Center for Genomics and Bioinformatics. Raw sequenced reads were demultiplexed by index and barcode by the sequencing facility, and only samples with matched index pairs were retained, thereby eliminating index mis-assignment. Sequencing produced 268.4 million raw 150 bp sequences across 46 libraries (Hawaiʻi Island = 21; Maui Nui = 25) of Mobula alfredi. We recovered 9,872 to 16,825,862 reads (5.8 ± 4.2 million; mean ± sd) per individual library (Table S1).
Mitogenome assembly and haplotyping
Due to the high prevalence of the GATC cut site in mitochondrial genomes combined with stochastic fragmentation, whole mitogenomes can be assembled from ezRAD libraries [40, 121, 124, 125]. We assembled the complete M. alfredi mitogenome (including the control region) from individual M03 from Maui (OP562409.1, described in ), which we used as the reference for aligning and haplotyping. Raw reads were subject to QA/QC, adapter trimming, and mate-pairing validation as described below for nuclear SNPs. Cleaned paired reads for 38 individual M. alfredi (Hawaiʻi Island = 18, Maui Nui = 20) were aligned to this reference mitogenome using BWA v.0.7.17 , Samtools v.1.6 , and Bamtools v.2.5.2 . Freebayes v.1.2.0  was used to call variant sites on all haplotype alleles (including reporting on all monomorphic sites) with a minimum of 4x coverage to call. These settings included: ploidy = 1, polymorphic site must be variable in at least 2 reads, minimum quality score 20, minimum mapping quality 15, minimum coverage of 4x. After filtering for low depth individuals (missing calls > 60% of sites), 34 individuals remained in the dataset. On average, those individuals had 82.6x coverage across 77.1% of the mitogenome, mean read depth from 16x to 366x, and high-quality base calls (≥ 4x) between 62.4 and 99.4% of the mitogenome (Additional file 4). Before filtering, we identified 62 potential variable sites (at ≥ 4x) across the mitogenomes: 11 of which were shared among two or more individuals, while the majority (51) were present in a single individual (singletons). Many of the singleton sites had low read coverage and/or multiple alleles present per individual (i.e., appear as heterozygotes), which strongly indicate sequencing or mapping errors when haploids. We evaluated the relationship between the number of variable sites and the minimum read depth to determine a variant site filtering threshold (Fig. S1). We observed that the number of shared variable sites was relatively constant across read depths while the number of singleton sites decreased with increasing coverage and plateaus at a minimum average read depth of ~ 10x. Therefore, in selecting mitogenome variant sites for analysis we applied a minimum average coverage threshold of 10x, which resulted in 12 potential variant sites. We applied one additional filter using mapping discrepancies between the reference and the alternate alleles, which filtered three singletons with high coverage of both reference and alternate alleles in a single individual (i.e., appeared as heterozygotes) and are likely due to sequencing/mapping errors in haploids. The application of these two filters resulted in a final set of nine variant sites across the mitogenome, including eight shared mutations (parsimony informative) and one singleton site. We then used a less stringent coverage threshold of 4x to call alleles at those nine variant sites. We then extracted mitogenome haplotypes into fasta files for 34 individuals using VCFTools v.0.1.12a  and masked sites with less than 4x coverage with N’s.
Mitogenome data analysis
We calculated molecular diversity indices using DnaSP v.6.12.01 . We imputed missing data in genodive v.2.0b27 , using population specific allele frequencies. There was a total of 41 missing allele calls (13.4%) out of 306 alleles (34 individuals x 9 variable sites). Imputed missing data were not used for any other tests or inference outside of estimating diversity indices in DnaSP. Neutrality tests, Tajima’s D  and Fu’s Fs  were performed using arlequin v.220.127.116.11.  to infer demographic history. Analyses of molecular variance (AMOVA) were conducted in arlequin, to test for genetic structure among individuals from different islands and estimate the degree of genetic differentiation (ΦST) among islands, using a Tamura-Nei model of nucleotide evolution (selected by jModelTest 2 ). Using 10,000 permutations, we calculated pairwise ΦST for the whole mitogenome and separately for five individual genes with variable sites including four coding regions (16S, NADH4, NADH5, cytochrome b) and the control region. A haplotype network was constructed for the mitogenome with network v.18.104.22.168  using a median-joining algorithm  and default settings.
Nuclear SNP discovery and genotyping
Our nuclear SNP discovery and genotyping workflow included five major steps: (1) quality filtering (adapter cleaning and mate-pair validation), (2) de-novo pseudo-reference generation, (3) alignment to pseudo-reference RAD contigs, (4) variant calling, and (5) SNP quality control filtering. For quality filtering, raw reads were assessed using FastQC  and reads with sequencing adapters present were filtered out using Cutadapt v1.11  with maximum error rate (10%), minimum overlap 12 bases, and search algorithm (-b; found anywhere in the sequence). Orphaned reads were removed leaving only mate-paired reads. A total of 134,119,542 mate-paired reads remained after adapter cleaning and quality filtering (Table S1).
De-novo pseudo-reference generation
A pseudo-reference of ezRAD loci for Mobula alfredi was generated de novo using Seanome . In short, we input sequences from 4 individual reef manta rays (two from each island group) with the highest sequence read count following cleaning (each with 4–6 million reads). Properly mated paired reads were merged using Pear v.0.9.10  using default settings and minimum assembly length of 100 bases. Then commonly shared regions (CSRs) were found using minimum length of 100 bases for a shared region (contig), and minimum similarity of 95% to be clustered together. A total of 22,098,270 merged reads from four individuals were assembled into 872,500 commonly shared regions between 100 and 650 bases (mean = 202 bases; 99% between 100 and 300 bases) with an average of 25x coverage. We then used Usearch v8.1  to remove duplicate and reverse complement contigs and cluster unique sequences with 95% similarity (centroids) into 734,031 clusters. We then used mothur v.1.4.3  to exclude contigs: with homopolymers ≥ 8 bases, ambiguities (i.e., no Ns), lengths < 200 bases, and that do not start and end with restriction enzyme cut site (GATC). These filtering steps produced a final output of 359,756 reference contigs (RAD loci) with an average length of 220 bases (range 200–529) and a total length of 79,380,335 bases. We then used local BLASTn to query pseudo-reference contigs (359,756) to the M. alfredi mitogenome (GenBank Accession: OP562409). Five contigs matched well to the mitogenome reference (98–100% identity), all with Scores > 400, very low E-values (1e-150) and overlapping lengths between 271 and 337 bases. We checked if any of these five contigs matched with the final post-filtering SNP set contigs (2048 SNPs on 459 contigs) and none were present. Therefore, this confirms the final SNP set contains only nuclear loci, with no hits to the mitogenome reference > 100 bases and an e-score > 1e-20.
Alignment, variant calling, and SNP filtering
The dDocent pipeline  was used to map reads with BWA and call variants with FreeBayes. We indexed our reference contigs using Samtools, BWA, and Picard tools v.1.102 . Fastq sequences were adapter-cleaned, mate-paired, and all filtered to minimum length of 150 bases. Eight libraries with too few reads (< 150,000 cleaned-mated reads) were excluded from alignments and all downstream analyses (Table S1). Sequences from the 38 remaining individuals were aligned to our pseudo-reference with BWA (BWA-mem, paired sequences) with settings based on dDocent recommendations. Using Samtools we extracted only properly paired mappings and excluded unmapped reads. We performed variant calling with FreeBayes using the following settings: minimum mapping quality (5), minimum quality score (15), read-max-mismatch fraction (0.2), mismatch-base-quality threshold (10), read indel limit (5), min-alt total (10), read mismatch limit (20), and max 4 alleles.
SNPs in nuclear RAD loci were extensively filtered before analysis using a workflow (see Additional file 3 for detailed SNP filtering workflow) based on the dDocent protocol  and following recommendations from O’Leary et al.  and Portnoy et al. . In summary, the initial dataset of 49,028 SNPs output from FreeBayes was first filtered to remove all genotypes with < 5 reads per individual, quality scores < 25, and loci called in < 50% of individuals. SNPs were then filtered to meet the following criteria: called in 85% of all individuals (i.e., < 15% missing data per SNP); minor allele count > 2, minor allele frequency > 5% across all individuals, and conforming to expectations of Hardy-Weinberg equilibrium (HWE; P < 0.01 in both populations). This included the removal of SNPs with: low quality-to-depth ratios (< 0.25), discrepancies between mapping qualities and properly paired status of reference and alternate alleles, mean read depth < 20 and > 120x; and the removal of RAD loci (contigs) with excess SNPs (> 24) and those identified as possible paralogs. The final nuclear SNP set contained 2048 SNPs across 459 RAD loci (mean 4.5 SNPs per RAD locus) with 38 individuals (Table 7) containing high quality genotypes called in > 95% of all sites (i.e., < 5% missing genotypes per individual). Raw variants were filtered sequentially using VCFtools, VcfLib v.1.0.3 [129, 148], Rad Haplotyper v.1.1.9  and dDocent bash scripts from Puritz et al.  detailed steps in Additional file 3. VCF and other file format conversions were executed using PGDSpider v.22.214.171.124 . Relatedness of individuals was assessed in VCFtools, using the statistic of Yang et al. . The average relatedness among individuals was − 0.02 and no pairwise comparison had a relatedness greater than 0.15. We conclude that the final genotype set has no closely related individuals.
Outlier loci detection and annotation
We assessed population structure using three nuclear SNP datasets: all loci (2048), neutral loci (2038), and outlier loci (10), which allowed us to make inferences regarding selection. To identify neutral and outlier loci we used the R package OutFLANK [152, 153], which estimates a null distribution of FST for loci unlikely to be under strong positive selection. We ran OutFLANK implementing a left and right trim factor of 0.05, a minimum heterozygosity of < 10%, and a false discovery rate of 5% (q = 0.05). The contigs (RAD loci) containing outlier SNPs identified using OutFLANK were used as queries (e-value < 10− 6) against the nr database at the National Center for Biotechnology Information (NCBI), using Blast2GO Pro  to find homologous sequences, map and annotate Gene Ontology (GO) terms (e-value < 10− 6, annotation cut-off > 55 and a GO weight > 5).
Characterizing island populations
Genodive was used to generate genetic diversity indices for all three datasets, as well as to test for population structure. Genetic structure among sample locations was evaluated with an analysis of molecular variance (AMOVA) in arlequin. Deviations from null distributions were tested with non-parametric permutation procedures (N = 9999). Pairwise FST statistics were generated to assess genetic structure between locations. False discovery rates were controlled for and maintained at α = 0.05 among all pairwise tests [155, 156]. Genetic partitioning was assessed for all datasets using structure v.2.3.2 , a Bayesian method that estimates ancestry and categorizes individuals into discrete populations. Sampling location was provided for each individual and set as priors; the admixture model and allele frequency correlated model were implemented in structure to assess ancestry. The simulation was run for 1 million generations with the first 100,000 discarded as burn-in. Five replicates of each simulation from K = 1 to 5 genetic clusters were run. We determined the most likely number of genetic clusters (K) indicated from the Evanno method  and selecting the clusters inferred from delta K vs. K in structure harvester v.0.6.93 . Structure results were analyzed and visualized using the on-line tool CLUMPAK , which integrates the program clumpp v.1.1.2  and minimizes the variance across all iterations. We also tested population structuring for the neutral and outlier loci datasets using a discriminant analysis of principal components (DAPC) in Adegenet v4.0.2 for R . We imposed the number of clusters (K) of two and ran analyses without any priors for both neutral and outlier loci. Each DAPC was performed using one discriminant function, which is the maximum when K = 2, and the optimal number of principal components (12 for neutral loci and 1 for outliers) was chosen using the a-spline optimization procedure . Contemporary effective population size estimates for each island were calculated for each of the datasets using the molecular co-ancestry method of  as implemented in NeEstimator v.2.1 . Here, we use the point random mating Linkage Disequilibrium (LD) method and report estimates for critical values (the criterion for excluding rare alleles = 0.05). We applied a physical linkage correction factor following Venables et al.  and report both the raw and adjusted estimates of Ne, and the parametric 95% confidence levels.
Direction and magnitude of migration among islands
To examine the direction and magnitude of migration between Maui Nui and Hawaiʻi Island, migration rates were calculated using MIGRATE-N v.4.4.3 [165, 166] separately on both the neutral nuclear loci (2038 SNPs) and mitochondrial haplotypes. Several test runs were conducted to determine the appropriate prior values for the parameters θ (four times effective population size multiplied by mutation rate per site per generation, 4Neµ) and M (immigration rate divided by the mutation rate, m/µ). In the final analyses the mean prior values for θ and for M were set in both directions (i.e., between islands). After checking for data convergence, the mode and 95% percentiles of θ and M were used to calculate the effective number of migrants per generation (Me) between populations to determine the direction and magnitude of migration. A relative migration network was constructed in R using the divMigrate function in the diveRsity package [167, 168] and implemented using Jostʻs D  statistical method, based on 1000 bootstrap replicates. We took the inverse migration rate (1/Me) to calculate the number of generations needed to achieve 1 effective migrant between populations, then multiplied by the estimated generation time for M. alfredi (29 years; [10, 25, 28]) to provide an approximate number of years every 1 individual effectively migrates between populations.
The datasets supporting the conclusions of this article are available on Zenodo https://doi.org/10.5281/zenodo.7312296 or are included within the article (Additional files). Mitogenome reference was deposited in GenBank (Accession: OP562409; ). Raw sequences were deposited in NCBI’s Sequence Read Archive: https://www.ncbi.nlm.nih.gov/sra/PRJNA885488; SRA Accessions: SRR21884470-SRR21884507; BioProject https://www.ncbi.nlm.nih.gov/bioproject/885488.
Marshall AD, Compagno LJ, Bennett MB. Redescription of the genus Manta with resurrection of Manta alfredi. Zootaxa. 2009;2301:1–28.
Last P. In: Stevens J, editor. Sharks and Rays of Australia. 2nd ed. Harvard University Press; 2009.
Stewart JD, Jaine FRA, Armstrong AJ, Armstrong AO, Bennett MB, Burgess KB, et al. Research Priorities to support effective Manta and Devil Ray Conservation. Front Mar Sci. 2018;5:314.
Couturier LIE, Marshall AD, Jaine FRA, Kashiwagi T, Pierce SJ, Townsend KA, et al. Biology, ecology and conservation of the Mobulidae. J Fish Biol. 2012;80:1075–119.
Dewar H, Mous P, Domeier M, Muljadi A, Pet J, Whitty J. Movements and site fidelity of the giant manta ray, Manta birostris, in the Komodo Marine Park, Indonesia. Mar Biol. 2008;155:121–33.
Braun CD, Skomal GB, Thorrold SR, Berumen ML. Movements of the reef manta ray (Manta alfredi) in the Red Sea using satellite and acoustic telemetry. Mar Biol. 2015;162:2351–62.
Kessel ST, Elamin NA, Yurkowski DJ, Chekchak T, Walter RP, Klaus R, et al. Conservation of reef manta rays (Manta alfredi) in a UNESCO World Heritage Site: large-scale island development or sustainable tourism? PLoS ONE. 2017;12:e0185419.
Marshall AD, Dudgeon CL, Bennett MB. Size and structure of a photographically identified population of manta rays Manta alfredi in southern Mozambique. Mar Biol. 2011;158:1111–24.
Deakos M, Baker J, Bejder L. Characteristics of a manta ray Manta alfredi population off Maui, Hawaii, and implications for management. Mar Ecol Prog Ser. 2011;429:245–60.
Couturier LIE, Dudgeon CL, Pollock KH, Jaine FRA, Bennett MB, Townsend KA, et al. Population dynamics of the reef manta ray Manta alfredi in eastern Australia. Coral Reefs. 2014;33:329–42.
Peel LR, Stevens GMW, Daly R, Daly CAK, Collin SP, Nogués J, et al. Regional movements of reef Manta Rays (Mobula alfredi) in Seychelles Waters. Front Mar Sci. 2020;7:558.
Andrzejaczek S, Chapple T, Curnick D, Carlisle A, Castleton M, Jacoby D, et al. Individual variation in residency and regional movements of reef manta rays Mobula alfredi in a large marine protected area. Mar Ecol Prog Ser. 2020;639:137–53.
Harris JL, Hosegood P, Robinson E, Embling CB, Hilbourne S, Stevens GMW. Fine-scale oceanographic drivers of reef manta ray (Mobula alfredi) visitation patterns at a feeding aggregation site. Ecol Evol. 2021;11:4588–604.
Clark T. Abundance, home range, and movement patterns of manta rays (Manta alfredi, M. birostris) in Hawai‘i. PhD Dissertation ed. University of Hawaii; 2010. pp. 1–149.
Setyawan E, Sianipar AB, Erdmann MV, Fischer AM, Haddy JA, Beale CS et al. Site fidelity and movement patterns of reef manta rays (Mobula alfredi: Mobulidae) using passive acoustic telemetry in northern Raja Ampat, Indonesia. Nat Conserv Res. 2018;3.
Germanov ES, Pierce SJ, Marshall AD, Hendrawan IG, Kefi A, Bejder L, et al. Residency, movement patterns, behavior and demographics of reef manta rays in Komodo National Park. Peerj. 2022;10:e13302.
Couturier L, Newman P, Jaine F, Bennett M, Venables W, Cagua E, et al. Variation in occupancy and habitat use of Mobula alfredi at a major aggregation site. Mar Ecol Prog Ser. 2018;599:125–45.
Jaine F, Rohner C, Weeks S, Couturier L, Bennett M, Townsend K, et al. Movements and habitat use of reef manta rays off eastern Australia: offshore excursions, deep diving and eddy affinity revealed by satellite telemetry. Mar Ecol Prog Ser. 2014;510:73–86.
Heinrichs S, O’Malley M, Medd H, Hilton P. Manta ray of hope: global threat to manta and mobula rays. Manta Ray of Hope Project. 2011;(www.mantarayofhope.com).
Pate J, Marshall A. Urban manta rays: potential manta ray nursery habitat along a highly developed Florida coastline. Endanger Species Res. 2020;43:51–64.
Croll DA, Dewar H, Dulvy NK, Fernando D, Francis MP, Galván-Magaña F, et al. Vulnerabilities and fisheries impacts: the uncertain future of manta and devil rays. Aquat Conserv Mar Freshw Ecosyst. 2016;26:562–75.
O’Malley M, Townsend KA, Hilton P, Heinrichs S, Stewart JD. Characterization of the trade in manta and devil ray gill plates in China and South-east Asia through trader surveys. Aquat Conserv Mar Freshw Ecosyst. 2017;27:394–413.
Marshall AD, Bennett MB. Reproductive ecology of the reef manta ray Manta alfredi in southern Mozambique. J Fish Biol. 2010;77:169–90.
Deakos MH. The reproductive ecology of resident manta rays (Manta alfredi) off Maui, Hawaii, with an emphasis on body size. Environ Biol Fish. 2012;94:443–56.
Stevens G. Conservation and population ecology of manta rays in the Maldives. PhD Dissertation, University of York, United Kingdom; 2016.
Lawson JM, Fordham SV, O’Malley MP, Davidson LNK, Walls RHL, Heupel MR, et al. Sympathy for the devil: a conservation strategy for devil and manta rays. Peerj. 2017;5:e3027.
Dulvy NK, Pardo SA, Simpfendorfer CA, Carlson JK. Diagnosing the dangerous demography of manta rays using life history theory. PeerJ. 2014;2:e400–19.
Marshall A, Barreto R, Carlson J, Fernando D, Fordham S, Francis M et al. Mobula alfredi. The IUCN Red List of Threatened Species 2019. 2019;e.T195459A68632178.
Kashiwagi T, Marshall AD, Bennett MB, Ovenden JR. The genetic signature of recent speciation in manta rays (Manta alfredi and M. birostris). Mol Phylogenet Evol. 2012;64(1):1–29.
White WT, Corrigan S, Yang L, Henderson A, Bazinet A, Swofford D, et al. Phylogeny of the manta and devilrays (Chondrichthyes: Mobulidae), with an updated taxonomic arrangement for the family. Zool J Linn Soc. 2018;182:50–75.
Hosegood J, Humble E, Ogden R, de Bruyn M, Creer S, Stevens GMW, et al. Phylogenomics and species delimitation for effective conservation of manta and devil rays. Mol Ecol. 2020;29:4783–96.
Hinojosa-Alvarez S, Walter RP, Díaz-Jaimes P, Galván-Magaña F, Paig-Tran EM. A potential third Manta Ray species near the Yucatán Peninsula? Evidence for a recently diverged and novel genetic Manta group from the Gulf of Mexico. PeerJ. 2016;4:e2586–20.
Hosegood J. Genomic tools for conservation and management of Manta and Devil Rays (Mobula Spp). United Kingdom: Bangor University; 2020. PhD Dissertation.
Venables SK, Marshall AD, Armstrong AJ, Tomkins JL, Kennington WJ. Genome-wide SNPs detect no evidence of genetic population structure for reef manta rays (Mobula alfredi) in southern Mozambique. Heredity. 2020;126:308–19.
Lassauce H, Dudgeon C, Armstrong A, Wantiez L, Carroll E. Evidence of fine-scale genetic structure for reef manta rays Mobula alfredi in New Caledonia. Endanger Species Res. 2022;47:249–64.
Clark TB. Population structure of Manta birostris (Chondrichthyes: Mobulidae) from the Pacific and Atlantic Oceans. Masters Thesis, Texas A&M University, College Station, USA; 2002.
Randall JE. Reef and shore fishes of the Hawaiian Islands. University of Hawaii Press, Honolulu, Hawaii.; 2007.
Manta Pacific Research Foundation. http://www.mantapacific.org. Accessed 1 Oct 2022.
Maui Photo Identification Catalog. Hawaii Association for Marine Education and Research Webpage. www.hamerinhawaii.org.
Whitney JL, Coleman RR, Deakos MH. The complete mitochondrial genome of the reef Manta Ray, Mobula alfredi, from Hawaii. Mitochondrial DNA Part B. 2023;8:197–203.
Carpentier AS, Berthe C, Ender I, Jaine FRA, Mourier J, Stevens G, et al. Preliminary insights into the population characteristics and distribution of reef (Mobula alfredi) and oceanic (M. birostris) manta rays in french polynesia. Coral Reefs. 2019;38:1197–210.
Nielsen EE, Hemmer-Hansen J, Larsen PF, Bekkevold D. Population genomics of marine fishes: identifying adaptive variation in space and time. Mol Ecol. 2009;18:3128–50.
Allendorf FW, Hohenlohe PA, Luikart G. Genomics and the future of conservation genetics. Nat Publishing Group. 2010;11:697–709.
Portnoy DS, Puritz JB, Hollenbeck CM, Gelsleichter J, Chapman D, Gold JR. Selection and sex-biased dispersal in a coastal shark: the influence of philopatry on adaptive variation. Mol Ecol. 2015;24:5877–85.
Melnick DJ, Hoelzer GA. Differences in male and female macaque dispersal lead to contrasting distributions of nuclear and mitochondrial DNA variation. Int J Primatol. 1992;13:379–93.
Dudgeon CL, Blower DC, Broderick D, Giles JL, Holmes BJ, Kashiwagi T, et al. A review of the application of molecular genetics for fisheries management and conservation of sharks and rays. J Fish Biol. 2012;80:1789–843.
Feutry P, Kyne PM, Pillans RD, Chen X, Marthick J, Morgan DL, et al. Whole mitogenome sequencing refines population structure of the critically endangered sawfish Pristis pristis. Mar Ecol Prog Ser. 2015;533:237–44.
Green M, D’Anastasi B, Hobbs J, Feldheim K, McAuley R, Peverell S, et al. Mixed-marker approach suggests maternal philopatry and sex-biased behaviours of narrow sawfish Anoxypristis cuspidata. Endanger Species Res. 2018;37:45–54.
Roycroft EJ, Port AL, Lavery SD. Population structure and male-biased dispersal in the short-tail stingray Bathytoshia brevicaudata (Myliobatoidei: Dasyatidae). Conserv Genet. 2019;20:717–28.
Smith K, Feldheim K, Carlson J, Wiley T, Taylor S. Female philopatry in smalltooth sawfish Pristis pectinata: conservation and management implications. Endanger Species Res. 2021;45:85–98.
Hueter RE, Heupel MR. Evidence of Philopatry in Sharks and Implications for the management of Shark Fisheries. North Am J Fish Manag. 2004;35:1–11.
Daly-Engel TS, Seraphin KD, Holland KN, Coffey JP, Nance HA, Toonen RJ, et al. Global phylogeography with mixed-marker analysis reveals male-mediated dispersal in the Endangered Scalloped Hammerhead Shark (Sphyrna lewini). PLoS ONE. 2012;7:e29986.
Portnoy DS, Heist EJ. Molecular markers: progress and prospects for understanding reproductive ecology in elasmobranchs. J Fish Biol. 2012;80:1120–40.
Momigliano P, Harcourt R, Robbins WD, Jaiteh V, Mahardika GN, Sembiring A et al. Genetic structure and signatures of selection in grey reef sharks (Carcharhinus amblyrhynchos). 2017;119:142–53.
Phillips NM, Devloo-Delva F, McCall C, Daly-Engel TS. Reviewing the genetic evidence for sex-biased dispersal in elasmobranchs. Rev Fish Biol Fisher. 2021;31:821–41.
Castillo-Páez A, Sosa-Nishizaki O, Sandoval-Castillo J, Galván-Magaña F, Blanco-Parra M-P, Rocha-Olivares A. Strong Population structure and shallow mitochondrial phylogeny in the Banded Guitarfish, Zapteryx exasperata (Jordan y Gilbert, 1880), from the Northern Mexican Pacific. J Hered. 2014;105:91–100.
Chevolot M, Ellis JR, Hoarau G, Rijnsdorp AD, Stam WT, Olsen JL. Population structure of the thornback ray (Raja clavata L.) in british waters. J Sea Res. 2006;56:305–16.
Flowers K, Ajemian M, Bassos-Hull K, Feldheim K, Hueter R, Papastamatiou Y, et al. A review of batoid philopatry, with implications for future research and population management. Mar Ecol Prog Ser. 2016;562:251–61.
Kitchen-Wheeler A-M, Ari C, Edwards AJ. Population estimates of Alfred mantas (Manta alfredi) in central Maldives atolls: North Male, Ari and Baa. Environ Biol Fish. 2012;93:557–75.
Germanov ES, Marshall AD. Running the Gauntlet: Regional Movement patterns of Manta alfredi through a Complex of Parks and Fisheries. PLoS ONE. 2014;9:e110071.
Knochel AM, Hussey NE, Kessel ST, Braun CD, Cochran JEM, Hill G, et al. Home sweet home: spatiotemporal distribution and site fidelity of the reef manta ray (Mobula alfredi) in Dungonab Bay, Sudan. Mov Ecol. 2022;10:22:1–17.
Setyawan E, Erdmann MV, Mambrasar R, Hasan AW, Sianipar AB, Constantine R, et al. Residency and use of an important Nursery Habitat, Raja Ampat’s Wayag lagoon, by juvenile reef Manta Rays (Mobula alfredi). Front Mar Sci. 2022;9:815094.
Armstrong AO, Stevens GMW, Townsend KA, Murray A, Bennett MB, Armstrong AJ, et al. Reef manta rays forage on tidally driven, high density zooplankton patches in Hanifaru Bay, Maldives. Peerj. 2021;9:e11992.
Doty MS, Oguri M. The island mass effect. ICES J Mar Sci. 1956;22:33–7.
Gove JM, McManus MA, Neuheimer AB, Polovina JJ, Drazen JC, Smith CR, Merrifield MA, Friedlander AM, Ehses JS, Young CW, Dillon AK, Williams GJ. Near-island biological hotspots in barren ocean basins. Nat Commun. 2016;7:10581:1–8.
Harris JL, Stevens GMW. Environmental drivers of reef manta ray (Mobula alfredi) visitation patterns to key aggregation habitats in the Maldives. PLoS ONE. 2021;16:e0252470.
Whitney JL, Gove JM, McManus MA, Smith KA, Lecky J, Neubauer P, Phipps JE, Contreras EA, Kobayashi DR, Asner GP. Surface slicks are pelagic nurseries for diverse ocean fauna. Sci Rep. 2021;11:3197.
Smith K, Whitney JL, Lecky J, Gove JM, Copeland A, Kobayashi DR, McManus MA. Accumulation in ocean surface lines varies by physical driver in coastal Hawaiian waters. Cont Shelf Res. 2021;104558.
Baird RW, Gorgone AM, McSweeney DJ, Ligon AD, Deakos MH, Webster DL, et al. Population structure of island-associated dolphins: evidence from photo‐identification of common bottlenose dolphins (Tursiops truncatus) in the main Hawaiian Islands. Mar Mammal Sci. 2009;25:251–74.
Martien KK, Baird RW, Hedrick NM, Gorgone AM, Thieleking JL, McSweeney DJ, et al. Population structure of island-associated dolphins: evidence from mitochondrial and microsatellite markers for common bottlenose dolphins (Tursiops truncatus) around the main Hawaiian Islands. Mar Mammal Sci. 2011;28:E208–32.
Andrews K, Karczmarski L, Au WWL, Rickards SH, Vanderlip CA, Bowen BW, et al. Rolling stones and stable homes: social structure, habitat diversity and population genetics of the Hawaiian spinner dolphin (Stenella longirostris). Mol Ecol. 2010;19:732–48.
Baird RW, Webster DL, Mahaffy SD, McSweeney DJ, Schorr GS, Ligon AD. Site fidelity and association patterns in a deep-water dolphin: Rough‐toothed dolphins (Steno bredanensis) in the Hawaiian Archipelago. Mar Mammal Sci. 2008;24:535–53.
Albertson GR, Baird RW, Oremus M, Poole MM, Martien KK, Baker CS. Staying close to home? Genetic differentiation of rough-toothed dolphins near oceanic islands in the central Pacific Ocean. Conserv Genet. 2017;18:33–51.
Courbis S, Baird RW, Cipriano F, Duffield D. Multiple populations of Pantropical spotted Dolphins in Hawaiian Waters. J Hered. 2014;105:627–41.
Lassauce H, Chateau O, Erdmann MV, Wantiez L. Diving behavior of the reef manta ray (Mobula alfredi) in New Caledonia: more frequent and deeper night-time diving to 672 meters. PLoS ONE. 2020;15:e0228815.
Setyawan E, Stevenson BC, Erdmann MV, Hasan AW, Sianipar AB, Mofu I, et al. Population estimates of photo-identified individuals using a modified POPAN model reveal that Raja Ampat’s reef manta rays are thriving. Front Mar Sci. 2022;9:1014791.
Armstrong AO, Armstrong AJ, Bennett MB, Richardson AJ, Townsend KA, Dudgeon CL. Photographic identification and citizen science combine to reveal long distance movements of individual reef manta rays Mobula alfredi along Australia’s east coast. Mar Biodivers Rec. 2019;12:14.
Couturier LIE, Jaine FRA, Townsend KA, Weeks SJ, Richardson AJ, Bennett MB. Distribution, site affinity and regional movements of the manta ray, Manta alfredi (Krefft, 1868), along the east coast of Australia. Mar Freshw Res. 2011;62:628–37.
Setyawan E, Erdmann MV, Lewis SA, Mambrasar R, Hasan AW, Templeton S, et al. Natural history of manta rays in the Bird’s Head Seascape, Indonesia, with an analysis of the demography and spatial ecology of Mobula alfredi (Elasmobranchii: Mobulidae). J Open Sci Foundation. 2020;36:49–83.
Venables S, van Duinkerken D, Rohner C, Marshall A. Habitat use and movement patterns of reef manta rays Mobula alfredi in southern Mozambique. Mar Ecol Prog Ser. 2020;634:99–114.
Hirschfeld M, Dudgeon C, Sheaves M, Barnett A, MacNeil A. Barriers in a sea of elasmobranchs: from fishing for populations to testing hypotheses in population genetics. Global Ecol Biogeogr. 2021;30:2147–63.
Vignaud T, Clua E, Mourier J, Maynard J, Planes S. Microsatellite analyses of blacktip reef sharks (Carcharhinus melanopterus) in a fragmented environment show structured clusters. PLoS ONE. 2013;8:e61067.
Pazmiño DA. Genome-wide SNPs reveal low effective population size within confined management units of the highly vagile Galapagos shark (Carcharhinus galapagensis). Conserv Genet. 2017;18:1151–63.
Arauz R, Chávez EJ, Hoyos-Padilla EM, Marshall AD. First record of the reef manta ray, Mobula alfredi, from the eastern Pacific. Mar Biodivers Rec. 2019;12:3.
Stewart JD, Beale CS, Fernando D, Sianipar AB, Burton RS, Semmens BX, et al. Spatial ecology and conservation of Manta birostris in the Indo-Pacific. Biol Conserv. 2016;200 C:178–83.
Harris JL, McGregor PK, Oates Y, Stevens GMW. Gone with the wind: seasonal distribution and habitat use by the reef manta ray (Mobula alfredi) in the Maldives, implications for conservation. Aquat Conserv Mar Freshw Ecosyst. 2020;30:1649–64.
Whiteley AR, Coombs JA, Hudy M, Robinson Z, Nislow KH, Letcher BH. Sampling strategies for estimating brook trout effective population size. Conserv Genet. 2012;13:625–37.
Waples RS, Antao T, Luikart G. Effects of overlapping generations on linkage disequilibrium estimates of Effective Population size. Genetics. 2014;197:769–80.
Deakos MH. Ecology and social behavior of resident manta ray (Manta alfredi) population off Maui, Hawai‘i. PhD Dissertation, University of Hawaii, Honolulu, USA; 2010.
Palstra FP, Fraser DJ. Effective/census population size ratio estimation: a compendium and appraisal. Ecol Evol. 2012;2:2357–65.
Dudgeon CL, Ovenden JR. The relationship between abundance and genetic effective population size in elasmobranchs: an example from the globally threatened zebra shark Stegostoma fasciatum within its protected range. Conserv Genet. 2015;16:1443–54.
Kashiwagi T. Conservation biology and genetics of the largest living rays: manta rays. Australia: The University of Queensland; 2014. PhD Dissertation.
Marshall AD. Holmberg. Manta Matcher Photo-identification Library. https://www.mantamatcher.org.
Frankham R. Conservation Genetics. Annu Rev Genet. 1995;29:305–27.
Hedgecock D. Does variance in reproductive success limit effective population sizes of marine organisms? Genetics and evolution of aquatic organisms. London: Chapman & Hall; 1994. 122–34.
Waples RS, Luikart G, Faulkner JR, Tallmon DA. Simple life-history traits explain key effective population size ratios across diverse taxa. Proc Royal Soc B Biological Sci. 2013;280:20131339.
Waples RS, Do C. Ldne: a program for estimating effective population size from data on linkage disequilibrium. Mol Ecol Resour. 2008;8:753–6.
Tallmon D, Gregovich D, Waples R, Baker CS, Jackson J, Taylor BL, et al. When are genetic methods useful for estimating contemporary abundance and detecting population trends? Mol Ecol Resour. 2010;10:684–92.
Palstra FP, Ruzzante DE. Genetic estimates of contemporary effective population size: what can they tell us about the importance of genetic stochasticity for wild population persistence? Mol Ecol. 2008;17:3428–47.
Portnoy DS, McDowell JR, McCandless CT, Musick JA, Graves JE. Effective size closely approximates the census size in the heavily exploited western Atlantic population of the sandbar shark, Carcharhinus plumbeus. Conserv Genet. 2008;10:1697.
Blower D, Pandolfi J, Bruce B, Gomez-Cabrera M, Ovenden J. Population genetics of australian white sharks reveals fine-scale spatial structure, transoceanic dispersal events and low effective population sizes. Mar Ecol Prog Ser. 2012;455:229–44.
Homma K, Maruyama T, Itoh T, Ishihara H, Uchida S. Biology of manta ray, Manta birostris Walbaum, in the Indo-Pacific. In: S´eret, B., Sire, editors J-Y, editors. Proceedings of the Fifth Indo-Pacific Fish Conference. Noum´ea, 1997. Soci´et´e Francaise d’Ichthyologie; 1999. p. 209–16.
Peel L, Stevens G, Daly R, Daly CK, Lea J, Clarke C, et al. Movement and residency patterns of reef manta rays Mobula alfredi in the Amirante Islands, Seychelles. Mar Ecol Prog Ser. 2019;621:169–84.
Notarbartolo-di-Sciara G. Natural history of the rays of the genus Mobula in the Gulf of California. Fish Bull. 1988;86:45–6.
Hoelzel AR, Shivji MS, Magnussen J, Francis MP. Low worldwide genetic diversity in the basking shark (Cetorhinus maximus). Biol Lett. 2006;2:639–42.
Ahonen H, Harcourt RG, Stow AJ. Nuclear and mitochondrial DNA reveals isolation of imperilled grey nurse shark populations (Carcharias taurus). Mol Ecol. 2009;18:4409–21.
Feutry P, Kyne PM, Pillans RD, Chen X, Naylor GJ, Grewe PM. Mitogenomics of the Speartooth Shark challenges ten years of control region sequencing. BMC Evol Biol. 2014;14:373–9.
Karl SA, Castro ALF, Garla RC. Population genetics of the nurse shark (Ginglymostoma cirratum) in the western Atlantic. Mar Biol. 2012;159:489–98.
Martin AP, Naylor GJP, Palumbi SR. Rates of mitochondrial DNA evolution in sharks are slow compared with mammals. Nature. 1992;357:153–5.
Spielman D, Brook BW, Frankham R. Most species are not driven to extinction before genetic factors impact them. Proc Natl Acad Sci. 2004;101:15261–4.
Teixeira JC, Huber CD. The inflated significance of neutral genetic diversity in conservation genetics. Proc Natl Acad Sci. 2021;118:e2015096118.
Jakimska A, Konieczka P, Skóra K, Namieśnik J. Bioaccumulation of metals in tissues of marine animals, Part I: the role and impact of heavy metals on organisms. Pol J Environ Stud. 2011;20:1117–25.
Germanov ES, Bejder L, Chabanne DBH, Dharmadi D, Hendrawan IG, Marshall AD, et al. Contrasting Habitat Use and Population Dynamics of reef Manta Rays within the Nusa Penida Marine Protected Area, Indonesia. Front Mar Sci. 2019;6:215.
Rohner C, Pierce S, Marshall A, Weeks S, Bennett M, Richardson A. Trends in sightings and environmental influences on a coastal aggregation of manta rays and whale sharks. Mar Ecol Prog Ser. 2013;482:153–68.
Woodworth-Jefcoats PA, Polovina JJ, Drazen JC. Climate change is projected to reduce carrying capacity and redistribute species richness in North Pacific pelagic marine ecosystems. Global Change Biol. 2017;23:1000–8.
Murray A, Garrud E, Ender I, Lee-Brooks K, Atkins R, Lynam R, et al. Protecting the million-dollar mantas; creating an evidence-based code of conduct for manta ray tourism interactions. J Ecotourism. 2020;19:132–47.
Dudgeon CL, Pollock KH, Braccini JM, Semmens JM, Barnett A. Integrating acoustic telemetry into mark–recapture models to improve the precision of apparent survival and abundance estimates. Oecologia. 2015;178:761–72.
DiBattista JD, Feldheim KA, Garant D, Gruber SH, Hendry AP. Anthropogenic disturbance and evolutionary parameters: a lemon shark population experiencing habitat loss. Evol Appl. 2011;4:1–17.
White WT, Giles J, Dharmadi, Potter IC. Data on the bycatch fishery and reproductive biology of mobulid rays (Myliobatiformes) in Indonesia. Fish Res. 2005;82:65–73.
Deakos MH. Paired-laser photogrammetry as a simple and accurate system for measuring the body size of free-ranging manta rays Manta alfredi. Aquat Biology. 2010;10:1–10.
Toonen RJ, Puritz JB, Forsman ZH, Whitney JL, Fernandez-Silva I, Andrews KR, et al. ezRAD: a simplified method for genomic genotyping in non-model organisms. Peerj. 2013;1:e203.
Knapp I, Puritz J, Bird C, Whitney J, Sudek M, Forsman Z et al. ezRAD- an accessible next-generation RAD sequencing protocol suitable for non-model organisms_v3.2. protocols.io. 2016. https://doi.org/10.17504/protocols.io.bszgnf3w.
Fisher S, Barry A, Abreu J, Minie B, Nolan J, Delorey TM, et al. A scalable, fully automated process for construction of sequence-ready human exome targeted capture libraries. Genome Biol. 2011;12:R1.
Tisthammer KH, Forsman ZH, Sindorf VL, Massey TL, Bielecki CR, Toonen RJ. The complete mitochondrial genome of the lobe coral Porites lobata (Anthozoa: Scleractinia) sequenced using ezRAD. Mitochondrial Dna Part B. 2016;1:247–9.
Antaky CC, Kitamura PK, Knapp IS, Toonen RJ, Price MR. The complete mitochondrial genome of the Band-rumped storm petrel (Oceanodroma castro). Mitochondrial Dna Part B. 2019;4:1271–2.
Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25:1754–60.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Barnett DW, Garrison EK, Quinlan AR, Strömberg MP, Marth GT. BamTools: a C + + API and toolkit for analyzing and managing BAM files. Bioinformatics. 2011;27:1691–2.
Garrison E, Marth G. Haplotype-based variant detection from short-read sequencing. Preprint at arXiv. 2012. https://doi.org/10.48550/arXiv.1207.3907.
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.
Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, et al. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34:3299–302.
Meirmans PG, van Tienderen PH. Genotype and genodive: two programs for the analysis of genetic diversity of asexual organisms. Mol Ecol Notes. 2004;4:792–4.
Tajima F. Evolutionary relationship of DNA sequences in finite populations. Genetics. 1983;105:437–60.
Fu Y-X. Statistical tests of Neutrality of mutations against Population Growth, hitchhiking and background selection. Genetics. 1997;147:915–25.
Excoffier L, Lischer HE. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10:564–7.
Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9:772–2.
NETWORK Software. http://www.fluxus-engineering.com/network_terms.htm. Accessed 1 May 2021.
Bandelt HJ, Forster P, Rohl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Andrews S. FastQC: A Quality Control Tool for High Throughput Sequence Data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. 2010.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnetJournal. 2011;17:1–3.
Belclaid EC. M. Seanome. https://github.com/UH-Bioinformatics/Seanome.git.
Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: a fast and accurate Illumina paired-end read merger. Bioinformatics. 2014;30:614–20.
Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26:2460–1.
Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: Open-Source, Platform-Independent, community-supported Software for describing and comparing Microbial Communities. Appl Environ Microb. 2009;75:7537–41.
Puritz JB, Hollenbeck CM, Gold JR. dDocent: a RADseq, variant-calling pipeline designed for population genomics of non-model organisms. Peerj. 2014;2:e431.
Broad. Picard Tools. Broad Institute, GitHub repository. http://broadinstitute.github.io/picard/; 2013.
O’Leary SJ, Puritz JB, Willis SC, Hollenbeck CM, Portnoy DS. These aren’t the loci you’e looking for: principles of effective SNP filtering for molecular ecologists. Mol Ecol. 2018;27:3193–206.
Gabor EG. M. VcfLib v.1.0.3. https://github.com/vcflib/vcflib.
Willis SC, Hollenbeck CM, Puritz JB, Gold JR, Portnoy DS. Haplotyping RAD loci: an efficient method to filter paralogs and account for physical linkage. Mol Ecol Resour. 2017;17:955–65.
Lischer HEL, Excoffier L. PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics. 2012;28:298–9.
Yang J, Benyamin B, McEvoy BP, Gordon S, Henders AK, Nyholt DR, et al. Common SNPs explain a large proportion of the heritability for human height. Nat Genet. 2010;42:565–9.
Whitlock MC, Lotterhos KE. Reliable detection of loci responsible for local adaptation: inference of a null model through trimming the distribution of FST. Am Nat. 2015;186:24–36.
OutFLANK Software. https://github.com/whitlock/OutFLANK.
Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics. 2001;29:1165–88.
Narum SR. Beyond Bonferroni: less conservative analyses for conservation genetics. Conserv Genet. 2006;7:783–7.
Pritchard JK, Stephens M, Donnelly P. Inference of Population structure using multilocus genotype data. Genetics. 2000;155:945–59.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.
Earl DA, vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4:359–61.
Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I. Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour. 2015;15:1179–91.
Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23:1801–6.
Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5.
Nomura T. Estimation of effective number of breeders from molecular coancestry of single cohort sample. Evol Appl. 2008;1:462–74.
Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR. NeEstimatorv2: re-implementation of software for the estimation of contemporary effective population size (ne) from genetic data. Mol Ecol Resour. 2014;14:209–14.
Beerli P. Comparison of bayesian and maximum-likelihood inference of population genetic parameters. Bioinformatics. 2006;22:341–5.
Beerli P, Felsenstein J. Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. Proc Natl Acad Sci. 2001;98:4563–8.
Keenan K, McGinnity P, Cross TF. diveRsity: an R package for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol Evol. 2013;4:782–8.
Sundqvist L, Keenan K, Zackrisson M, Prodöhl P, Kleinhans D. Directional genetic differentiation and relative migration. Ecol Evol. 2016;6:3461–75.
Jost L. GST and its relatives do not measure differentiation. Mol Ecol. 2008;17:4015–26.
Special thanks to Lee James with Ultimate Whale Watch Adventure, Craig Venema, and Amy Miller for providing a boat and valuable time in the field collecting data and retrieving satellite tags. Thanks to Keller Laros with the Manta Pacific Research Foundation for assistance with the West Hawaii field work. Thanks to Doug Perrine and Bob Gladen who donated their boats and their time during field work in Kona. We thank Joey Lecky for his assistance with graphics in Fig. 1. Thanks to Stephen Karl, Brian Bowen, Rob Toonen, Jamison Gove, and Ryan Rykaczewski for their support of this research. Thanks to the State of Hawai‘i Division of Aquatic Resources especially Catherine Gewecke. Special thanks to Patricia Rosel, Matthew Craig, and Chelsey Young for reviewing earlier drafts of the manuscript.
This research was supported by donations from Sven and Krista Lindblad to the Hawaiʻi Association for Marine Education and Research, the Cooperative Institute for Marine and Atmospheric Research (CIMAR), NOAA’s West Hawai‘i Integrated Ecosystem Assessment Program, the University of Texas at Austin’s Provosts Early Career Fellowship, and NOAA’s Pacific Islands Fisheries Science Center. The funding bodies played no role in the design of the study and collection, analysis, interpretation of data, and in writing the manuscript.
The authors declare no competing interests.
Ethics approval and consent to participate
All methods were performed in accordance with institutional, state, and national guidelines and regulations regarding the collection of biopsy samples and were approved by the University of Hawai‘i Animal Care & Use Committee, Protocol No. 08-591-2, and Assurance number A3423-01. This study is reported in accordance with ARRIVE guidelines (https://arriveguidelines.org). This study complied with the International Union for Conservation of Nature (IUCN) policies for research involving species that are threatened or at risk of extinction, the Convention on Biological Diversity and the Convention on the Trade in Endangered Species of Wild Fauna and Flora.
Consent for publication
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
Supplementary tables S1-S6.
(1) Step-by-step protocol for library prep benchwork and (2) Detailed steps of SNP filtering workflow
Mitogenome Assembly & Alignment Summary Tables
About this article
Cite this article
Whitney, J.L., Coleman, R.R. & Deakos, M.H. Genomic evidence indicates small island-resident populations and sex-biased behaviors of Hawaiian reef Manta Rays. BMC Ecol Evo 23, 31 (2023). https://doi.org/10.1186/s12862-023-02130-0
- Hawaiian Islands
- Devil rays
- Genetic connectivity