Genetic Structure and Natal Origins of Immature Hawksbill Turtles (Eretmochelys imbricata) in Brazilian Waters

Understanding the connections between sea turtle populations is fundamental for their effective conservation. Brazil hosts important hawksbill feeding areas, but few studies have focused on how they connect with nesting populations in the Atlantic. Here, we (1) characterized mitochondrial DNA control region haplotypes of immature hawksbills feeding along the coast of Brazil (five areas ranging from equatorial to temperate latitudes, 157 skin samples), (2) analyzed genetic structure among Atlantic hawksbill feeding populations, and (3) inferred natal origins of hawksbills in Brazilian waters using genetic, oceanographic, and population size information. We report ten haplotypes for the sampled Brazilian sites, most of which were previously observed at other Atlantic feeding grounds and rookeries. Genetic profiles of Brazilian feeding areas were significantly different from those in other regions (Caribbean and Africa), and a significant structure was observed between Brazilian feeding grounds grouped into areas influenced by the South Equatorial/North Brazil Current and those influenced by the Brazil Current. Our genetic analysis estimates that the studied Brazilian feeding aggregations are mostly composed of animals originating from the domestic rookeries Bahia and Pipa, but some contributions from African and Caribbean rookeries were also observed. Oceanographic data corroborated the local origins, but showed higher connection with West Africa and none with the Caribbean. High correlation was observed between origins estimated through genetics/rookery size and oceanographic/rookery size data, demonstrating that ocean currents and population sizes influence haplotype distribution of Brazil's hawksbill populations. The information presented here highlights the importance of national conservation strategies and international cooperation for the recovery of endangered hawksbill turtle populations.


