Exploring the genetic diversity and population structure of Mobula birostris in two key aggregation zones in the Eastern Tropical Pacific

: The giant manta ray Mobula birostris is the largest ray species in the world. Little is known about its genetic composition in key aggregation sites such as the Galapagos Islands and Isla de la Plata, near the province of Manabi in mainland Ecuador. This study aimed to determine the genetic diversity and population structure of M. birostris in these 2 locations to better understand its connectivity and distribution in Ecuadorian oceanic waters and to assist in its conservation and appropriate management. A total of 127 samples from mainland Ecuador (2013−2018) and 21 samples from Galapagos (2019) were collected and analyzed using 8 microsatellite loci. Results showed a moderately high level of genetic diversity for giant manta rays from both sites (mainland Ecuador H e = 0.72; Galapagos H e = 0.66). Population structure analyses suggests the presence of 2 different populations in the Galapagos and mainland Ecuador. The different genetic compositions found for each location could be associated with the displayed resident behavior, linked to the formation of upwelling systems caused by oceanic currents that bring nutrient-rich waters to both sites year-round. Our genetic connectivity analysis confirmed low gene flow between these 2 locations, further rejecting the hypothesis of a single panmictic population of M. birostris in Ecuador. Taken together, these results provide valuable information about the genetic composition and diversity of the giant manta ray, an Endangered species which has been scarcely studied in the Eastern Tropical Pacific.


