Population recovery changes population composition at a major southern Caribbean juvenile developmental habitat for the green turtle, Chelonia mydas

Understanding the population composition and dynamics of migratory megafauna at key developmental habitats is critical for conservation and management. The present study investigated whether differential recovery of Caribbean green turtle (Chelonia mydas) rookeries influenced population composition at a major juvenile feeding ground in the southern Caribbean (Lac Bay, Bonaire, Caribbean Netherlands) using genetic and demographic analyses. Genetic divergence indicated a strong temporal shift in population composition between 2006–2007 and 2015–2016 (ϕST = 0.101, P < 0.001). Juvenile recruitment (<75.0 cm straight carapace length; SCL) from the north-western Caribbean increased from 12% to 38% while recruitment from the eastern Caribbean region decreased from 46% to 20% between 2006–2007 and 2015–2016. Furthermore, the product of the population growth rate and adult female abundance was a significant predictor for population composition in 2015–2016. Our results may reflect early warning signals of declining reproductive output at eastern Caribbean rookeries, potential displacement effects of smaller rookeries by larger rookeries, and advocate for genetic monitoring as a useful method for monitoring trends in juvenile megafauna. Furthermore, these findings underline the need for adequate conservation of juvenile developmental habitats and a deeper understanding of the interactions between megafaunal population dynamics in different habitats.

turtles 2 . Mixed aggregations in key developmental habitats present unique challenges to conservation and management because anthropogenic impacts potentially affect the future generations of multiple populations 2 . The conservation risk of population overlap is exacerbated for vulnerable populations, such as small populations of threatened species. Understanding the dynamics and composition at key developmental habitats, such as juvenile feeding grounds, is thus of paramount concern to conservation and management 2 .
However, little is known about the temporal dynamics at key developmental habitats even though recruitment is unlikely to remain constant [3][4][5] . Dispersal of marine species is governed by different processes at different life stages. Ocean currents influence dispersal of early life stages and temporal fluctuations in ocean currents can produce changes in the population composition at juvenile feeding grounds 3,4 . In addition, juvenile recruitment to feeding grounds can fluctuate over time due to environmental stochasticity affecting the production and survival of offspring 4 as well as long-term declines in reproductive success 5 . Human activities too can cause negative and positive changes, both through destructive activities as well as conservation measures. In the extreme case that human activities result in the local extinction of breeding populations, these extirpated populations will no longer produce recruits. By contrast, increases in reproductive success, for example due to conservation measures, can enhance recruitment 6 and possibly change the dynamics at juvenile feeding grounds.
A prime example of where the dynamics and composition at juvenile feeding grounds could be altered by conservation measures are green turtles (Chelonia mydas). Past human exploitation and habitat degradation has decimated sea turtle populations globally [7][8][9] . In the Caribbean, the present abundance of green turtles is at less than one percent of pre-exploitation levels based upon historical data 9 . Recent studies have highlighted that some rookeries are showing signs of population recovery following past sea turtle conservation measures [10][11][12][13][14] . Protection of nesting females has led to substantial increases in nesting trends, i.e. the number of nests produced, in several rookeries 10,12,13 . Increasing nesting trends as high as 14% per year have been reported in rookeries in Florida 10,14 . Rookeries in the eastern Caribbean seemed to have recovered at a slower pace compared to the north-western Caribbean 13,14 . Adult female abundance at the Aves Island, Venezuela, rookery has been increasing by approximately 5% per year during recent decades 13 . Little is known, however, what effects this population recovery may have on the dynamics and composition at juvenile feeding grounds.
In sea turtles, hatchlings disperse from natal rookeries to oceanic developmental habitats through a combination of active swimming and passive drifting in ocean currents 15,16 . Juvenile sea turtles later recruit to coastal feeding grounds 17 shared by multiple rookeries. Juvenile dispersal is influenced by factors, such as; ocean currents 18 , distance among rookeries and feeding grounds 18 , natal homing 19 and adult female abundance at rookeries 20 . Given the influence of abundance 20 and reproductive output at rookeries 4,5 upon sea turtle dispersal, it is possible that the population composition at juvenile feeding grounds changes as a result of differences in rookery recovery rates.
The hypothesis that the population composition at juvenile feeding grounds changes due to differential recovery rates can be tested by studying temporal genetic heterogeneity at feeding grounds using mitochondrial DNA (mtDNA). There is strong population structure in mtDNA diversity among sea turtle rookeries as a result of natal homing in adult females 21 . MtDNA markers can therefore be used as genetic tags to estimate juvenile recruitment into a feeding ground by assessing the relative proportion different rookeries contribute to a juvenile feeding ground (e.g. 18,20 ). If recovery rates differ among genetically diverged rookeries, temporal genetic heterogeneity is expected to increase with time, i.e. reflecting a directional change in mtDNA diversity. However, investigating changes in recruitment over longer timeframes is warranted because juveniles can spend up to 15 years at feeding grounds before moving to deeper feeding grounds (e.g. 22,23 ), which could lead to autocorrelation in mtDNA diversity between successive years. Autocorrelation in mtDNA diversity can be further diminished by investigating changes in recruitment using only recent recruits, i.e. small juveniles that presumably arrived recently at a feeding ground.
The present study assessed the relationship between differential recovery rates and changes in the population composition at a major juvenile green turtle feeding ground 24 located in Lac Bay (Bonaire, Caribbean Netherlands) in the southern Caribbean. Decadal changes in population composition were investigated by estimating temporal genetic heterogeneity and assessing changes in juvenile recruitment between 2006-2007 and 2015-2016. Decadal changes in juvenile recruitment were correlated with rookery recovery rates that were estimated as the product of population growth rates and adult female abundance to account for variation in population size among rookeries. The genetic and demographic analyses were conducted at both the level of individual rookeries and rookeries grouped into three regions reflecting mtDNA stocks recognized in sea turtle conservation and management 25 .