Introduction
After hatching, sea turtles often present an epipelagic stage characterized by wide dispersal, many times over national boundaries and even oceans [1]. This stage is generally followed by recruitment to coastal areas with adequate conditions for feeding, resting and development [2]. Understanding how populations connect and how animals disperse from rookeries to feeding areas is a challenging task, but is essential for setting priorities and defining management strategies for conservation [3]. In this context, molecular genetic data have been fundamental in obtaining relevant information on interpopulational connectivity, migrations and natal origins of marine turtle feeding populations [4][5][6][7].
Feeding grounds are generally composed of individuals originating from a mixture of sources, being therefore known as ''mixed stocks''. Mitochondrial DNA (mtDNA) haplotype frequencies can be used to determine connections between sea turtles at a feeding ground to their areas of origin (rookeries) by Mixed Stock Analysis (MSA), which uses Markov Chain Monte Carlo (MCMC) sampling to estimate the rookery origins of individuals [8]. ''Many-to-many'' MSA simultaneously estimates origins and destinations of a metapopulation of multiple feeding grounds and rookeries [9], providing a more comprehensive understanding of how areas are linked, and the possible movements that animals undertake. It is usually accepted that these movements are influenced by ocean currents: dispersal at the initial epipelagic phase is thought to be mainly shaped by currents due to the low swimming capability of hatchlings [1,10]; as animals grow and attain more autonomous movement, it is believed that the influence of currents on migrations weakens, but still occurs [11]. The dispersal of post-breeding sea turtles could in fact reflect their previous drift scenarios as hatchlings, with prevailing ocean currents around nesting areas possibly determining how adult turtles select their foraging sites [12]. Multidisciplinary approaches using oceanographic and genetic information are being increasingly applied to studies of the dispersal and migration patterns of sea turtles (e.g. [7,[13][14][15]).
Understanding migratory pathways and origins of animals at their non-reproductive stages is important for determining how impacts such as direct harvest, habitat degradation, and bycatch in fisheries can be shared by seemingly separate populations that are in fact connected [3]. This is of special significance for the critically endangered [16] hawksbill turtle (Eretmochelys imbricata), since the tortoiseshell commerce is still debated by the Convention on International Trade in Endangered Species (CITES) [17]. Due to the uncertainty of how a harvest/other impacts in one area will affect other rookeries and feeding areas, effective conservation and/or sustainable use strategies must be based on the extent of connectivity between hawksbill populations; this type of research is therefore a priority for this highly migratory and endangered species [18].
Important hawksbill nesting and feeding areas exist along the Brazilian coast [19][20][21][22][23], but the connections between them are still largely unknown [24]. Here, we decrease this gap by describing 740 bp mtDNA haplotypes of hawksbills from five Brazilian feeding aggregates, assessing genetic structure, and estimating natal origins through many-to-many MSA analysis. In addition, we verify the influence of ocean currents on connectivity, inferring origins by analyzing surface drifter data for the Atlantic.

Ethics statement
According to Normative Instruction 154/March 2007, all capture, tagging, sampling and transport of biological samples of wild animals for scientific purposes must have approval from Instituto Chico Mendes de Conservação da Biodiversidade (ICMBio) SISBIO committees. This study was approved by the Instituto Chico Mendes de Conservação da Biodiversidade, and conducted under SISBIO licenses #225043, #14122, and #159622. All animal handling was performed by trained personnel, following widely accepted and ethical protocols described in [25]. When capturing live turtles, the following measures were taken to alleviate stress: 1) turtles were kept out of the water for a maximum of ten minutes; 2) work was performed in a shaded area; and 3) animals were released at the same location of capture.

Sampling and laboratory analyses
Skin samples were obtained from 157 immature hawksbills at five areas in Brazilian waters ranging from equatorial to temperate latitudes: (1) São Pedro and São Paulo Archipelago (SPSP, n = 12); (2) Abrolhos National Marine Park (Abrolhos, n = 65), (3) Ceará state coast (Ceará, n = 23); (4) Bahia state coast (Bahia, n = 32); and (5) South Brazil region, which combines the Arvoredo Marine Biological Reserve (n = 6) and Cassino Beach (n = 19). The Arvoredo and Cassino areas were grouped due to small sample size and geographical isolation regarding the other feeding grounds. Two skin samples with approximately 4-5 mm were collected from the fore flippers of each turtle using disposable biopsy punches, and preserved in absolute alcohol. Samples were taken from turtles hand-captured by divers at SPSP, Abrolhos, and Arvoredo, as well as from individuals incidentally caught in fishing nets or stranded on beaches (alive or dead) at Ceará, Bahia, and Cassino Beach. All animals were measured (Curved Carapace Length -CCL), and live turtles were tagged with Inconel Tags (National Tag and Band Co.) using standard techniques [26,27] prior to release. Sampling locations are shown in Figure 1.
Samples were macerated and kept at 37uC in a lysis buffer until complete digestion. DNA was extracted using Genomic DNA Extraction Kits (Norgen Biotek) or a standard phenol:chloroform method [28]. mtDNA control region fragments of approximately 850 bp were amplified via Polymerase Chain Reaction (PCR) using primers LCM15382/H950 [29], as follows: 59 at 94uC; 36 cycles of 300 at 94uC, 300 at 50uC, 19 at 72uC; 109 at 72uC. Samples were purified with Illustra GFX purification kits (GE Healthcare) and sequenced in both directions through capillary electrophoresis using the Applied BiosystemsH 3130 Genetic Analyzer (Valid Biotechnology, Universidade Federal de Minas Gerais, Brazil).

Haplotypes and diversities
Sequences were aligned and cropped to 740 bp using Clustal X 2.0 [30], and classified according to GenBankH and the Atlantic Ocean hawksbill haplotype database (under construction, Abreu-Gobrois pers. comm.). We used 740 bp haplotypes for genetic structure analyses between Brazilian populations, but for comparison with all other Atlantic Ocean feeding grounds and rookeries described in literature, a second cropping to 382 bp (original classification, [31]) was necessary. Arlequin 3.5 [32] was used to assess haplotype (h) and nucleotide (p) diversities of the study areas.

Genetic structure
Genetic divergence between feeding populations was verified in Arlequin 3.5 through pairwise fixation indices F-st and Q-st, respectively using haplotype frequencies only and a Tamura-Nei model of nucleotide substitution (as determined through jMo-delTest 0.1.1, [33]). Besides our study areas, the following feeding aggregations were included in genetic structure analyses: (1)  Figure 1). Analysis of Molecular Variance (AMOVA) was implemented to verify ocean basin structure by grouping Atlantic feeding grounds into Brazilian, African, Caribbean and Gulf of Mexico regions. To identify regional structure within Brazilian waters, we divided feeding grounds into two groups: areas influenced mainly by the South Equatorial/ North Brazil Current (SPSP, Noronha and Ceará) and those influenced by the Brazil Current (Bahia, Abrolhos and South Brazil).

