A total of 1,111 turtle blood or tissue samples were collected from 20 Chelonia mydas rookeries located across the East, Central, and West Pacific ( Figure 1 , Table 1 , Supplementary Table 1 ). Here, a rookery (or deme) generally includes one or more beach nesting sites that are in close geographic proximity (<50 km) and are likely to be used by females from the same breeding population with a shared demographic history. Tagging studies show that females often nest on multiple nearby beaches (Alvarado and Figueroa, 1991; Dethmers et al., 2006; Fonseca et al., 2018) and this supported combining some adjacent nesting beaches into rookery units (Dutton et al., 2014b; Jensen et al., 2019) ( Supplementary Table 1 ).

Samples, including some previously used in FitzSimmons et al. (1997a; 1997b); Dethmers et al. (2006); Roden et al. (2013), and Dutton et al. (2014a; 2014b) were collected between 1989 and 2017 from nesting females and dead hatchlings. Samples were preserved in 70% ethanol, a NaCl saturated 20% DMSO solution, or saturated NaCl solution as outlined in Dutton (1996). All samples were archived, stored (-20°C), processed, and analyzed at the Southwest Fisheries Science Center in La Jolla, California. Genomic DNA was extracted from tissues using multiple methods appropriate to tissue type, storage method, and condition (FitzSimmons et al., 1997b; Roden et al., 2013; Dutton et al., 2014b).

All samples in the study were genotyped using polymerase chain reaction (PCR) with ten microsatellite loci characterized in Dutton and Frey (2009); Frey et al. (2013), and Roden et al. (2013). Loci included A6, B108, B116, C102, D1, D108, D115, D2, D102, and D107 and each forward primer was labeled with one fluorescent tag, HEX or FAM. Microsatellite loci were amplified using approximately 2-15 ng template DNA and reaction conditions outlined in Dutton and Frey (2009) with the following adjustments: 0.1 mM each dNTP, 0.2 mg/ml bovine serum albumin (BSA), and 0.06 U/μl Taq DNA polymerase.

All PCR products, including positive and negative controls, were checked for amplification using 2% agarose gels with ethidium bromide staining (Sambrook et al., 1989). Amplified products were diluted and pooled in pairs of two loci per capillary (B116/D2, D115/C102, D107/D1, D108/D102, A6/B108). Pairs were chosen to avoid overlapping size ranges and genotyped with an ABI Prism 3730 DNA Analyzer using a ROX 500 fluorescent size standard (ABI, USA) and accompanying GeneMapper software v4.0 (Applied Biosystems) (see Supplementary Table 2 for genotype data). To facilitate consistency and normalization of allele sizes, standard reference samples with known genotypes were included within each PCR reaction plate to serve as positive controls. Approximately 20% of all samples were PCR amplified multiple times to replicate genotypes and assess the genotyping error rate across the microsatellite loci.

Genotypes were generated for 1,111 individuals (representing 20 rookeries) across 10 microsatellite loci. After QC tests to remove potential duplicates and samples with greater than 30% missing genotype data, 1,064 samples were analyzed. No signs of linkage disequilibrium were detected in the data. The overall genetic diversity per locus ranged from 3-24 alleles, with 1.60 to 8.78 effective alleles ( Table 2 ) with an average genotyping error rate of 0.02 across all loci. The observed heterozygosity among populations ranged from 0.658 to 0.795 ( Table 1 ). Four loci (D108, D1, B116, C102) showed signs of heterozygote deficiency and consequently departed from HWE ( Table 2 ). D108 deviated significantly from HWE in 16 of 20 rookeries, D1 in 7, B116 in 7, and C102 in 4 out of 20 rookeries (p < 0.05; see Supplementary Table 3 ). Significant heterozygote deficits were estimated across loci in 14 out of the 20 rookeries ( Table 1 ). Microchecker results indicated that null alleles may be present in 15 of the rookeries for D108, and to a lesser extent for D1 (seven), B116 (four) and C102 (four) and null.all() produced similar results, with D108 the primary locus with potential null alleles. There was no evidence of large allelic dropout at any loci across all rookeries. G _IS values in the four affected loci were all positive ( Table 2 ) indicating heterozygote deficiencies, but G _IS was not correlated with G _ST making HWE departures due to Wahlund effects unlikely (Waples, 2015). Null alleles were not likely to be driving heterozygote deficits, since the number of homozygotes was not correlated with missing data in any locus (Pearson’s r = -0.08 – 0.39, p = 0.08 – 0.86) ( Wang et al., 2012; Dharmarajan et al., 2013; Dąbrowski et al., 2015) and since G_IS was not correlated with G_ST (De Meeûs, 2018). In order to further evaluate whether artifacts introduced by null alleles, especially in D108, might affect the results of standard population structure analysis, we conducted pairwise F_ST tests with and without the potentially problematic loci, as well as D108 allele frequencies adjusted for null alleles (see below).