Results
Genetic diversity and temporal genetic structure. The sequence of 474 base pairs of the mtDNA control region was determined in 332 juvenile green turtles (30.0-75.0 cm maximum straight carapace length (SCL); mean ± SD = 50.66 ± 9.71 cm) sampled in Lac Bay between 2006 and 2016 (Tables S1-S2). Re-sequencing revealed one discrepant mtDNA sequence, corresponding to a sequencing consistency rate at >99%. Nineteen mtDNA haplotypes were detected among the 332 mtDNA sequences (Tables S1-S2), including a previously unreported mtDNA haplotype (GenBank accession number: MN481527; Fig. S1) designated 'CM-A76' in accordance with the commonly employed Atlantic green turtle mtDNA haplotype nomenclature, e.g. 26 . Adding 41 recaptured individuals increased the final sample size to 373 mtDNA control region sequences. The most common mtDNA haplotypes were CM-A03 (52%) and CM-A05 (27%) followed by CM-A01 (9%) and CM-A08 (3% ;  Table S1). Haplotypes CM-A01 and CM-A03 increased in frequency during the study period while the frequency of CM-A05 decreased ( Fig. 1 insert). The number of haplotypes and haplotype diversity varied among years without any apparent temporal trend (Table S2) Table S5) and time (P = 0.012, r 2 = 0.21). However, the elevated genetic divergence could be due to sampling variance and smaller sample sizes (Table S3). temporal changes in juvenile recruitment. Temporal changes in juvenile recruitment were observed for the eastern (ΔC = −0.22) and north-western Caribbean (ΔC = 0.26; Fig. 2A; Table 1). The contribution from the north-western Caribbean region to the Lac Bay juvenile feeding ground increased from 12% in 2006-2007 to 38% in 2015-2016. By contrast, the contribution from the eastern Caribbean decreased from 40% (2006)(2007) to 18% (2015)(2016). Juvenile recruitment from the south-western Caribbean region appeared constant throughout the study period (ΔC = 0.02; Fig. 2A; Table 1). Overall recruitment from the southern Atlantic was low  Table 1). The mean contribution of the Cayman Islands was higher than expected (2006-2007: 7%; 2015-2016: 1%) given the low abundance of adult females at this rookery 27 , but this is partially explained by the unweighted prior used in the mixed stock analysis.
Including only <50.0 cm SCL juveniles in the mixed stock analysis resulted in stronger increased recruitment from the north-western Caribbean (10% to 50%; ΔC = 0.41) and decreased recruitment from the eastern Caribbean (52% to 13%; ΔC = −0.  Table 1). These trends were also reflected in contributions estimated at the level of rookeries (Fig. 2D). Overall recruitment from the south-western Caribbean was stable but lower for <50.0 cm juveniles (ΔC = 0.02; 31% in 2006-2007 and 33% in 2015-2016; Fig. 2D). However, there was also greater uncertainty in regional-and rookery contribution estimates. correlation with rookery recovery trends. The highest population growth rate (r = 0.136) and reproductive output (Nr = 3,944) was estimated for the north-western Caribbean region (Table 1). Lower population growth rates were estimated for the eastern Caribbean (r = 0.066), south-western Caribbean (r = 0.017) and southern Atlantic (r = 0.035) regions. Nesting trends in Suriname were not linear but seemed to increase mainly during the last decade (Table S6). The expected reproductive output of the south-western Caribbean (Nr = 2,240) was higher than that estimated for the eastern Caribbean (Nr = 1,448) due to the high adult female abundance at the Tortuguero rookery in Costa Rica ( Table 1). The scarcity of rookery-specific data necessitated assuming regional population growth rate estimates for the Buck Island, Guadeloupe, Fernando de Noronha, Trindade Island, Bioko Island and São Tomé and Príncipe rookeries.
Spearman's correlation coefficient ρ was estimated at 1.00 (P = 0.083) at the regional level for both <75.0 cm SCL and <50.0 cm SCL juveniles. The monotonic association between ΔC and Nr, i.e. the ranking of ΔC and Nr matched perfectly at the regional level but was nearly significant, which was probably the consequence of only four data points that resulted from grouping rookeries into regions. At the level of rookeries, the monotonic association between ΔC and Nr was non-significant for both <75.0 cm juveniles (ρ = 0.46, P = 0.16) and <50.0 cm juveniles (ρ = 0.13, P = 0.71). However, logistic regression suggested Nr was a nearly significant predictor for whether recruitment increased or decreased for <75.0 cm SCL juveniles (P = 0.09), but not for <50.0 cm SCL juveniles (P = 0.25).
A linear model incorporating the natural logarithm of adult female abundance N (P = 0.047) and geographical distance D (P = 0.094) was found to best describe mean rookery contribution estimates in 2006-2007 www.nature.com/scientificreports www.nature.com/scientificreports/ for <75.0 cm SCL juveniles (Table S7) (Table S7).