Natal origins
A many-to-many MSA using the ''mixstock'' package in R [9,39] was performed to estimate source contributions to the study areas in Brazil (feeding ground-centric), as well as the destinations of animals hatched at Brazilian rookeries (rookery-centric). We included all above-cited Atlantic foraging grounds (18 areas, n = 1361 individuals) and thirteen rookeries (n = 875) in this metapopulational analysis: (1) Bahia/Sergipe coasts (referred to as Bahia, n = 92); (2) Pipa Beach (referred to as Pipa, n = 27);  [4,5,35,38,[40][41][42][43][44], locations shown in Figure 1). Nesting population size (estimated number of females per nesting season) of each area was included in this analysis as an ecological covariate, following [9]. One MCMC was implemented for each rookery, with chain lengths of 20000, and one half chain discarded as ''burn-in''. Gelman-Rubin convergence factors varied from 1.0 to 1.04 (average 1.02), indicating convergence. Hybrid haplotype EixCc BR3, from the Bahia rookery, was maintained in the MSA since it was also observed at feeding areas; orphan haplotypes were excluded.
Lagrangian drifter data was also used to infer on the possible role of ocean currents on the dispersal of hatchlings, following [15]. For such, we downloaded surface drifter data freely available from NOAA's Global Drifter Program (www.aoml.noaa. gov/ envids/gld). Areas with size 4u64u (latitude and longitude) were delineated around all rookeries considered in MSA (n = 13), and the number of drifters that passed through them and reached the Brazilian foraging grounds was counted (drifters with less than three months of transmission were excluded). Feeding areas with similar drifter patterns were grouped into target regions for simplification: 1) SPSP; 2) Noronha and Ceará; and 3) Bahia, Abrolhos and South Brazil. Based on these drifter counts and rookery population size information (females/nesting season), the probability that a turtle arriving at the target region originated from a previously characterized rookery was calculated under a Bayesian framework (for method details refer to [15]).
Finally, natal origin estimates inferred by genetic/rookery size and drifters/rookery size data were compared to better assess the influence of ocean currents on feeding grounds aggregations in Brazil. We evaluated this correlation through a Mantel test implemented with the package ''vegan'' in R [39,45], and performed linear regression by regressing log-transformed proportions between genetic and drifter profiles.

Genetic structure
We detected overall significant genetic difference between and within nesting and feeding populations in the Atlantic Ocean (p,0.001). F-st and Q-st values show that feeding aggregations are not homogenous throughout the Atlantic (Table 3), with significant regional differences: AMOVA confirmed a pronounced structure between Brazilian, African and Caribbean feeding populations (F-st = 0.618, p,0.001). Structure was also high when dividing feeding areas into Brazil, Africa, Caribbean and Gulf of Mexico groups (F-st = 0.553, p,0.001). Within Brazilian feeding aggregations, some differences were observed (see Table 3), with two areas being more differentiated from the remainder: the oceanic islands of Fernando de Noronha/Rocas Atoll, off the northeastern coast of Brazil, and the southernmost hawksbill occurrence area, South Brazil. Structure was significant (Fst = 0.581, p,0.001) between areas under influence of the South Equatorial/North Brazil Current (SPSP, Noronha and Ceará) and those influenced by the Brazil Current (Bahia, Abrolhos and South Brazil).