F’_ST and G’’_ST based on the full 10 loci yielded results that were similar to standard F_ST, but Jost’s D estimator sometimes indicated a lack of differentiation between rookeries where that would be unlikely, such as between Lanyu, Taiwan and Pacific Costa Rica (see Supplementary Table 4 ), which was inconsistent with other results. All but eight of 190 pairwise rookery comparisons indicated significant differentiation based on pairwise F_ST tests (p < 0.05; Figure 2 , see Supplementary Table 4 ). Insignificant pairwise structuring was mostly found in the Southwest Pacific along the continental shelf of eastern Australia, or involved rookeries with sample sizes < 20 ( Figure 2 ). Furthermore, population structure between the eastern and western Pacific was clear, as seen in the PCoA ordinations of pairwise F_ST values ( Figure 3 ).

Pairwise F_ST results based on nine loci (excluding the most potentially problematic locus, D108), yielded a similar outcome with regard to the pairwise comparisons that were not significantly different with 10 loci. However, with D108 removed, PCR and MICH as well as CS and sGBR were no longer significantly differentiated ( Supplementary Table 5 ). When all 10 loci were included with allele frequencies for D108 adjusted for null alleles, the only difference was significant differentiation between nGoC and sGoC with the adjusted D108 dataset ( Supplementary Table 6 ). Removal of all four potentially problematic loci (D108, B116, C102, and D1) resulted in further lack of significant differentiation between several populations in the West and Central Pacific. The broad geographic patterns of structure were still evident with the six loci dataset, as well as most of the finer scale patterns within the East Pacific ( Supplementary Figure 2 ). In sum, while some loci displayed deviations from expectations, extensive exploration and testing of the data supported retaining every locus for subsequent analysis in order to optimize the ability to detect any finer scale structure that might exist.

Results of the POWSIM analysis for all 10 loci showed that the power to detect weak structure (F_ST = 0.001) was poor with sample sizes <50, but better for F_ST < 0.0025 with n >10 (see Supplementary Figure 1A ). Power was robust for detecting weak structure (F_ST= 0.0025) with sample sizes >50 within the Central North and East Pacific (see Supplementary Figure 1B ), and with sample sizes >20 in the Western and South Pacific (see Supplementary Figure 1C ).

All DAPC results retained 100 principal components as predictors for discriminant analysis, and 98-99% of the conserved variance. The broad scale pattern for the DAPC analysis was consistent with pairwise population structure ( Figures 2 , 3 ), showing three distinct clusters for East, Central, and West Pacific rookeries ( Figure 4A ), and also reflected much the same spatial patterns within the East Pacific ( Figure 4C ), Northwest Pacific ( Figure 4D ), and Southcentral/Southwest Pacific ( Figure 4B ) as the pairwise results.

The most information-rich locus in the microsatellite panel based on the DAPC loading values was C102, which was also the most polymorphic locus, having 24 alleles ( Table 2 ). With the exception of the Southwest Pacific, C102 alleles exhibited the highest loading values in all of our DAPC analyses, underscoring the value of retaining this locus (Pearman et al., 2022). In the broad-scale analysis, a single C102 allele accounted for 14% of the genetic variance between the East and West Pacific (axis 1, Figure 4A ). Another C102 allele accounted for approximately 17% of the variance on axis 2 of the same plot, differentiating Hawaii from the rest of the Pacific ( Figure 4A ). The least important locus in terms of loading values appeared to be D102, which was also the least polymorphic, with only three alleles. Summing the variance contributions from all alleles for each locus revealed a positive and roughly linear relationship between loading values and the total number of alleles for all DAPC analyses (data not shown) suggesting that the informativeness of each locus was associated with allelic diversity.