Discussion
The present study assessed the changes in population composition at a major juvenile green turtle feeding ground located in Lac Bay, Bonaire, during the last decade. Genetic and demographic analyses suggested an increase in the proportion of juvenile green turtles in Lac Bay from rapidly recovering rookeries in the north-western Caribbean during the study period. In the north-western Caribbean, recovery of rookeries has previously been associated with increases in juvenile abundance at local, proximate feeding grounds 6 . However, juvenile green turtle abundance did not increase significantly in Lac Bay during the last decade 28 . If abundance were stable, the observed temporal changes in juvenile recruitment to the Lac Bay feeding ground would have reflected changes in the abundance of juveniles contributed by different Caribbean rookeries. If this presumption is correct, past sea turtle conservation measures in the north-western Caribbean could have resulted in increased juvenile abundance even at distantly located feeding grounds. By contrast, fewer juveniles originated from the eastern Caribbean and southern Atlantic. The decreased juvenile recruitment from the eastern Caribbean could be an early warning signal that reproductive output is declining in that region 29 . In the northern Great Barrier Reef, increased tidal inundation and rainwater flooding have been associated with reduced hatching success at the Raine Island green turtle rookery and decreased juvenile recruitment to local feeding grounds 5 . Rookeries in the southern Great Barrier Reef have been recovering during the last few decades 10 also suggested a potential effect of population recovery on juvenile recruitment 5 , though a lack of nesting trend data from rookeries outside the Great Barrier Reef precluded investigating this hypothesis 5 .
However, differences in nesting trends at rookeries and recruitment to juvenile feeding grounds are difficult to interpret 29 . Juvenile sea turtles form mixed aggregations at feeding grounds comprised of individuals from many different rookeries 18,20 that can differ in nesting trends 10,13,14 . Furthermore, adult females within a rookery will often utilize different geographically disparate feeding grounds 30 . Differential feeding habitat use can contribute to skewed contributions to reproductive success 31 . Hatchling mortality and nest production vary considerably between successive years and can lead to interannual fluctuations in juvenile recruitment to feeding grounds 4,32 . Recruitment to juvenile feeding grounds may also depend on population densities relative to carrying capacities 29 . Dynamics at juvenile feeding grounds are also expected to lag behind changes in nesting trends due to the time between hatching and recruiting to coastal feeding grounds 29 . However, more long-term studies will be www.nature.com/scientificreports www.nature.com/scientificreports/ required to further understand the time lag between hatching and recruitment, and its effect on metapopulation dynamics.
The decreased recruitment from the eastern Caribbean can possibly be explained by a difference in the timing of nesting trends. Nesting trends were stable in Suriname, the largest green turtle rookery in the eastern Caribbean, between the 1970's and the 2000's and started to increase during the last decade, while nesting trends in the Archie Carr National Wildlife Refuge green turtle rookery in Florida increased since the 1990's 10 . Green turtles spend approximately 3 to 5 years in oceanic habitats 17 before recruiting to coastal feeding grounds at a SCL of 25 to 35 cm 23 . Juveniles arrive in Lac Bay at 35 to 40 cm (SCL) at an age of 6 to 10 years depending on their rate of growth 6,33 . Changes in juvenile recruitment to Lac Bay are therefore expected to lag up to a decade behind changes in nesting trends. Given that nesting trends in Suriname have started to increase during the last decade, it is possible that juvenile recruitment from the eastern Caribbean will increase in the near-future. This does not explain why recruitment from the eastern Caribbean decreased despite long-term stability in nesting trends. However, nesting trends may be a poor proxy for reproductive output. For example, hatching success may be very sensitive to environmental fluctuations, such as changes in sand temperature at nesting beaches 34 . Nesting trends may therefore appear stable, but nonetheless result in few recruits due to low hatching success.
Sea turtles may shift to other feeding grounds between different size classes as a result of juvenile natal homing 19,35,36 . This was recently demonstrated in a study of juvenile green turtles in Japanese feeding grounds, where the contribution of local rookeries was higher for larger (i.e. 50-70 cm SCL) juveniles 35 . In the present study, recruitment differed between the two size partitions (i.e. <50.0 cm and <75.0 cm) within 2006-2007 and 2015-2016. Recruitment from the south-western Caribbean was higher overall when all juveniles were analysed, while recruitment from the north-western and eastern Caribbean was higher for small juveniles. These differences are potentially due to juvenile natal homing, reflecting the emigration of larger juveniles of north-western or eastern Caribbean origins from Lac. In addition, this implies that studying the recruitment dynamics of small juveniles is warranted to understand the link between population dynamics at feeding grounds and rookeries if larger juveniles shift between feeding grounds.
The stable abundance of green turtles 28 raises the possibility that the Lac Bay feeding ground is at or near carrying capacity. Species can respond to increased densities by dispersing to other feeding grounds 37 and it has been argued that dispersal in green turtles may in part be density-dependent 33 . Increased abundances at feeding grounds may lead to an increased propensity to disperse as local abundance nears carrying capacity and competition for resources increases 6,33 . For example, the exceeding of local carrying capacity in a juvenile feeding ground in the Bahamas resulted in net emigration of individuals over subsequent years until abundance stabilized 33 . If a feeding ground is at carrying capacity and dispersal propensity is equal for all individuals within a feeding ground, the equilibrium population composition is expected to be determined by the relative levels of recruitment from various source populations. Population composition may therefore change as a result of recruitment from one rookery outweighing recruitment from another rookery. In other words, rookeries contributing a larger number of offspring to a feeding ground may increasingly dominate shared feeding grounds over time, akin to a 'displacement effect' . This displacement effect could in part explain decreased recruitment despite stable nesting trends in the eastern Caribbean. Further monitoring of juvenile recruitment and abundance at feeding grounds in relation to population dynamics and reproductive output at rookeries is warranted to investigate whether such a displacement effect occurs.
Short-term fluctuations in mtDNA diversity between successive years in sea turtles can be a result of stochasticity in reproductive output 4 . By contrast, if long-term effects such as declines 5 or increases in hatching success, lead to differential juvenile recruitment to feeding grounds over time, genetic signals of 'directional change' in mtDNA diversity over multiple years, e.g. a decade, are expected. Short-term and long-term effects can therefore possibly be disentangled by assessing the direction and nature of changes in genetic heterogeneity during longer time intervals, as was done in the present study. The degree of temporal genetic heterogeneity observed during this study correlated strongly with time and suggested a directional change in mtDNA diversity. The temporal genetic heterogeneity we observed in Lac Bay was similar to levels recorded among different green turtle feeding grounds (e.g. φ ST = 0.168 between Barbados and Almofala in northern Brazil 38 ). These findings suggest the observed temporal genetic heterogeneity is in part explained by long-term changes in juvenile recruitment though the presence of short-term effects could not be rejected. However, short-term and long-term effects are not necessarily mutually exclusive 4 .
Genetic changes within rookeries, e.g. via genetic drift 39 , during the study period could, in theory, lead to a false-positive signal of juvenile recruitment. However, this is unlikely to explain the findings in the present study since no temporal heterogeneity in mtDNA has been observed at any rookeries so far [40][41][42] , although these studies only lasted two to three years. Temporal heterogeneity in mtDNA diversity has been reported at a loggerhead sea turtle rookery in Florida but could simply be due to sampling variance 43 . Second, possible changes in mtDNA diversity at the regional level during the study would require gene flow among distant rookeries which is unlikely given the high degree of natal homing to rookeries observed in nesting female sea turtles 21 .
Genetic assignment methods such as mixed stock analysis have reduced statistical power to determine the origin of individuals when putative sources are genetically similar 44 . Green turtle rookeries in the Caribbean differ substantially in mtDNA haplotype composition at the regional level but less so within regions. Accordingly, the estimated contribution from individual Caribbean rookeries was subject to a high degree of uncertainty in this study and complicated making inferences at the level of individual rookeries. Sequencing longer fragments of the mtDNA control region 45 , mitochondrial short tandem repeats 46 or mitochondrial genomes 47 could possibly reduce uncertainties in mixed stock analysis. However, this would require generating novel genetic data for all major Caribbean green turtle rookeries, which was beyond the scope of this study.
Effective conservation measures are key in an era characterized by an accelerated loss of biodiversity driven by anthropogenic activities. Sea turtle conservation and management has typically focused upon protecting nesting (2019) 9:14392 | https://doi.org/10.1038/s41598-019-50753-5 www.nature.com/scientificreports www.nature.com/scientificreports/ beaches 48 , even though nesting beaches constitute only a small, albeit vital, part of the life history of sea turtles. The increasing numbers of juveniles originating from rookeries that are showing signs of recovery observed in the present study are encouraging and highlights the success of current attention to sea turtle nesting beaches. In addition, the present study demonstrates genetic monitoring may represent a useful method for monitoring trends in juveniles, which may provide an early warning signal for declining reproductive success 29 and improve our understanding of sea turtle metapopulation dynamics 48 . However, sufficient feeding habitat quality and quantity is required for juveniles to mature and contribute to future generations 48 , and fortunately feeding habitats are increasingly receiving attention 49 . Climate change, invasive species and habitat degradation continue to threaten developmental habitats of sea turtles [50][51][52] . In Lac Bay, the invasive seagrass Halophila stipulacea has expanded rapidly during recent years 53 and a recent study demonstrated green turtles facilitate the expansion of H. stipulacea through selectively grazing on the native seagrass Thalassia testudinum 28 . If developmental habitat quality and quantity are not maintained, increased abundances of juveniles may elevate intra-specific competition for resources and increase risks of overconsumption, e.g. overgrazing of seagrass meadows 54 , and habitat collapse 55 . An understanding of the interactions between nesting trends at rookeries, recruitment and dynamics at feeding grounds as well as an understanding of the ecological interactions within feeding grounds is required to ensure adequate protection of both adult breeding and juvenile developmental habitats in endangered marine megafauna. A 474 base pair fragment of the mitochondrial DNA control region was amplified by nested PCR (polymerase chain reaction) amplifications 57,58 . An initial PCR amplification was conducted using the primers CM15412F (forward; 5′-AAAGCATTGGTCTTGTAAACC-3′) and CM16333R (reverse; 5′-TATGTCAGTTTGGTCAGTCTC-3′) followed by a PCR amplification using the primers CM15791F (forward; 5′-CAACCATGAATATTGTCACAGT-3′) and CM15984R (reverse; 5′-CATTCAACCAAAGGCCTTTTA-3′). PCR amplifications were conducted in a 10 μL reaction volume containing 1 μM of each primer, 1X standard  59 . Cycle-sequencing was performed with primers CM15412F (for the first fragment) and CM15791F (for the second fragment) using an ABI BigDye ® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems Inc.) following the manufacturer's instructions but using 1/16 of the BigDye ® Ready Reaction Mix.