Natal origins
Many-to-many MSA estimates suggest that the feeding aggregations in Brazil are mostly composed of animals originating from domestic rookeries, but also present some distant origins (Figure 3a). The Brazilian rookeries Bahia and Pipa were the predominant sources, with respective contributions of approximately 12% and 28% for SPSP, 30% and 22% for Noronha, 22% and 18% for Ceará, 36% and 17% for Bahia, 21% and 28% for Abrolhos, 34% and 30% for South Brazil. Rookeries from Africa and the Caribbean showed generally lower but noteworthy contributions: Principe, Barbados and Cuba contributed respectively 8%, 10% and 14% to Noronha; Barbados, Puerto Rico and Cuba 12%, 9% and 13% to SPSP; Barbados and Cuba 12% and 15% to Ceará, 9% and 10% to Bahia, 8% and 13% to Abrolhos.
Rookery-centric MSA results for the Brazilian rookeries of Bahia and Pipa shows that these nesting areas contribute mainly to the domestic feeding populations: respectively 6% and 16% to SPSP; 12% and 9% to Ceará; 20% and 10% to Bahia; 11% and 17% to Abrolhos; 18% and 18% to South Brazil; and 10% and   12% to Noronha (Figure 3b). All other contributions were lower than 5%. MSA outputs for overseas feeding grounds and rookeries (not discussed in this study) are presented in Figures S1 and S2.
Of the available drifter data, a total of 469 drifters passed through the Atlantic rookeries, of which 388 transmitted for over three months. Of these, 37 drifters arrived at our three target areas in Brazil, originating from the Bahia, Pipa and Principe rookeries (trajectories shown in Figure 4). Displacements from rookeries to Brazilian feeding areas were most likely by means of the North Brazil Current for drifters leaving the Pipa rookery, North Brazil/ Brazil Current for those leaving Bahia, and South Equatorial Current for drifters from Principe. No drifters from the Caribbean arrived at the Brazilian coastline (see Figure S3). When estimating natal origins of hawksbills at our target areas through drifter/ population size data, significant contributions were: Pipa (58%), Bahia (9%) and Cuba (11%) for SPSP; Pipa (62%), Bahia (14%) and Cuba (10%) for Noronha/Ceará; and Bahia (72%) for Abrolhos/Bahia/South Brazil (Table 4). Correlation between natal origins as calculated by MSA (which included genetic/ rookery size data) and drifters/rookery size information was significant (Mantel test, r = 0.791, p,0.05; linear regression r = 0.560, p,0.01; Figure 5).

Haplotypes and diversities
As illustrated in Figure 2, several haplotypes observed at Brazilian feeding areas are common at rookeries located in the Caribbean (A01/A, A09/F, A11/F, A24/Q), Brazil (A32, A62, EixCcBR3) and Africa (EATL). We also detected two orphan haplotypes (i.e. not seen at rookeries): A76 and A92 (the latter previously undescribed); as noted in [35,58], the presence of orphan haplotypes indicates the need for further sampling of rookeries. Mean haplotype diversity of feeding areas in Brazil was 0.413, which is higher than the mean diversity of the two described feeding grounds off West Africa (h = 0.335), but lower than the diversity of ten populations in the Caribbean (h = 0.647). These low diversities could be a result of natal origins from fewer, less diverse sources when compared to the Caribbean, where a large number of rookeries with relatively high diversities are present. SPSP and Noronha, located respectively around 350 and 1000 km off northeast Brazil, presented the highest haplotype diversities within the sampled Brazilian feeding grounds. Bass et al. (2006) suggest that feeding areas located at the confluence of several current systems present higher haplotype diversities [59]; Vilaça et al. (2013) propose that Noronha's location, influenced by both the North Brazilian and South Equatorial Currents, is a possible explanation for its elevated diversity [24]. Our results support the hypothesis that ocean currents influence diversity since SPSP, located near Noronha, showed highest diversity values among all Brazilian feeding grounds.