EEMS analyses were performed on three subsets of the data (Northwest Pacific, Southwest Pacific, and East Pacific; Figures 5A–C , 6A–C , respectively), corresponding to the subsets used for DAPC, because the relationship between genetic dissimilarity and geographic distance was not uniform across regions. The East Pacific rookeries had the most IBD-like dissimilarity interpolations and the best migration matrix was a circular, one-dimensional stepping-stone model ( Figure 5C ). The coastal distance between MICH and PCR was a relatively good predictor of genetic dissimilarity between these two continental-shelf rookeries, which appear more similar to each other than to the oceanic islands of the East Pacific ( Figure 5C ). The most genetically diverse rookery in this region was GAL ( Figure 6C ). This diversity is the amount of genetic dissimilarity observed between individuals within rookeries, and not an expression of heterozygosity or the number of effective alleles, which did not identify this higher level of diversity ( Table 1 ).

In the West Pacific, geographic distance was a poor predictor of genetic dissimilarity between rookeries, above and below the equator ( Figures 5A, B ). In the Northwest Pacific, the migration matrix that produced the best results was a full island model. In the Southwest Pacific the best matrix was a nearly full island model, but the model fit was slightly better when AS and FP weren’t directly connected to the Australian continental-shelf rookeries. The Australian rookeries had the highest within-deme (rookery) genetic dissimilarity (e.g. diversity) in the entire Pacific ( Figure 6B ). Coincidentally, Australian rookeries (CS, nGBR, sGBR) also had the highest genetic dissimilarity between demes ( Figure 5B ), suggesting that genetic distance increased more than geographic distance among these rookeries.

Our microsatellite dataset for Pacific green turtles achieved an ocean-basin scale nDNA population analysis for an endangered marine species and showed meaningful genetic structure at regional and population-level spatial scales. The microsatellite genotypes used herein were produced from many of the same tissue and DNA extractions as the mtDNA analysis of Jensen et al. (2019), making the two studies highly comparable and complementary. On the whole, the nDNA corroborated population patterns recovered from mtDNA, with the strongest feature of both datasets being a deep genetic partition between the East and West Pacific. The genetically intermediate Hawaiian rookery was also more related to the East Pacific in the nDNA as observed with mtDNA (Dutton et al., 2014a; Jensen et al., 2019; Hamabata et al., 2020). Additionally, nDNA segregated most rookeries from each other ( Figure 2 ), demonstrating that it too captures a natal homing signal, similar to mitochondrial genetic structure. Examples include the Northwest and Central Pacific, where most rookeries indicated significant population differentiation. This finding dispels the notion that widespread male-mediated gene flow occurs and reinforces the existence of natal homing, at least on a regional level, by males. Our results from the EEMS analysis indicating elevated within-deme (rookery) diversity in the Australian rookeries compared to the East Pacific are not surprising considering past mtDNA studies that identified the Southwest Pacific as a diversity hotspot for mtDNA. Central American rookeries display a shallower demographic history, which reflects a relatively recent recolonization of the East Pacific from Hawaii (Dutton et al., 2014a). At the same time, Australian turtle populations have much deeper lineages given that the West Pacific has sustained green turtles for much longer than the East Pacific (Jensen et al., 2019).

While the East and Central North Pacific (Hawaii) rookeries were extensively sampled, coverage in the Central, South and West Pacific was less comprehensive and included some rookeries with relatively small sample sizes. The French Polynesia rookery, for example, had only 12 samples and this may not have been adequate for detecting the level of population structure between it and Palau, American Samoa, and New Caledonia. POWSIM results indicate that our current assay is poor at detecting weak differentiation with sample sizes < 20, and a lack of structure between some minimally sampled rookeries and others could be due to a lack of power. Considering that the insularity of most sea turtle breeding populations is consistently confirmed by genetic analysis, instances where pairwise genetic structure is lacking, and the sample size is small, results should be interpreted cautiously.