Methods
The PCR conditions consisted of 25 cycles each with 10 seconds at 96 °C, 5 seconds at 50 °C and 4 minutes at 72 °C. Excess nucleotides and primers were removed by ethanol/EDTA precipitation. The cycle-sequencing products were re-suspended overnight in deionized formamide and the order of cycle-sequencing products was resolved by capillary electrophoresis on an Applied Biosystems 3730xl DNA Analyzer ™ (Life Technologies Inc.) at the University of Groningen.
MtDNA control region consensus sequences were assembled from one or more DNA sequences per DNA extraction using a custom software pipeline (Palsbøll, unpublished). Individual mtDNA sequences from each sample were assembled, the consensus sequence estimated from the Phred quality scores and subsequently aligned to a mtDNA control region green turtle reference sequence (GenBank accession number: JN632497.1) and truncated to 474 bp using MIRA ver. 4.9.5.2 60 . Assembled sequences were aligned with MUSCLE ver. 3.8.31 61 in SEQOTRON ver. 1.0.1 62 using the default parameter settings. Aligned sequences were visually inspected and manually corrected. A total of 96 randomly chosen DNA extractions were re-sequenced to assess sequencing consistency rates.
The data generated in the present study were combined with mtDNA sequence haplotype frequency data from juvenile green turtles (SCL < 75. MtDNA haplotype diversity 64 and nucleotide diversity 64  www.nature.com/scientificreports www.nature.com/scientificreports/ Population genetic analyses were repeated using only recent recruits (i.e. small juveniles with SCL-max <50.0 cm).
Temporal changes in juvenile recruitment were inferred estimating the contribution of Caribbean rookeries (Fig. 1) to the Lac Bay turtle feeding ground in 2006-2007 and 2015-2016 using the Bayesian mixed stock analysis approach implemented in BAYES ver. 11/23/11 69 . Published mtDNA control region haplotype (also 474 bp) data from Atlantic green turtle rookeries 40,47,[70][71][72][73][74][75][76][77] were used as source populations (Appendix 1). Two mixed stock assessments were conducted: (1) rookeries were grouped into regions (Fig. 1); 'north-western Caribbean' , 'south-western Caribbean' , 'eastern Caribbean' 20 and 'southern Atlantic' and (2) individual rookeries representing source populations. Additional mixed stock assessments were conducted using only small juveniles (<50.0 cm SCL) that presumably represent recent recruits and are more informative of recent changes in recruitment dynamics. A uniform prior was used where the prior contribution was equal among each source, i.e. prior contribution was set to 1/k where k was the number of putative sources, to avoid bias in subsequent demographic analyses. The employed MCMC settings are listed in Table S8. A Gelman-Rubin shrink factor below 1.2 was inferred as MCMC chain convergence 69 . Temporal changes in juvenile recruitment were estimated as: where C i was the mean contribution estimated for 2006-2007 and C j was the mean contribution estimated for 2015-2016.
Annual population growth rates at north-western, south-western Caribbean rookeries were estimated by calculating the weighted mean of annual population growth rates at nesting sites (data from Mazaris et al. 14 ) within rookeries (Tables S9-S10). Nest site-specific annual population growth rates were weighted using estimates of adult female abundance at nesting sites 27,73 . For the eastern Caribbean, annual population growth rates at the French Guiana and Suriname rookeries were estimated using temporal trends in the number of recorded nests (Table S6) using the approach adopted from Mazaris et al. 14 where r is the annual population growth rate, N L is the mean number of recorded nests in the last three years of the time series, N F the mean number of recorded nests in the first three years of the time series and n the length of the time series. The estimate of the annual population growth rate at the Aves Island rookery was obtained from García-Cruz et al. 13 . Southern Atlantic rookery annual population growth rates were obtained from Mazaris et al. 14 (Bioko Island and Rocas Atoll) and Weber et al. 11 (Ascension Island). Regional annual population growth rates were estimated as weighted mean annual population growth rates for rookeries within regions.
The product of the annual population growth rate rand the adult female abundance N (Nr) was used as a proxy for reproductive output in order to account for variation in adult female abundance among rookeries. The association between ΔC and Nr, was tested using nonparametric Spearman rank-order correlation in R ver. 3.5.3 78 . In addition, we tested how well Nr predicted whether recruitment increased or decreased over time at the level of rookeries (i.e. as a binary response variable) using logistic regression in R ver. 3.5.3. Multiple linear regression was used to test the effect of N, geographic distance (D; estimated as the shortest distance between source populations and the study site using the geosphere R package ver. 1.5-10) and Nr on mean rookery contribution estimates during 2006-2007 and 2015-2016 for both size classes. Model selection was performed via a stepwise algorithm using the Akaike Information Criterion (AIC) in R ver. 3.5.3.

Data Availability
The DNA sequence of the novel mtDNA haplotype CM-A76 has been deposited in GenBank (accession number: MN481527).