INTRODUCTION
The Mobulidae family comprises all manta and devil rays distributed globally throughout tropical and temperate oceanic waters (Couturier et al. 2015). Mobulids are large pelagic filter-feeding rays whose diet is mainly composed of planktonic animals (Burgess 2017, Bray 2018). These species are characterized by low reproduction rates (one pup every 1−3 yr) as well as long maturation times (3−6 yr for Mobula alfredi, 10 yr for M. mobular, unknown for M. birostris) (Dulvy et al. 2014, Stewart et al. 2018, characteristics that lead to slow recovery rates and increased vulnerability for the group (Marshall et al. 2011a,b, O'Malley et al. 2013. The main threats for mobulids include targeted and bycatch fishing as well as habitat degradation and climate change (Couturier et al. 2012).
The giant manta ray M. birostris is the largest mobulid species, with individuals that measure up to 7 m in disc width and weigh up to 2 t (Marshall et al. 2009, McClain et al. 2015. This species shows a preference for pelagic, offshore habitats in subtropical oceanic waters (Kashiwagi et al. 2011). Due to its large size and seasonal sightings, M. birostris has been proposed to be a highly nomadic and migratory species (Hearn et al. 2014, Couturier et al. 2015. However, site fidelity has been reported in specific locations such as the Revillagigedo Islands in Mexico and the Raja Ampat Archipelago in Indonesia (Stewart et al. 2016). The first studies conducted on M. birostris (Clark 2002, Marshall 2008, Dewar et al. 2008 pre-date its taxonomic separation from M. alfredi (Marshall et al. 2009), so in these publications, the information on these 2 species could be misleading. However, in this paper, when we refer to giant manta rays (M. birostris), we have only used studies published after the taxonomic separation.
The attraction of M. birostris to highly productive tropical and subtropical areas where some targeted species such as tuna aggregate, coupled with the giant manta ray's large size and distribution in epipelagic zones, have made it vulnerable to bycatch in fisheries (Stewart et al. 2018). This has led to population declines and the categorization of the species as Endangered in the IUCN Red List of Threatened Species (IUCN 2020). Although important conservation steps have been made, such as the inclusion of M. birostris in Appendix II of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) (Marshall et al. 2011a,b), research is still needed to understand population structure and dynamics at local and regional scales, especially in key aggregation sites such as Isla de La Plata near the Ecuadorian coastline (Dulvy et al. 2014, Hearn et al. 2014, Kashiwagi et al. 2015 and the Galapagos Islands (Hearn et al. 2014, Lezama-Ochoa et al. 2019a. In highly mobile marine species, the identification of aggregation sites has proven useful for understanding critical biological processes (Derous et al. 2007, Lieber et al. 2020. In mobulids, aggregation sites have been linked to areas with high primary productivity, where individuals congregate mainly to forage and reproduce (Armstrong et al. 2016, Stewart et al. 2018. Characterization of the populations that occur in these areas could provide insights into the giant manta ray's distribution, movements, and behavior (Stewart et al. 2018).
Genetic studies of M. birostris are limited worldwide (Stewart et al. 2016, Hosegood 2020, and none of these reports include samples from Ecuadorian oceanic waters. World information about giant manta rays has been predominantly gathered through diver counts, photo-identification surveys, and acoustic and satellite tags of targeted populations at aggregation sites (Stewart et al. 2018). Although this data has contributed to an understanding of the population dynamics of the species, genetic analyses are still needed to assess population structure and identify areas of overlap between mobulid hotspots and fisheries, which could ultimately help reduce their bycatch rates (Stewart et al. 2018). Moreover, genetic information (e.g. reduced genetic diversity) could provide insights into the risks and vulnerability associated with inbreeding or demographic changes and contribute to creating efficient conservation plans based on genetic assessments rather than directing efforts based on political or geographic boundaries (Milton & Shaklee 1987, Allendorf et al. 2010, Bestervan der Merwe & Gledhill 2015, Domingues et al. 2018, Stewart et al. 2018. In this sense, evaluating the genetic diversity of M. birostris could help prioritize vulnerable populations for conservation and management actions (Stewart et al. 2018).
Population structure and genetic diversity in marine environments can vary temporally due to several factors, such as climatic change, inner habitat characteristics, oceanic currents, geographic distances, biotic interactions, and early life history (Reynolds et al. 2017, Lassauce et al. 2022. Im portantly, giant manta rays display site fidelity linked to resource availability, which is another key process that can shape population structure by limiting migration rates (Stewart et al. 2016, Lezama-Ochoa et al. 2019a. Specifically, planktivores rely on conditions influencing their environment, which shape the abundance and availability of their food sources (Burgess 2017, Lassauce et al. 2022. This, in turn, has been hypothesized as a driver for giant manta ray movements and migration (Stewart et al. 2016). Understanding how these factors may affect the populations present off the coast of mainland Ecuador and in the Galapagos is key to comprehending the ecology and behavior of the individuals present at these sites.
Given the favorable oceanic conditions (i.e. high productivity) during the summer off mainland Ecuador (July to October) and in the Galapagos (December to May) and the ability of M. birostris to perform long-distance migration (Marshall et al. 2011a, Hearn et al. 2014, Andrzejaczek et al. 2021, we hypothesized that M. birostris individuals could migrate between these 2 locations during seasons of high primary productivity and that they might constitute a panmictic population. Therefore, the objective of this study was to assess the genetic diversity, population structure, and connectivity of M. birostris individuals sampled off mainland Ecuador and in Galapagos throughout the course of several sighting seasons.

Sample collection and DNA extraction
A total of 148 Mobula birostris muscle/skin tissue samples were collected from free-swimming giant manta rays, using a biopsy tip mounted on a spear pole while SCUBA or free diving. Of these samples, 127 were collected off the coast of mainland Ecuador during the 2013−2018 sighting seasons from August to September (2013: n = 24;2014: n = 26;2015: n = 26;2016: n = 11;2017: n = 14;2018: n = 26); 121 were collected at Isla de la Plata (located ~30 km from Manabi Province in mainland Ecuador) and 6 were collected from Bajo Cope (~65 km from Isla de la Plata). Both locations (Isla de la Plata and Bajo Cope) are hereafter referred to collectively as mainland Ecuador. The remaining 21 samples were collected from the Galapagos during the 2019 sighting season from March to June. Samples were stored in 96% ethanol until DNA extraction.
Sampling details and collection sites are provided in Table S1 in the Supplement at www.int-res.com/ articles/suppl/m699p075_supp.pdf. Genomic DNA was extracted using a Proteinase K protocol described by Broderick et al. (2011). Maps of the georeferenced sampling sites on mainland Ecuador and Galapagos were drawn using ArcGIS Desktop 10.2 (ESRI 2020) (see Fig. 1).

Simple-sequence repeat amplification and genotyping
A total of 8 simple-sequence repeat (SSR) polymorphic microsatellite markers designed for M. alfredi (Kashiwagi et al. 2012) were used in this study. PCR reactions had a final volume of 25 μl, which included 1× PCR buffer (Invitrogen), 1.5−2.5 mM MgCl 2 (Invitrogen), 0.4 mM dNTPs (Invitrogen), 0.15 μM of tail A modified forward primer, 0.5 μM of reverse primer, 0.2 μM of tail A fluorescent primer, 0.1−1 mg ml −1 BSA, 1 U of Taq DNA polymerase (Invitrogen), and 10−20 ng of DNA. The thermocycling conditions consisted of an initial denaturation step at 94°C for 5 min, followed by 37 cycles (30 s at 94°C, 45 s at 55−65°C, 72°C for 1.5 min) with a final extension step of 10 min at 72°C (Kashiwagi et al. 2012 with modifications). PCR products from the 2014−2019 seasons were labeled with one of 4 fluorescent dyes (6-FAM, VIC, PET, or NED), using universal primers in a 3primer system described by Blacket et al. (2012). Labeled amplified PCR products were genotyped by Macrogen on an ABI 3130 Genetic Analyzer (Ap plied Biosystems) automatic capillary sequencer, using 500LIZ as a size standard. Genotyping results were analyzed using the software GeneMarker v.2.4.0 (Softgenetics).
Samples for the 2013 season were genotyped via polyacrylamide gel electrophoresis: PCR products were separated using 6% polyacrylamide gels, and gels were revealed using the silver staining protocol as described by Benbouza et al. (2006). Amplicon size for each allele was calculated using linear regression with a 100 bp ladder (Invitrogen). For standardizing the allele sizes obtained by polyacrylamide electrophoresis with the ones obtained by fluorescence genotyping, 15 bp corresponding to the universal 5' tail was added to the final amplicon size.

Relatedness analysis
Relatedness between individuals and within seasons was examined using the R package 'related' v.1.0 (Pew et al. 2015) in R v.3.6.2 (R Core Team 2019). This program calculates relatedness from micro satellites using 7 different estimators (Pew et al. 2015). To test the performance of the estimators, we calculated the Pearson correlation coefficient between observed and expected relatedness for each case. Based on these results, we chose to use the Wang estimator for relatedness analysis (Wang 2002).
To account for the possibility that the same individual might have been resampled multiple times within a season, we performed pairwise comparisons among individuals and then randomly eliminated one individual from each pairwise comparison that showed high relatedness values (R > 0.7). Based on the allele frequencies of the data set, we performed simulations of 100 individuals with a predefined degree of relatedness. We then calculated the pairwise relatedness values among all sampled individuals and clas-sified them into specific relatedness groups. Relatedness indexes among seasons were statistically compared by 1-way ANOVA. We also conducted 1000 Monte Carlo simulations (1000 iterations) to calculate the expected relatedness distribution in each season.

Genetic diversity estimation
Genetic diversity parameters such as the number of alleles per locus (N a ), polymorphic information content ( (Raymond & Rousset 1995).

Evaluation of population structure
F ST was calculated by running an analysis of molecular variance (AMOVA) in ARLEQUIN v.3.5.2.2 (Excoffier & Lischer 2010) with significance tested against 10 000 permutations. The AMOVA significance level was adjusted to α = 0.05. A principal coordinate analysis (PCoA) was performed using the 'ape' v.5.4-1 package (Paradis & Schliep 2019) and plotted with ggplot2 v.3.3.3 (Wickham 2009) to visualize the genetic structure between Galapagos and mainland Ecuador and among seasons within mainland Ecuador. Population structure was inferred through a Bayesian individual-based clustering ap proach using STRUC-TURE v.2.3.4 (Pritchard et al. 2000). The number of clusters was determined using an admixture model, where K values between 1−10 were evaluated, considering 5 iterations for each K value. Each run consisted of 1 000 000 Markov chain Monte Carlo (MCMC) iterations with a 100 000 step burn-in period. The optimal K value was obtained using the Evanno method (Evanno et al. 2005) implemented in STRUC-TURE HARVESTER (Earl & von Holdt 2012), and information corresponding to all in dividual runs was compiled using CLUMPP (Jakobsson & Rosenberg 2007). The final STRUCTURE plot was obtained using the software 'Distruct' (Rosenberg 2003).

Migration analysis
To investigate directional gene flow and its relative magnitude among the different seasons of M. biro stris in mainland Ecuador and the Galapagos, we used the R package 'diveRsity' v.1.9.90 (Keenan et al. 2013). This method computes migration rates be tween all sites and then normalizes them to obtain re lative migration rates that vary between 0 and 1. Computations were performed using effective number of migrants (N mAlcala ), which is a statistic introduced by Alcala et al. (2014) (Sundqvist et al. 2016). This statistic uses the genetic differentiation measures G ST and D to estimate N mAlcala , from which the relative migration rate is calculated.

Compilation of environmental oceanographic data
To further describe the variables associated with both sampling locations, information corresponding to oceanographic conditions during the sampling periods was collected from 2 online databases. The data was downloaded from the EU Copernicus Marine Environment Monitoring Service (CMEMS; http:// marine.copernicus.eu) and Environmental Research Division's Data Access Program (ERDDAP; https:// www.ncei.noaa.gov/erddap/index.html). CMEMS is a platform that gathers information from multiple databases on the ocean's state, variability, and dynamics from marine ecosystems around the world (von Schuck mann et al. 2016, Lezama-Ochoa et al. 2019b). ERDDAP is a data server that grants access to scientific data sets on oceanographic data collected from satellites and buoys around the world (Simons 2020). The following environmental variables were downloaded: daily sea surface temperature (SST; °C), monthly chlorophyll a (chl a; mg m −3 ), and monthly oxygen concentration (O 2 ; mmol m −3 ). All variables were collected at a 0.25° spatial resolution using python routines and motuclient from CMEMS and a 0.025° spatial resolution from ERDDAP. These variables were selected since they are considered indicators of primary productivity levels associated with aggregation areas for giant manta rays (Lezama-Ochoa et al. 2019a).

Relatedness
A total of 148 Mobula birostris samples were collected and genotyped from mainland Ecuador (n = 127) during 6 consecutive seasons (2013−2018), and from southwest Galapagos (n = 21) during one season (2019) (Fig. 1). In order to avoid biased estimations of genetic diversity and population structure, a relatedness analysis was conducted to exclude related individuals from further analysis. The results from this analysis showed 2 pairs of duplicate individuals (R = 1) and 10 pairs of highly related individuals (R > 0.7). One individual from each pair was randomly removed from the data set (4 from 2013, 1 from 2014, 1 from 2015, 1 from 2016, 2 from 2018, and 3 from 2019). Therefore, all remaining analyses were performed on a data set of 136 individuals (118 from mainland Ecuador and 18 from Galapagos).
Monte Carlo simulations indicated that the relatedness observed for individuals for the 2013, 2014, 2018, and 2019 seasons was significantly higher (p < 0.001) than would be expected by random mating, meaning that individuals from these seasons are more closely related than expected by chance (Fig. S1). According to the Monte Carlo simulations, giant manta rays from the same seasons (2013,2014,2018, and 2019) presented a higher median Wang relatedness value (middle line in the box plot in Fig. S2) compared to individuals from 2015, 2016, and 2017. Significant differences (p < 0.05) in relatedness indexes among the aforementioned seasons were found by 1way ANOVA.
Analysis among seasons indicated that the Wang median value was higher when considering only individuals from mainland Ecuador (seasons 2013− 2018) than when comparing all individuals from both localities (2013−2019). This finding indicates that there was a higher number of related individuals in the samples collected only in mainland Ecuador than when samples from Galapagos were included (Fig. S2).

Genetic diversity
The data set of 136 individuals was genotyped at 8 microsatellite loci. The mean N o frequency was 10.13%. Paired t-tests indicated no significant difference (p > 0.05) between uncorrected and corrected F ST values, which indicates that N o likely had no effect on the genetic analyses. All 8 loci were used in subsequent analyses. Allelic dropout was distributed between 0 and 20% across the 8 analyzed loci (0%: MA14; 20%: MA15). All markers deviated significantly (p < 0.05) from HWE by having less H o than expected, except for locus MA09 (p = 0.64) (Sharma et al. 2016, Chen et al. 2017. A total of 94 alleles, ranging from 7−25 alleles per locus, were found in the data set. Furthermore, the PIC indicated that MA30 was the most informative marker, with a total of 25 alleles and a PIC value of 0.90. H e ranged from 0.49 (MA21) to 0.90 (MA30) per locus, with a mean of 0.74. A summary of the overall genetic diversity per locus across the 136 individuals is provided in Table 1.
Overall genetic diversity varied among the 6 seasons in mainland Ecuador. The average N a per season ranged from 47 (2013) to 58 (2015), and N pa ranged from 1 allele (2016)  Allelic richness estimates were also higher for mainland Ecuador (7.15) than Galapagos (5) ( Table 2).

Population differentiation and population structure
Genetic differentiation based on pair wise F ST estimates between the 6 sighting seasons for mainland Ecuador indicated that the greatest differentia-     (Table 3). When comparing mainland Ecuador to Galapagos, the F ST index revealed significant genetic differentiation (F ST = 0.08, p < 0.05).
A PCoA based on genetic distances suggested that the first 2 principal components accounted for 10.35 and 7.48% of the genetic variation in the 136 M. birostris individuals, together explaining 17.83% of the total variation. The PCoA plot revealed 2 clusters, one grouping the Galapagos samples and the other grouping samples from mainland Ecuador. However, a few individuals were included in both clusters (Fig. 2). A PCoA performed only on the samples from mainland Ecuador comparing seasons (Fig. S3) suggested no clear population structure (2013−2018). Overall, the PCoA plotting analysis comparing mainland Ecuador and Galapagos illustrated a pattern consistent with the obtained pairwise F ST values, suggesting the possibility that these are 2 genetically different populations.
The STRUCTURE clustering analysis for all samples (2013−2019) tested K values ranging from 1−10 (Fig. S4). Results from K = 2 to K = 10 indicated differentiation between lineages from mainland Ecuador (2013−2018) and Galapagos (2019) (Fig. S4). In addition, the K = 3 and K = 5 graphs showed different genetic lineages between mainland individuals from the 2013 and 2018 seasons as well as the rest of the mainland Ecuador seasons (2014−2017) (Fig. S4), which was also supported by F ST values (Table 3). According to the Evanno method, the optimal ΔK value is 3, supporting the presence of 3 genetic lineages: 2 corresponding to mainland Ecuador (purple and blue bars in Fig. 3a) and one corresponding to Galapagos (orange bars in Fig. 3a).
The STRUCTURE analysis for the different seasons in mainland Ecuador (2013−2018) also tested K values from 1−10 (Fig. S5). In this case, the K = 2 and K = 4 graphs once again indicated a difference in the composition of the genetic lineages between individ-

Gene flow and relative migration between mainland Ecuador and Galapagos
The migration networks revealed different levels of migration rates and gene flow within 2013−2018 in mainland Ecuador, and between these years and the single sighting season (2019) in Galapagos. The relative migration network (Fig. 4) (Fig. S6). Together, these results (F ST indexes, PCoA, STRUCTURE, and migration analyses) support the presence of 2 separate populations: one in mainland Ecuador and one in the Galapagos, refuting the single panmictic population hypothesis.

Environmental characteristics at the sampling locations associated with giant manta ray aggregation
To understand how environmental characteristics could be contributing to the presence of 2 distinct populations in Ecuadorian waters, oceanographic variables were gathered for mainland Ecuador and Galapagos for the 2013−2019 seasons (  (Table 4). All ranges observed for these variables are congruent with conditions previously reported for the occurrence,

DISCUSSION
In this study, we examined the genetic composition of Mobula birostris in 2 important aggregation zones of the Eastern Tropical Pacific and shed light on its population structure and distribution -key factors for the conservation of this Endangered marine species. Based on F ST , population structure, and migration analyses, our results suggest the presence of 2 different populations of giant manta rays, with reduced/limited gene flow be tween mainland Ecuador and Galapagos.

Genetic status of M. birostris in 2 aggregation zones in the Eastern Pacific
Moderately high levels of genetic diversity were found in mainland Ecuador and the Galapagos (H e = 0.72 and H e = 0.66, respectively). To the best of our knowledge, no previous studies have evaluated the genetic diversity of M. birostris using SSR markers, but these values are comparable to or higher than those reported for other elasmobranchs, such as whale sharks collected worldwide (H e = 0.68) (Schmidt et al. 2009), white sharks from South Africa (H e = 0.66) (Pardini et al. 2000), and closely re lated species such as the reef manta ray from Japan (H e = 0.48) (Kashiwagi et al. 2012).
The only other study evaluating gen etic diversity in M. birostris was conducted by Hosegood (2020), who used single nucleotide polymorphism (SNP) markers from individuals sampled from 7 locations worldwide and reported low levels of genetic diversity (H e = 0.06). However, the data from Hosegood (2020) is not comparable to our study due to both the inner characteristics of the molecular markers used (SNPs: genome-wide analyzed loci; SSRs: specific analyzed loci) (Fischer et al. 2017) and the difference in sampling strategy: our sample size was larger (Hosegood 2020: n = 99; this study: n = 136) and we sampled manta rays from a much smaller geographic area that had not been previously studied.
The moderately high genetic diversity we found could be the result of the aggregation of a large number of individuals in the studied locations. Isla de la Plata off mainland Ecuador hosts one of the largest populations of M. birostris in the world (a total of 2400 individuals), compared to other populations reported in Pacific Mexico (n = 715) and Raja Ampat in Indonesia (n = 588) (Beale et al. 2019, Palomino et al. 2020. In general, larger populations tend to show higher levels of genetic diversity, mainly because of a higher likelihood of maintaining a robust genetic pool (Yun et al. 2020). In addition, the Galapagos and mainland Ecuador have been reported as important habitat and aggregation sites for this species (Vergara-Chen et al. 2015, Lezama-Ochoa et al. 2019a. Aggregation sites have been found to foster encounters and interbreeding between individuals that may come from genetically distinct stocks. Such events could help explain the increased genetic diversity indexes obtained for giant manta rays at these locations, as has been previously reported for green sea turtles in the Atlantic (Vásquez-Carrillo et al. 2020).
For species such as M. birostris, where population declines have been reported, it is important to determine parameters such as genetic diversity and connectivity (Beale et al. 2019, Sandoval-Castillo 2019, Hosegood 2020. A decreasing population usually loses diversity due to inbreeding, excessive genetic drift, and population fragmentation, which in turn    It is important to note the limitations of the present study with respect to estimation of genetic diversity and inference of population structure due to the number of samples analyzed in some seasons (2016 = 10 samples; 2017 = 14 samples) and because of the type and number of markers used (8 SSR markers). It has been suggested that at least 30 individuals per population are necessary for population genetic analysis (Putman & Carbone 2014, Danusevicius et al. 2016; however, this number cannot always be reached, especially when working with species such as M. birostris, where sampling is challenging. In ad dition, Putman & Carbone (2014) suggested that ideally > 50 SSRs would be needed to obtain good resolution of a species' population structure, whereas other reports state that correct assignment of individuals to their specific clusters can be obtained with as few as 8 markers (Arthofer et al. 2018). Nevertheless, further studies should consider increasing the number of loci analyzed and the number of samples per locality to obtain a more accurate estimate of the genetic diversity and population structure of this species.

Population structure and connectivity
A combination of F ST indexes, PCoA analysis, Bayesian population structure modeling, and migration networks suggest the presence of 2 discrete populations of M. birostris in mainland Ecuador and Galapagos (~1000 km apart). Similarly, Stewart et al. (2016) used SNP markers and revealed wellstructured giant manta ray subpopulations in Coastal Mexico (Bahía de Banderas), offshore Mexico (Revillagigedo Islands), and Sri Lanka. However, Hosegood (2020) reported a lack of differentiation and population structure between distant locations from around the world (Mexican Pacific, Mexican Caribbean, South Africa, Sri Lanka, Phillipines, and Peru), using SNPs. Differences between both studies can be explained by the contrasting analytic strategies used with their data sets. Importantly, the study by Hosegood (2020) did not include samples from Ecuador. The presence of discrete groups in a previously unstudied area highlights the need to fill geographic gaps and to understand the genetic structure of M. birostris globally.
Migration patterns of M. birostris in the Eastern Pacific have been analyzed recently using satellite tags and acoustic devices and suggest few exchanges of individuals and sporadic large-scale movements (Hearn et al. 2014, Palomino et al. 2020, Andrzejaczek et al. 2021. Palomino et al. (2020) tagged 46 giant manta rays at Isla de la Plata (near Manabi Province in mainland Ecuador) with acoustic devices (n = 30) and satellite tags (n = 16). Most of the giant manta rays displayed high residency within the study area (moving less than 95.3 km), whereas a few individuals migrated south towards Peruvian waters (387 km). More recently, Andrzejaczek et al. (2021) found that 2 tagged individuals moved north along the coast of Ecuador from Tumbes to other sites such as Isla de la Plata (263 km), and a third individual performed a large-scale movement from Peru to northern Ecuador and then traveled in a westerly direction to the Galapagos Islands (approximately 1300 km). Importantly, the satellite data for this third individual shows that it stayed within the Galapagos Marine Reserve for less than 1 mo in September 2018 before moving in a southeastern direction. This differs from local observations in southern Galapagos, where giant manta rays are commonly seen from December to June (D. Pazmiño pers. comm.). Hearn et al. (2014) also reported that one out of 9 tagged individuals moved from Isla de la Plata to the Galapagos Islands. Despite evidence of single individuals migrating be tween locations, the limited gene flow found in this study between mainland Ecuador and Galapagos (Figs. 4 & S6) is not surprising since large populations usually have low genetic drift and require very low migration rates to maintain genetic similarity even if they have become physically separated (Veríssimo et al. 2017, Marandel et al. 2018. The 2 tagged individuals that migrated between mainland Ecuador and Galapagos in 2012(Hearn et al. 2014, Andrzejaczek et al. 2021) could indicate genetic connectivity among these populations, since in order to establish this kind of connectivity it would be necessary to find at least one migrant for every 10 generations. For M. birostris, which has an approximate generation time of 25 yr, these 2 migration events could be sufficient to establish genetic connectivity between these 2 populations (Lowe & Allendorf 2010, Marshall et al. 2011a, Marandel et al. 2018. Notwithstanding, without knowing the reproductive status of these 2 marked individuals it would be difficult to establish any conclusions about the connectivity between these populations. Thus, the limited connectivity found in this study may be associated with a resident behavior of M. birostris at the sampling locations. Such behavior has also been re-ported for M. birostris in Mexican waters, where a defined population structure was found between Bahia de Banderas and Revillagigedo Island (distance ~ 600 km) (Stewart et al. 2016). Palomino et al. (2020) and Andrzejaczek et al. (2021) suggested that the residence patterns and spatial− temporal variations in vertical movement observed for giant manta rays moving from Peru to Ecuador and vice versa could be associated with high primary productivity (in southeast Pacific coasts) and higher zooplankton biomass during the austral spring months. Stewart et al. (2016) postulated that with a yearround food source, a suitable juvenile habitat of M. birostris overlapping with adult habitat might eliminate incentives for long-range migratory behavior. Other highly migratory elasmobranchs, such as tiger sharks, have shown philopatric and residency behavior in the Galapagos associated with year-round reliable food sources (Acuña-Marrero et al. 2017). Thus, favorable oceanographic conditions found in Ecu adorian waters can influence spatial patterns of residency be tween mainland Ecuador and Galapagos populations by providing suitable habitat and a continuous food source in both locations, and by reducing the need for long-distance travel to feed.
In mainland Ecuador, the cold Humboldt Current induces upwelling of cold, nutrient-rich waters throughout the year. In the Galapagos, the convergence of the Panama, Humboldt, and Cromwell currents bring nutrient-rich waters to the northern, central, and western parts of the archipelago (Seminoff et al. 2008, Lezama-Ochoa et al. 2019a. The distinct oceanic currents in the Galapagos and around mainland Ecuador result in high productivity in both areas (Lezama-Ochoa et al. 2019a Anderson et al. 2011, Graham et al. 2012, Lezama-Ochoa et al. 2019a, Putra et al. 2020, Farmer et al. 2022, supporting the hypothesis that given favorable oceanographic conditions, giant manta rays prefer to stay rather than travel long distances. We suggest that the genetic differentiation found among seasons (2013−2018 in relation to 2014−2017) in mainland Ecuador could be partially explained by fluctuations in temperatures and prey availability during phenomena like the El Niño−Southern Oscillation (ENSO) climate pattern (Osgood et al. 2021).
Seasons 2013 and part of 2018 were characterized by a strong coastal La Niña event, in contrast to the El Niño event reported for 2014−2016 (IMARPE 2021). La Niña is an oceanic and atmospheric phenomenon that is the colder counterpart of El Niño. During El Niño, winds that blow warm water from the Ecuadorian Pacific towards the Asian Pacific die down; at the same time, the influx of the Antarctic Humboldt Current, which brings cold water to the equator, is reduced (Montecino & Lange 2009, Edgar et al. 2010, Grados et al. 2018. The combination of these 2 events prevents the upwelling of cold and nutrient-rich deep ocean water towards the surface. This, in turn, starves primary levels of the food chain, hindering ocean production in general (Echevin et al. 2011, Grados et al. 2018). This process is particularly im portant for mobulid rays, as it has been suggested that gradual shifts in temperature could influence distributional changes since mobulids tend to track spatial changes in prey availability and are prone to migrate to areas with high primary productivity ( It is important to highlight the need for further studies to better understand the influence of ENSO on the genetic composition of M. birostris.

CONCLUSIONS
Our study is a significant step toward a better understanding of the genetic diversity and population structure of Mobula birostris in important aggregation sites in mainland Ecuador and Galapagos. Results suggest a low migration rate between these 2 locations, which could be explained by a lack of long-distance (>1000 km) movement by M. birostris individuals given the favorable oceanographic conditions at each site. In particular, the permanent upwelling systems associated with the presence of Humboldt, Pa nama, and Cromwell currents in these 2 areas may influence the spatial and genetic composition of M. biros tris, a filter-feeding organism that is known to chase ephemeral boosts of productivity (Hearn et al. 2014).
The well-defined genetic structure found in giant manta rays between Galapagos and mainland Ecuador suggests that individuals from these 2 locations form separate populations, and the limited gene flow between the 2 aggregation sites does not support the idea of a panmictic population in Ecuadorian waters. This information is important to delimit appropriate management units that should be monitored and managed separately in M. birostris conservation (IUCN 2020). Moreover, this information can help guide policies to control and eliminate illegal fisheries by implementing self-enforcement and including local communities in small-scale management actions, as has been previously successful in the Raja Ampat Shark and Ray Sanctuary in Indonesia (Stewart et al. 2016). These efforts will help prevent de clines in populations of Endangered species such as the giant manta ray in Ecuadorian waters.