Nevertheless, results of the EEMS analysis when considered together with results of the pairwise tests, provide insights that allow a more nuanced interpretation of the population structure analysis results by themselves for the under-sampled Australian rookeries. The fixation indices and DAPC results showed the nGBR and sGBR as genetically similar, and the CS rookery as differentiated from the nGBR and sGBR ( Figures 2 , 4B ), and yet the EEMS results showed that all these three rookeries (CS, nGBR, sGBR) had the highest genetic dissimilarity between demes ( Figure 5B ), suggesting that genetic distance increased more than geographic distance among these rookeries, a scenario more consistent with discontinuity rather than homogeneity. Samples from the sGBR and CS were few (n = 20 and n = 18, respectively) and these small sample sizes might not adequately characterize the elevated genetic diversity that the EEMS analysis revealed for the Australian rookeries ( Figure 6B ). If so, then the genetic dissimilarity between these rookeries may more or less reflect the dissimilarity within them, making it difficult to delineate population boundaries. Analyses such as DAPC combine the within- and among-deme genetic variation into the same plot, and intentionally minimize within-deme variation to better represent population structure visually. EEMS analysis separates these two different variance components instead and presents them individually in continuous geographic space trying to interpolate, rather than discriminate, the boundary between defined populations. Therefore, when comparing DAPC and EEMS results it is important to consider the different ways these analyses handle the same genetic variation.

Heterozygote deficiencies in some individual loci, and in global tests including all loci in several of the sampled rookeries, are most likely features of green turtle populations that have been subjected to a number of disruptive processes, and not the result of genotyping errors. For instance, local extirpation and recolonization, driven partly by Pleistocene sea-level disturbances to nesting beaches (Dethmers et al., 2006), could have led to more genetic drift than expected under stable conditions and HWE (Dutton et al., 2014a; Horne, 2014). In addition, recent human activities that have negatively impacted green turtles (Seminoff et al., 2015) may have eroded genetic diversity even in rookeries that remain relatively large, or reduced populations to the point where they are too small, or too inbred, to obtain a truly random sample (Dharmarajan et al., 2013). Nevertheless, the idea that genetic drift is the strongest signal in our data is supported by the fact that pairwise F_ST values (and analogs) tended to be more geographically intuitive than pairwise Jost’s D, which estimates connectivity based on the sharing of alleles rather than fixation (Jost et al., 2018). Moreover, fixation indices are sensitive to all causes of genetic drift, such as founder effects, population bottlenecks, as well as infrequent migration (Whitlock, 2011). Lastly, the use of DAPC, a method that does not assume HWE, recovered the same general population patterns as the fixation indices.

Overall, the nDNA population patterns are in agreement with previous mtDNA studies, with a few noteworthy exceptions. The first was a signal of genetic differentiation found between CNMI and Guam that has not been previously detected using mtDNA control region sequences (Dutton et al., 2014b). However recent whole mitochondrial genome sequencing has provided greater resolution to assess green turtle population structure and shows evidence of highly significant genetic differentiation between CNMI and Guam (Frey et al., in prep) in agreement with this study. This signal of differentiation is surprising due to the proximity of Guam and CNMI and may be the result of genetic drift due to small population size or a founder effect similar to the findings in the Hawaiian Islands by Frey et al. (2013). Guam and CNMI are both very small nesting populations with less than 50 females estimated in each (Seminoff et al., 2015). It is also noteworthy that the foraging populations around CNMI and Guam consist of mostly juveniles, with adults tending to be observed around the nesting season (Martin et al., 2016; Summers et al., 2017; Gaos et al., 2020). The significant differentiation that we found between Guam and CNMI with nDNA data suggests that males (presumably homing to their natal beaches) are likely to be mating with nesting females that are reproductively receptive during the breeding season in the vicinity of the nesting beaches, rather than in mixed courtship areas, thus limiting the opportunity for male-mediated gene flow (Frey et al., in prep).