Genetic structure
Our study suggests regional genetic structuring within the Atlantic, between (i) foraging areas off Brazil, Africa and the Caribbean, a pattern that has been previously noted [24], and (ii) feeding areas in the Gulf of Mexico and the remainder within the Caribbean. Although turtles disperse widely, they can be constrained to certain regions due to factors such as geographical distance, rookery population size, animal behavior and surface ocean currents. For instance, the Gulf of Mexico presents a distinct current pattern (as shown in Figure S3) that could retain animals within the region. In a general manner differentiation between Brazilian feeding grounds was not pronounced, but significant structuring was observed when grouping Brazilian populations into groups influenced by different ocean currents. This suggests that migration patterns are somewhat limited, and that genetic composition could be affected by the different factors noted above. We highlight here the importance of analyzing longer haplotype sequences, which present higher resolution and better detect population structure [40,60]. Although longer haplotypes are being increasingly analyzed, past hawksbill turtle studies mostly characterize 382 bp fragments; it is important that future studies focus on larger fragments to improve our understanding on structure and connectivity of hawksbill populations. Table 3. Pairwise F-st (above diagonal) and Q-st (below diagonal) values between feeding populations in the Atlantic .   SP  CE  BA  AP  SB  NR  PP  CV  UV  PR  DR  TC  BH  CU  CY  MX  TX  Natal origins Many-to-many MSA showed that feeding aggregations in Brazil are generally composed of animals from limited sources, being mostly linked to the proximal Brazilian rookeries. Results also indicated that in some cases connectivity extends to rookeries located in the Caribbean and West Africa. Origins calculated from drifter data support main contributions from Bahia and Pipa, with almost all drifters arriving from national rookeries to the Brazilian feeding grounds, via the North Brazil and Brazil currents. The only other work to describe origins of hawksbills at a Brazilian feeding ground (Noronha/Rocas) also shows that main origins are from national sources, with lower contributions from other regions [24].
Our rookery-centric MSA, including the five populations investigated in this study plus the Noronha/Rocas feeding ground, decreased contributions of the Bahia rookery to ''unknown'' feeding grounds from almost 60% (reported in [24]) to less than 10%, and showed strong connections to feeding grounds in national waters. This is likely due to the detection at our feeding areas of a hybrid hawksbill x loggerhead haplotype (EixCc3 BR3), present in 18% of samples from the Bahia rookery [41]. This haplotype had not been registered at feeding grounds sampled by [24], but we report its occurrence at Ceará (4% frequency) and in the temperate waters of South Brazil (16% frequency) [61]. This demonstrates the importance of thorough genetic analyses of populations and the use of metapopulational approaches (ocean basin scale) in animal populations that are demographically connected between distant geographic localities.
Natal origin estimates also showed connection between West Africa and Brazil, with 8% contribution from Principe to Noronha. This is expected since the EATL haplotype, present with 100% frequency at the Principe rookery [35], was encountered at this feeding area. However, drifter tracks actually show a higher connection with West Africa by means of the South Equatorial Current; for example, some drifters at SPSP originated from Principe, but MSA showed no contribution from this rookery to this region as haplotype EATL was not detected. The displacement to the opposite direction via the Equatorial Counter Current is also likely as shown through drifter tracks. Furthermore, hawksbill tag returns have repeatedly evidenced transatlantic movements between these regions [62,63]. Increasing the number of analyzed samples from SPSP and other feeding grounds along the coast would possibly lead to the detection of the EATL haplotype and further evidence the connection between African and Brazilian hawksbill turtle populations.
MSA showed that Caribbean rookeries contributed to all Brazilian feeding grounds except South Brazil. Highest contributions were observed for SPSP, Noronha, and Ceará, which aggregated animals from up to three Caribbean rookeries. It is hypothesized that as sea turtles grow their feeding habitats become preferentially established near their natal beach [4]. Immature hawksbills sampled at SPSP and Noronha presented larger size classes ( Table 2) than those at the other feeding areas in Brazil, possibly indicating that some animals are moving closer to the Caribbean for the onset of reproduction.
The connection between Caribbean and Western South Atlantic green and hawksbill turtle populations has been extensively shown through MSA [13,15,24,64] and tag returns [65]. However, oceanographic data have yet to bring additional evidence to these links. In this study no drifters from Caribbean sources arrived in Brazilian waters, although the opposite displacement occurred. Drifter data are not abundant for every region and the life span of drifters could be too short (,291-450 d; [66]) to detect this connection. Particle modeling could solve the temporal limitation of drifters, but a study using particle tracking models (dispersal time = one year) also did not reveal any virtual drifter from Caribbean hawksbill rookeries arriving around Brazil [7]; a 5-year particle tracking study from the green turtle rookery Costa Rica also did not show this connection [13]. Contribution from the Caribbean to Brazil could be overestimated by MSA due to shared haplotypes between rookeries of these regions; on the other hand, MSA estimates might be correct but animals are undergoing movements undetected by drifters/ particle models (e.g. transport by coastal currents in shallow areas or active swimming). Additional genetic characterization of rookeries and feeding areas, as well as further description of animal behavior and oceanographic features on the routes between them, is necessary for resolving this ambiguity.
Although surface drifter patterns did not corroborate MSA origins in all cases, origins calculated by combining drifter data with population size of rookeries (see Table 4) were strongly correlated with MSA estimates. This indicates that ocean currents influence how patterns of hawksbill population genetic structure are shaped in Brazil. Ocean currents apparently play a substantial role in shaping populations in other regions as well; for example, a particle tracking study found significant correlation between particle distribution patterns and natal origin estimates for foraging aggregations in the Caribbean [7]. Restricting MSA to only small turtles, more likely to have recently recruited than larger juveniles that may have already conducted developmental migrations [67], could increase the correlation between natal origins and ocean currents. For allowing this type of meta-analysis, researchers should publish their datasets containing size/haplotype of individual turtles. In light of this suggestion, data from the present work was published in Figshare [68].
The integration of genetic analysis and oceanographic data is a valuable approach to understanding origin and distribution of immature animals; however, there are some caveats inherent to the analyses applied here. Mixed stock analysis may not adequately detect connections between some areas since not all source rookeries are thoroughly characterized, while those that are do not always present highly differentiated haplotype frequencies, leading to high uncertainty [58]. Furthermore, drifter trajectories might not represent the precise pathways of sea turtles, as analyses cannot be limited to hatching seasons due to large reduction of data, and drifters generally have short life spans and do not consider surface wind drag and turtle behavior/swimming [66]. Despite these caveats, the type of information described above has proven to be quite useful for indicating possible movements and migrations between areas (e.g. [10,15,69,70]). Ocean circulation models are also being increasingly used, and allow choosing the number and lifespan of particles, release periods, and incorporating additional factors such as wind drag and swimming power of hatchlings [71,72]. However, they are also limited due to uncertainties associated with simulation of turtle behavior and spatio-temporal resolution unable to capture fine-scale features [73]. An ideal scenario would integrate several connectivity indicators; for achieving this, other means of linking populations should continue to be investigated. Stable isotope analysis, for example, is currently being used to distinguish populations and unveil migratory origins by comparing isotopic signatures of animals to predominant regional isoscapes [74][75][76].