Within the East Pacific, a greater degree of connectivity was found among Central American and continental Mexican rookeries in comparison to previous mtDNA population studies. With nDNA microsatellites, we found evidence for greater connectivity between the mainland rookeries of Costa Rica and Mexico than with the island rookeries of the Galapagos and Revillagigedo. This is interesting given that telemetry and tag recaptures show green turtles from Galapagos rookeries sharing foraging grounds in Central and South America with those from Mexico and Costa Rica (Seminoff et al., 2008; Hart et al., 2015). The genetic dissimilarity of Galapagos turtles could be partly explained by differences in migration routes between males and females that nest on the continent vs. those that nest on oceanic islands. If continental green turtles maintain coastal migratory behaviors (Blanco et al., 2012; Hart et al., 2015), then Mexican and Costa Rican turtles could be using overlapping migratory pathways along the coastal shelf when moving between nesting and foraging areas, and may opportunistically interbreed where migratory pathways intersect with courtship areas when the turtles become reproductively active in these areas just prior to the nesting season. This opportunistic mating is likely driven by the behavior of migrating males given that they are reproductively active for longer periods than the females, which tend to be reproductively receptive for shorter durations at the breeding sites (Ulrich and Parkes, 1978; Crowell Comuzzie and Owens, 1990; Limpus, 1993; Owens, 1997). Furthermore, green turtle mating activity generally decreases once nesting seasons begin but in Michoacán, mating persists into the nesting season (Alvarado and Figueroa, 1991; Delgado-Trejo and Alvarado-Figueroa, 2012). In contrast, Galapagos green turtles that forage in continental shelf waters would not migrate along the same coastal routes. Instead, they would head offshore toward the high seas, where they would be less likely to have chance mating encounters (Seminoff et al., 2008). Given the great distances involved (necessitating longer migration times), it is also unlikely that Galapagos turtles would be reproductively active or encounter reproductively receptive females before they depart continental foraging grounds to start breeding migrations. The same would be true for Islas Revillagigedo green turtles that forage throughout the Baja Peninsula and Southern California, USA (Juarez-Ceron et al., 2003; Dutton et al., 2019). Additionally, lenient natal homing has been found to be more prevalent in continental nesters returning to nesting sites along continuous coastlines than for nesters homing to isolated insular coastlines (Levasseur et al., 2019; Shamblin et al., 2020), where the former may share or have very similar features, such as magnetic fields used as a cue to identify the magnetic signature of a natal beach (Brothers and Lohmann, 2018).

Conversely, the lack of nDNA structure between the three Australian breeding regions of nGBR, sGBR and sGoC ( Figures 2 , 4B ) was somewhat unexpected because these three regions have nearly fixed mitochondrial control region differences and the phylogeographic affinity between the GoC and Indian Ocean rookeries evident in mtDNA lineages (Dethmers et al., 2006; Jensen et al., 2019). Previous work comparing nuclear and mtDNA saw significant genetic differentiation for both markers across all three nesting populations, except for nDNA between the nGBR and sGBR (FitzSimmons et al., 1997b), suggesting that male-mediated gene flow acts to homogenize these two populations. However, mating typically occurs near the nesting beaches in the months before the start of the nesting season (Wood and Wood, 1980; Limpus, 1993), and there is strong evidence that male turtles in Australia are as philopatric as the females, returning to breed in the vicinity of their natal beaches (FitzSimmons et al., 1997a) and thereby minimizing the opportunity for gene flow across populations.

Hence, the nuclear DNA connectivity between the nGBR and sGBR observed in both studies is best explained by rare mating events where breeding migration routes overlap with courtship areas (FitzSimmons et al., 1997a; FitzSimmons et al., 1997b). The sGBR male and female turtles occupy foraging areas across the entire Great Barrier Reef, Torres Strait (between Australia and Papua New Guinea) and even into the Gulf of Carpentaria (Limpus, 2008). The timing of opportunistic mating is not well understood but it is likely to occur during their breeding migration, when sGBR male and female turtles traverse major nGBR foraging and courtship areas, creating the potential for nuclear but not mitochondrial gene flow, most likely by migrating sGBR males in breeding condition traveling south and receptive nGBR females already at their courtship areas. It is important to note that courtship areas may be at least a few hundred kilometers away from nesting beaches, as observed in green turtles that nest in Tortuguero, Costa Rica, and mate along migratory routes off Panama, 240 km distant (Meylan et al., 2013). Similarly, the large nGBR courtship areas in Torres Strait are 200 km (or more) northwest from the main nesting beach at Raine Island (for nesting distribution see (DES; CMS, n.d.). At the Howick Group feeding ground (south of nGBR) which is used by both nGBR and sGBR turtles, two vitellogenic female green turtles caught during mating were satellite-tagged and immediately migrated ~340 km north to nest at Raine Island (DES; CMS, n.d.). This suggests a broad geographic range in courtship areas used by the nGBR rookery as also suggested for the sGBR rookery where courtship activity has been recorded at sites up to hundreds of kilometers distant from the main courtship areas (C.J. Limpus, unpublished data). The well-monitored foraging ground at the Howick Group (Jensen et al., 2018) could also be a site to further investigate questions of male mediated gene flow raised here.

Previously, Bountiful Island (sGoC) was shown to be distinct at four microsatellite loci and that has been attributed to differences in nesting phenology between the Gulf of Carpentaria (Austral winter) and the two GBR populations (Austral summer), which minimizes the chance of interbreeding (FitzSimmons et al., 1997b). This contrasts with our study, which finds genetic homogeneity for nuclear markers among all three rookeries. Furthermore, in this study, the nGoC and sGoC rookeries revealed different patterns of differentiation with the GBR rookeries. While the peak nesting activity differs between the Gulf of Carpentaria (April-October) and the Great Barrier Reef (October-March), low-level nesting in the Gulf of Carpentaria and the nGBR green turtle populations occurs year-round (Limpus, 2008). Tagging data (Limpus, 2008) show that some females nesting in the nGBR and sGBR take up residency in foraging areas across the Gulf of Carpentaria. Recent satellite telemetry research has tracked nGBR nesters tagged at Raine Island (nGBR) to foraging areas around Bountiful Island and other green turtle nesting areas in the GoC (DES; CMS). Hence, a similar interpopulation (opportunistic) mating scenario could be happening between the GoC and GBR populations during the months when there may be some overlap in the timing of migration. Reasons for discrepancies between the studies may include low statistical power in both nDNA studies due to low sample size, and it remains unresolved what level of admixture occurs between nesting populations in the GoC and the GBR. More research will be needed to fully resolve green turtle population dynamics along the Australian continental shelf.

In the Southwest Pacific, some green turtles from rookeries in New Caledonia reside in foraging grounds along the GBR, and vice versa, requiring them to cross the Coral Sea during breeding migrations (Limpus, 2008; Read et al., 2014). Nevertheless, there is no evidence for male-mediated gene flow between New Caledonian rookeries and those on the Australian shelf, suggesting that male turtles from New Caledonia do not breed in the GBR but return to New Caledonia to do so, and vice versa. Instead, results from this study indicate that chance mating encounters may be happening between green turtles from New Caledonian and Coral Sea rookeries while sharing foraging or migration corridors. More research is clearly needed to better understand the reproductive physiology and the mechanisms of male mediated gene flow between demographically independent nesting populations that share foraging habitats or migrate through courtship areas utilized by other genetic stocks. Sperm storage, female receptiveness, and understanding the timing and duration of vitellogenesis (females) and spermatogenesis (males) prior to the onset of their breeding migration all play a vital role in understanding these mechanisms (Kitayama et al., 2021). Additionally, the levels to which opportunistic mating might happen are highly influenced by the geography of rookeries, foraging areas, and migratory routes.

Our results prompt reconsideration of how population units for conservation have been defined for green turtles in some cases based solely on patterns of mtDNA variation. The incongruence we found in some cases between our nDNA results and mtDNA highlights the need for a more nuanced approach in determining the extent to which rookeries represent demographically independent populations for conservation (e.g., Guam and CNMI see above). Failure to recognize Guam and CNMI as separate MUs could misinform conservation expectations; both populations are severely depleted, and given the lack of demographic connectivity, it is unlikely that any recovery actions on one of the islands will influence population trends on the other (Seminoff et al., 2015). On the other hand, the increasing abundance trend of the large East Pacific nesting population in Michoacán is the result of sustained conservation measures enacted over the last 50 years to protect nests and breeding adults (Delgado-Trejo and Alvarado-Figueroa, 2012). The Michoacán population is clearly defined as a separate MU from Costa Rica based on mtDNA, however the nDNA connectivity, and implied male-mediated gene flow among these populations established in this study suggests that demographic benefits from conservation at one location extend beyond the MU boundaries. These findings provide a broader regional framework for future assessment of whole genomic variation and potential for local adaptation between the continental (Mexico and Central America) vs. island (Galapagos and Revillagigedo Islands) populations using whole genome sequencing approaches (Álvarez-Varas et al., 2021).

Furthermore, the regional genomic connectivity among mtDNA genetic stocks (MUs) indicated in this study has conservation importance with regard to the ability of certain rookeries to adapt to environmental change and can help inform metapopulation modeling approaches for sea turtles (Jensen et al., 2022). Because warmer incubation temperatures produce more female offspring, warmer temperatures due to climate change may produce a disproportionate abundance of females in green turtle populations (Chaloupka, 2004; Hawkes et al., 2007; Hays et al., 2010; Hays et al., 2017). Studies have projected a lack of male turtles in some populations in the near future. In the nGBR, immature turtles from one feeding ground (the Howick Group) were found to be 99% female (Jensen et al., 2018), and less than 1% of hatchings born at Raine Island during the 2018 and 2019 nesting seasons were estimated to be male (Booth et al., 2020). At present, the nGBR boasts one of the largest green turtle rookeries in the world (Seminoff et al., 2015), but it is very vulnerable due to climate change-related feminization. Some models suggest that climate change could exacerbate the problem and even eliminate some turtle colonies (Hawkes et al., 2007; Hays et al., 2017; Jensen et al., 2022), but do not account for scenarios involving male mediated gene flow. These models could be improved by accounting for male-mediated gene flow and a better understanding of its role in maintenance of genetic diversity and ultimately population persistence. The within-region differences discovered in the present study will need to be resolved with a more extensive population genomic assessment involving additional genome-wide variation, so that loci deviating from HWE, such as were seen in our study, can be evaluated more effectively in a genomic context where signal to noise ratios can be more finely tuned (Pearman et al., 2022).

In the meantime, however, the findings echoed with multiple analytical approaches used in this study warrant consideration because the effects of metapopulation processes are often difficult to infer from F_ST values alone (Pannell and Charlesworth, 2000). This lesson is more broadly applicable to other widely distributed marine species with complex, and often cryptic life history stages.

Finally, our study provides an extensive Pacific-wide nDNA baseline dataset for population assignment analyses (“stock ID”) for individual green turtles of unknown population origin, such as in fisheries bycatch, strandings, and foraging grounds. The ability to conduct such assignments has previously been limited to Mixed Stock Analysis approaches using mtDNA data that often result in large uncertainty and require adequate sample sizes (Jensen et al., 2013). Assignment of individuals with multi-locus markers is particularly useful for identifying the origin of fisheries bycatch more precisely (Stewart et al., 2016). While West Pacific rookeries are under-represented relative to those in the Central and East Pacific in this baseline dataset, it nevertheless includes regional representation that provides for ocean-scale stock assignment (Horne et al. 2023).

Rapid technological advancements make it increasingly feasible and cost effective to generate high quality whole genome sequence data that also characterizes genes under selection (Bentley et al., 2023), and thus overcome the limitations of using a relatively small number of microsatellite loci described in our study. Equally important however is the requirement to adequately sample populations. There continue to be the logistical challenges of obtaining adequate and representative samples of species such as sea turtles due to their threatened status, extensive distribution throughout the open ocean and numerous remote nesting regions of many countries. Through a coordinated, comprehensive sampling effort by an expansive network of collaborators, we present the results of a landmark study with an extensive set of samples; many previously used in published mtDNA studies. This international collaborative framework will need to be further expanded and maintained to address geographic gaps, and also temporal gaps in sampling, since some of these samples were collected over a 10-22 year timeframe, over 30 years ago, and may not adequately reflect any temporal shifts in diversity that may have occurred. There are relatively few studies that evaluate whether the geographic patterns of genetic structure vary significantly over ecological time scales in marine species (Toonen and Grosberg, 2011; Horne et al., 2012). Given their long lifespans and life history traits that include nest site fidelity and overlapping generations, it is assumed that spatial genetic patterns of population structure generally persist over decades in sea turtles, which is within the reproductive lifespan of individuals (Bjorndal et al., 2005; Formia et al., 2007; Velez-Zuazo et al., 2008; Joseph and Nishizawa, 2016). However, severely depleted sea turtle populations that have been subjected to high anthropogenic mortality that has likely reduced reproductive lifespans of adults, may be more prone to shifts in genetic diversity within shorter decadal timeframes. Vargas et al. (2022) found a difference in mtDNA haplotype frequencies in the extremely small Brazilian leatherback nesting population over a 28-year period that could represent contemporary genetic drift. Our study provides a baseline for monitoring temporal changes for green turtle rookeries into the future. Nuclear DNA (nDNA) markers are a valuable complement to gain a comprehensive understanding of population boundaries and male-mediated gene flow by incorporating the male genetic connectivity component to analysis. With highly migratory, long-lived species like marine turtles, extensive collaboration with sampling, genetic marker application, and analysis techniques are central to understanding the connectivity of populations within oceans and worldwide.