Implications for hawksbill conservation
How populations connect, and how impacts on any particular population will affect others, are fundamental questions for the conservation of endangered, highly migratory animals. Hawksbill turtles are currently listed as critically endangered in the IUCN red list [16], and listed on Appendix I of the Convention on the International Trade in Endangered Species. Nevertheless, hawksbill products continue to be commercially valuable and in demand in some regions, and debate has arisen as to whether or not harvesting should be permitted, and at what levels such harvest would be sustainable [4,[77][78][79]. Some studies have tackled these issues and concluded, for example, that exploiting turtles in Caribbean feeding areas would affect rookeries throughout the entire region [4,79].
Although Brazilian populations have also been historically depleted due to harvest for consumption and fisheries bycatch, nesting populations are currently showing encouraging increases [22,23]. Our analyses indicate that mixed stocks are composed mostly of animals that originate from Brazilian rookeries, which is a fortunate scenario since a large sea turtle conservation project (Project Tamar-ICMBio) has been underway in Brazil since 1980 [80] and the need for transnational conservation policies is somewhat reduced. Despite this encouraging outlook, it is not yet known if immature feeding population numbers are also rising, and impacts such as bycatch [81] and ingestion of/entanglement in plastic debris [82] are increasingly causing sea turtle mortality. Furthermore, we demonstrate that some areas are demographically and genetically connected to rookeries in the Caribbean and West Africa. For instance, the SPSP and Noronha feeding grounds aggregate large animals from various national and international origins, and impacts at these developmental areas could affect several nesting populations throughout the Atlantic. Although Fernando de Noronha, Rocas Atoll and Arvoredo Reserve are protected areas where fishing is not allowed, other areas, including the waters around the São Pedro and São Paulo Archipelago, are important fishing grounds [83]. The Ceará, Bahia and Cassino Beach areas are also not integrally protected, and impacts at these areas will most likely influence the Brazilian nesting populations, and possibly some Caribbean rookeries.
In this study we reduce the gap in genetic information of hawksbill turtle feeding populations in Brazil, and illustrate how these populations are connected with proximal and distant rookeries in the Atlantic Ocean. The results presented here should be considered in future national and multinational strategies for mitigating bycatch/other impacts to increase protection and recovery of endangered populations of hawksbill turtles. Supporting Information