Proceedings of the Royal Society B: Biological Sciences
You have accessResearch articles

Green, yellow or black? Genetic differentiation and adaptation signatures in a highly migratory marine turtle

Rocío Álvarez-Varas

Rocío Álvarez-Varas

Departamento de Ciencias Ecológicas, Facultad de Ciencias, Universidad de Chile, Santiago, Chile

Núcleo Milenio de Ecología y Manejo Sustentable de Islas Oceánicas (ESMOI), Departamento de Biología Marina, Universidad Católica del Norte, Coquimbo, Chile

Qarapara Tortugas Marinas Chile NGO, Santiago, Chile

[email protected]

Google Scholar

Find this author on PubMed

,
Noemi Rojas-Hernández

Noemi Rojas-Hernández

Departamento de Ciencias Ecológicas, Facultad de Ciencias, Universidad de Chile, Santiago, Chile

Contribution: Data curation, Formal analysis, Methodology, Software, Validation, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Maike Heidemeyer

Maike Heidemeyer

Centro de Investigación en Biología Celular y Molecular (CIBCM), Universidad de Costa Rica, San José, Costa Rica

Contribution: Data curation, Funding acquisition, Investigation, Resources, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Cynthia Riginos

Cynthia Riginos

School of Biological Sciences, The University of Queensland, Brisbane, Australia

Contribution: Conceptualization, Investigation, Methodology, Software, Supervision, Validation, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Hugo A. Benítez

Hugo A. Benítez

Laboratorio de Ecología y Morfometría Evolutiva, Centro de Investigación de Estudios Avanzados del Maule, Universidad Católica del Maule, Talca, Chile

Google Scholar

Find this author on PubMed

,
Raúl Araya-Donoso

Raúl Araya-Donoso

School of Life Sciences, Arizona State University, Tempe, AZ, USA

Contribution: Formal analysis, Methodology, Software, Validation, Visualization, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Eduardo Reséndiz

Eduardo Reséndiz

Departamento Académico de Ciencias Marinas y Costeras, Universidad Autónoma de Baja California Sur, Mexico

Google Scholar

Find this author on PubMed

,
Mónica Lara-Uc

Mónica Lara-Uc

Departamento Académico de Ciencias Marinas y Costeras, Universidad Autónoma de Baja California Sur, Mexico

Contribution: Funding acquisition, Investigation, Resources, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Daniel A. Godoy

Daniel A. Godoy

Coastal-Marine Research Group, Institute of Natural and Mathematical Sciences, Massey University, Auckland, New Zealand

Contribution: Funding acquisition, Investigation, Resources, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Juan Pablo Muñoz-Pérez

Juan Pablo Muñoz-Pérez

Galapagos Science Center GSC (Universidad San Francisco de Quito USFQ-University of North Carolina at Chapel Hill UNC), Isla San Cristobal, Galápagos, Ecuador

University of the Sunshine Coast USC, 90 Sippy Downs Dr, Sippy Downs, Queensland 4556, Australia

Google Scholar

Find this author on PubMed

,
Daniela E. Alarcón-Ruales

Daniela E. Alarcón-Ruales

Galapagos Science Center GSC (Universidad San Francisco de Quito USFQ-University of North Carolina at Chapel Hill UNC), Isla San Cristobal, Galápagos, Ecuador

Google Scholar

Find this author on PubMed

,
Joanna Alfaro-Shigueto

Joanna Alfaro-Shigueto

ProDelphinus, Lima, Peru

Facultad de Biología Marina, Universidad Científica del Perú, Lima, Peru

Contribution: Funding acquisition, Investigation, Resources, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Clara Ortiz-Alvarez

Clara Ortiz-Alvarez

ProDelphinus, Lima, Peru

Contribution: Funding acquisition, Investigation, Resources, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Jeffrey C. Mangel

Jeffrey C. Mangel

ProDelphinus, Lima, Peru

Contribution: Funding acquisition, Investigation, Resources, Writing-review & editing

Google Scholar

Find this author on PubMed

,
Juliana A. Vianna

Juliana A. Vianna

Departamento de Ecosistemas y Medio Ambiente, Facultad de Agronomía e Ingeniería Forestal, Pontificia Universidad Católica de Chile, Santiago, Chile

Contribution: Conceptualization, Investigation, Supervision, Validation, Writing-review & editing

Google Scholar

Find this author on PubMed

and
David Véliz

David Véliz

Departamento de Ciencias Ecológicas, Facultad de Ciencias, Universidad de Chile, Santiago, Chile

Núcleo Milenio de Ecología y Manejo Sustentable de Islas Oceánicas (ESMOI), Departamento de Biología Marina, Universidad Católica del Norte, Coquimbo, Chile

[email protected]

Google Scholar

Find this author on PubMed

Abstract

Marine species may exhibit genetic structure accompanied by phenotypic differentiation related to adaptation despite their high mobility. Two shape-based morphotypes have been identified for the green turtle (Chelonia mydas) in the Pacific Ocean: the south-central/western or yellow turtle and north-central/eastern or black turtle. The genetic differentiation between these morphotypes and the adaptation of the black turtle to environmentally contrasting conditions of the eastern Pacific region has remained a mystery for decades. Here we addressed both questions using a reduced-representation genome approach (Dartseq; 9473 neutral SNPs) and identifying candidate outlier loci (67 outlier SNPs) of biological relevance between shape-based morphotypes from eight Pacific foraging grounds (n = 158). Our results support genetic divergence between morphotypes, probably arising from strong natal homing behaviour. Genes and enriched biological functions linked to thermoregulation, hypoxia, melanism, morphogenesis, osmoregulation, diet and reproduction were found to be outliers for differentiation, providing evidence for adaptation of C. mydas to the eastern Pacific region and suggesting independent evolutionary trajectories of the shape-based morphotypes. Our findings support the evolutionary distinctness of the enigmatic black turtle and contribute to the adaptive research and conservation genomics of a long-lived and highly mobile vertebrate.

1. Introduction

In the absence of physical barriers to gene flow, highly mobile marine species are expected to exhibit genetic homogeneity over long distances [1,2]. However, other studies provide evidence of strong population genetic structure, often accompanied by phenotypic variation [3,4]. Such patterns can be related to oceanographic factors, life-history traits and adaptation [14]. Sea turtles are migratory species with a complex life cycle that includes foraging grounds (FGs) geographically distant from breeding areas [5]. The green sea turtle (Chelonia mydas) has a global distribution and is characterized by strong female nesting site fidelity, moderate male-mediated gene flow and population overlap during migration and within FGs, leading to complex genetic structure patterns in these areas [69].

The formation of the Panama Isthmus separated Atlantic and Pacific C. mydas populations about 3.5 Ma [10,11]. In the Pacific Ocean, an intraspecific clade from the north-central and eastern Pacific regions (NC/EPGL) split from the south-central western Pacific lineage (SC/WPGL) around 0.34 Ma [9,11]. Recent studies using geometric morphometrics showed these main genetic lineages (SC/WPGL and NC/EPGL) vary significantly in their phenotypic traits, including the shape of the head, carapace, plastron and flipper [9,12]. These shape-based morphotypes (known as south-central/western and north-central/eastern Pacific morphotypes) [9,12] also seem to fit with the two widely recognized morphotypes based on coloration: yellow (south-central/western morphotype) and black/melanic (north-central/eastern morphotype) [1214]. Previous studies based on mitochondrial DNA (mtDNA) point to the NC/EPGL as a monophyletic group but included in the SC/WPGL (paraphyletic group) [9,11]. Specifically, the black or melanic morphotype, also known as ‘black turtle’, seems to be part of the C. mydas north-central/eastern Pacific clade restricted to rookeries of the eastern Pacific region (EPO) [9,12]. Although ecological and genetic studies have suggested spatial segregation of rookeries between the Pacific morphotypes [9,11,14], they feed sympatrically in certain FGs such as New Zealand, Costa Rica and Galapagos [8,9].

The eastern Pacific region (EPO) is located along the west side of the American continent between the north and south Pacific subtropical gyres; it is dominated by the presence of two cold water current upwelling systems, the California and Humboldt Currents [15]. The EPO comprises a unique marine ecosystem characterized by a high rate of inter-annual (i.e. El Niño Southern Oscillation) and multi-decadal climate variability, and resources that vary temporally and spatially [15,16]. The surface water salinity of this region is higher than in the western region due to a combination of factors including atmospheric convection, precipitation, evaporation and ocean dynamics [17]. This high climatic and oceanographic variability has a considerable effect on the population dynamics of a wide range of trophic levels within the food web [15]. Sea turtle species inhabiting the EPO have developed distinct biology, morphology and behaviour compared to their counterparts elsewhere [15,16].

There is considerable evidence indicating that besides having a distinctive body shape [9,12,18] and a darker coloration of carapace and plastron [13,14], the C. mydas EPO populations also show smaller size of nesting females [19], lower reproductive output (number of eggs per clutch) [15], different reproductive periodicity [16,19], lower somatic growth [14] and overwintering torpor during cold water periods [16,20]. Although this species has been widely recognized as herbivorous, EPO populations seem to have a predominantly omnivorous diet [16,21]. The specific oceanographic and environmental features of the EPO and the biological and ecological distinctiveness of the black turtle suggest they may be reproductively isolated from other C. mydas groups and uniquely adapted to the new colonized environment (EPO). Nevertheless, whether the black turtle is evolutionary distinct and whether its genome exhibits adaptation signatures is unresolved among the scientific community and sea turtle conservationists.

During recent years, adaptive research has experienced a growing increase owing to the development of new genome technologies and greater power of bioinformatics analysis [2,22]. Most genomic studies in marine systems have been oriented to species with short generational times [22,23], leading to a lack of information on long-lived and highly mobile marine organisms. Particularly in sea turtles, the field of genetics has been focused on the use of mitochondrial and microsatellite markers [5], and although there are recent studies based on genomic techniques (e.g. transcriptomics [24]; RAD-capture [25]; GBS [26]), to date there are no studies published on adaptive research in C. mydas.

Conservation and effective management of marine biodiversity require knowledge of the geographical extent of populations and understanding of the population connectivity and life histories [9]. To comprehend the processes that underlie the differences between shape-based morphotypes of the Pacific C. mydas (e.g. historical, oceanographic, ecological and adaptive processes, life-history traits), we assess the genetic structure among morphotypes using neutral loci, and also identify candidate loci to be under selection linked to biological functions. We tested whether there are genomic regions associated with cold resistance, pigmentation, osmoregulation, morphogenesis, digestive metabolism and reproductive activity that give insights into the adaptation of this species to contrasting environmental conditions (i.e. EPO versus the rest of the Pacific Ocean). Our approach contributes to conservation genomics based on preserving the genetic and adaptive evolutionary potential in endangered cosmopolitan species.

2. Material and methods

(a) Sampling and DNA extraction

A total of nine C. mydas FGs were sampled in this study, including two locations in the western Pacific Ocean (New Zealand and Fiji) and seven in the eastern Pacific (EPO; including Mexico, Costa Rica, Galapagos-Ecuador, Peru and Chile; figure 1; electronic supplementary material, table S1; see more detail in [9]). Five of these FGs host both shape-based morphotypes: New Zealand, Costa Rica, Galapagos–Ecuador, Peru and Easter Island in Chile [9]. Blood and skin samples were collected and preserved in 90% ethanol or saturated salt solutions. DNA was isolated using the salting-out protocol described by Aljanabi & Martinez [27]. DNA quality and concentration were determined with a specific fluorimetry method using a Qubit fluorimeter (Life Technologies, NY) and agarose gel electrophoresis (1% with GelRed Nucleic Acid Stain-Biotium). Samples were sequenced and genotyped using DArTseq TM genotyping technology at Diversity Arrays Technology, Australia [28].

Figure 1.

Figure 1. Map depicting the Chelonia mydas foraging grounds located in the Pacific Ocean included in this study. Colours in the pie graphs represent the proportions of each shape-based morphotype in foraging grounds. South-central/western Pacific morphotype: putative yellow morphotype and north-central/eastern Pacific morphotype: putative black/melanic morphotype; NZ, New Zealand; FI, Fiji; EI, Easter Island (Chile); ME, Mexico; CR, Costa Rica; GA, Galapagos (Ecuador); PE, Peru; AR, Arica (Chile); AT, Atacama (Chile). (Online version in colour.)

(b) Library preparation and sequencing

Samples were digested (150 ng of gDNA) with a combination of both a frequent and rare cutting methylation-sensitive restriction enzyme (PstI and SphI, respectively) [28]. Small fragments (less than 200 bp) of digested DNA were ligated to a barcoded adaptor (6–9 bp length) and amplified using PCR. PCR products were standardized in concentration and pooled for sequencing on a single HiSeq 2500 (Illumina) lane at a 2.5 million read depth (a plate with 94 samples was pooled in each lane). The mean read depth was 10.719, with a minimum value of 4.1 and a maximum of 76.1.

(c) Quality control and initial SNP calling

FASTQ files were further filtered to produce DArT scores and SNP reports using the DArTtoolbox bioinformatics pipeline [29]. Loci were identified as SNP or reference alleles according to the occurrence frequency. Further filtering was conducted using the dartR package implemented in R software [30]. This filtering used the following criteria: reproducibility (threshold = 1.00), minor allele frequencies (MAF) >0.01, individual call rate (proportion with non-missing scores for each individual, removing those individuals below a specified threshold) with threshold >0.90, locus call rate (proportion with non-missing scores for each locus, removing those loci below a specified threshold) with threshold >0.85, and discarding monomorphic markers. Fragments containing more than one SNP were filtered using the gl.filter.secondaries command.

(d) Dataset assignment

Given that C. mydas FGs host individuals with multiple origins [9], a sequential approach was used to define datasets in this study. First, we used an individual-based approach (without prior information) with all individuals and loci to detect how our samples were grouped (PCoA, DAPC and STRUCTURE; see below). Then we defined the genetic groups, and these groups were used for the detection of outliers. Neutral loci were used to estimate C. mydas genetic diversity and population structure, while outlier SNPs were used to evaluate adaptation signatures (candidate loci). The datasets were analysed separately. Samples were not grouped according to gender due to the difficulty of sexing juveniles in the field. Nevertheless, due to possible sex-biased dispersal in C. mydas, it is recommendable to integrate gender information in future studies to examine whether it affects the population genetic structure [31].

(e) Individual-based approach using all individuals and loci

A Euclidian distance-based principal coordinates analysis (PCoA) was carried out with the dartR package, in order to observe the genetic relationships among individuals. Discriminant analysis of principal components (DAPC) was also performed as an initial analysis of genetic structure using the find.clusters command of the ADEGENET package implemented in R software [32]. Finally, to find the most probable number of genetic clusters (K) and proportions of genotype admixture we used a Bayesian approach implemented in STRUCTURE V. 2.3.4 [33]. This analysis comprised three independent runs performed with a burn-in of 100 000 steps, followed by 200 000 additional Markov chain Monte Carlo (MCMC) iterations. An admixture ancestry model was assumed with independent allele frequencies and no population priors. LnP (D) values obtained for each K value were compared using Structure Harvester software [34].

Given that the results of the individual-based approach grouped samples according to the shape-based morphotypes previously defined by Álvarez-Varas et al. [9,12] (north-central/eastern and south-central/western Pacific morphotypes; electronic supplementary material, figure S1 and S2; see Results), all outlier analyses were carried out using these groups. For the purposes of this paper, the north-central/eastern Pacific morphotype was referred to as ‘black turtle’ and the south-central/western Pacific morphotype as ‘yellow turtle’.

FGs of Chelonia mydas are composed of mixed populations [9,12], which was corroborated in our individual-based analysis; grouping was not according to FGs (electronic supplementary material, figure S2b). For this reason, FGs were only considered as study sites (not as genetic groups/clusters). We considered the genetic groups as populations for the analyses described below.

(f) Genetic diversity and population structure among morphotypes based on neutral SNPs

After removing outliers as described below, we eliminated from the remaining neutral dataset all loci showing significant departures from Hardy–Weinberg equilibrium (HWE) using the gl.filter.hwe command (alpha = 0.05). Loci showing significant linkage disequilibrium (LD) (R2 = 0.20) were also eliminated using PLINK V. 1.9 [35].

Using this filtered neutral dataset (9473 SNPs), PCoA, DAPC and STRUCTURE analyses were carried out with the parameters mentioned above (i.e. the individual-based approach). Pairwise FST were calculated (between morphotypes) using the gl.fst.pop command of the dartR package; significance was assessed by 1000 bootstrap replicates. We inferred the past demography for each group separately using Stairway Plot 2 [36], considering a mutation rate of 1.2 × 10−8 substitutions per site per My [37]. Finally, for each group detected by these analyses, the expected heterozygosity (He) and observed heterozygosity (Ho) were calculated using Genetix V. 4.05 [38]; the number of private alleles (Ap) and the number of rare alleles (Ar; MAF < 5%) were estimated with the dartR package.

(g) Detection of outlier SNPs and regions potentially under selection

Although several tests for the identification of outlier loci have been developed, they may not account for population structure and its covariance with other variables (such as demography and mutation rates) accurately, which could lead to false positives [39]. To resolve this issue, three methods based on different models were run in the software BayeScan, Fsthet and PCAdapt. BayeScan compares empirical FST values with a null distribution, determining the likelihood of FST values to be greater or lower than the FST under neutrality [40]. The following parameters were used: 20 short pilot runs with 5000 integration, burn-in set to 5 × 104 and thinning interval and prior odds of 10. A false discovery rate (FDR) of 0.05 was used to select the candidate loci. The Fsthet package implemented in R software uses the raw empirical data to generate smoothed outlier plots for the FST–heterozygosity relationship [41]. A significance level of 0.05 was used to select the candidate loci by this method. Finally, the PCAdapt package, also implemented in R software, detects outlier loci based on PCoA by assuming that markers excessively related to population structure are candidates for potential adaptation [42]. An MAF = 0.01 and a q-value < 0.05 were used to select the candidate loci. Loci that are not significantly detected as outliers using these approaches should not be assumed as neutral loci, since failure of the outlier tests does not demonstrate neutrality [43]. Therefore, a conservative approach was used here leading to three datasets: (i) the outlier loci detected by at least two programmes, (ii) the non-outlier loci (assumed as neutral set), and (iii) a grey zone composed of outlier loci detected only by one programme (so not considered as either outlier or neutral loci). This last dataset was not considered for other analyses.

Reads containing the outlier SNPs were compared against the Chelonia mydas reference genome (CheMyd_1.0 [44]) using the BLASTn tool (National Center for Biotechnology Information, https://blast.ncbi.nlm.nih.gov/Blast.cgi). This analysis was used to detect if the reads corresponded to coding sequences (CDS), introns or intergenic regions (electronic supplementary material, table S2). In cases where the outlier was located within a CDS, we visually compared the nucleotide sequence and its corresponding amino acid translation to assess if the SNP corresponded to a synonymous or non-synonymous mutation (electronic supplementary material, table S2). For outlier sequences that corresponded to genic regions (i.e. introns and CDS), an overrepresentation analysis was performed with g:GOSt in g:Profiler [45] using the Benjamini–Hochberg method (FDR < 0.05) to determine which Reactome, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) biological processes (BP), molecular functions (MF) and cellular components (CC) categories were significantly overrepresented among C. mydas shape-based morphotypes. REViGO was used to visualize and analyse semantically statistically enriched GO terms, including the following parameters: medium similarity (0.7) for semantic overlap, clustered terms by p-value and the default whole UniProt database [46]. To facilitate data interpretation and discussion of results, putative general categories were defined based on the description associated with each of the GO terms (electronic supplementary material, table S3).

3. Results

(a) SNP calling and quality control

A total of 74 339 SNPs were obtained from 187 individuals of Chelonia mydas. Following the first genotype filtering, 10 637 SNPs from 158 individuals were retained. Twenty-nine individuals failed the ‘missing data’ filtering (great than 10%); 19 of which were from Fiji. The failure of these samples probably was due to their low DNA quality. After removing outlier loci and applying filters based on HWE and LD, 9473 loci and 158 individuals were retained. This reduced dataset comprised the neutral SNP dataset used for subsequent genetic diversity and structure analyses.

(b) Individual-based approach

The PCoA using all loci and retained individuals exhibited a clear segregation of the shape-based morphotypes in the first principal component, which accounted for 13.6% of the total genetic variance (electronic supplementary material, figure S1a). Higher dispersion of data was observed for the SC/WPGL in the second component; however, it only explained 1.8% of the total variance (electronic supplementary material, figure S1a). For the DAPC 110 PCs were retained, representing about 80% of the total variance. DAPC supported two well-differentiated groups discriminated by the first component (electronic supplementary material, figure S1b). The STRUCTURE analysis also revealed K = 2 as the most likely number of clusters, with a high probability of assignment of individuals to each cluster (electronic supplementary material, figure S2). These individuals corresponded to the different shape-based morphotypes described by Álvarez-Varas et al. [9,12].

(c) Genetic diversity and population structure

The analyses based on neutral data showed similar results to those obtained using all loci. Clear differentiation between shape-based morphotypes in the first principal component was observed in the PCoA (10.4% of the total genetic variance; figure 2a). Higher dispersion of data was observed for the south-central/western Pacific morphotype in the second component but it only explained 1.7% of the total variance; figure 2a). For the DAPC 110 PCs were also retained, representing about 80% of the total variance. DAPC supported two well-differentiated groups discriminated by the first component (figure 2b). The STRUCTURE analysis also indicated K = 2 as the most likely number of clusters (figure 2c). FST analysis detected significant differences between the groups (FST = 0.108 and p-value < 0.001). The Stairway Plot 2 models showed a tendency toward population decline over 80 kya in both shape-based morphotypes (electronic supplementary material, figure S3).

Figure 2.

Figure 2. Clustering analyses based on the neutral dataset (9473 SNPs) of Chelonia mydas, including individuals from eight Pacific foraging grounds. (a) Principal coordinate analysis (PCoA), (b) discriminant analysis of principal components (DAPC) and (c) STRUCTURE analysis according to the Pacific shape-based morphotypes at K = 2, K = 3 and K = 4. All analyses favoured K = 2 represented by both shape-based morphotypes. (Online version in colour.)

Values of genetic diversity were: expected heterozygosity (He): 0.1568 and 0.1594; observed heterozygosity (Ho): 0.1472 and 0.1446; private alleles (Ap): 1568 and 928, and rare alleles (Ar): 4768 and 4743 for black (north-central/eastern Pacific) and yellow (south-central/western Pacific) morphotypes, respectively.

(d) Candidate SNPs under selection

After initial filtering, 641 outlier loci were detected using Fsthet; 306 with PCAdapt and only six with Bayescan (figure 3a). Therefore, to establish the outlier SNPs dataset six loci shared by the three approaches and 61 shared by Fsthet and PCAdapt were maintained, totalling 67 SNPs (figure 3a). The rest of the outlier loci (detected by only one programme) were classified into the grey zone. Sixty-four of the total outlier SNPs aligned unambiguously against the C. mydas genome. Thirty-two SNPs were localized in genes (electronic supplementary material, table S2; figure 3b) and 35 in intergenic regions. Thirty SNPs situated in genes were in introns and two were localized in CDS; both of the latter were synonymous mutations (electronic supplementary material, table S2). The CREBBP gene exhibited a synonymous TTT/TTC substitution in phenylalanine (F) and the ZNF384 gene showed a CCA/CCG synonymous substitution in proline (P) (electronic supplementary material, table S2).

Figure 3.

Figure 3. Outlier SNPs recovered among Chelonia mydas Pacific shape-based morphotypes (67 SNPs; electronic supplementary material, table S2). (a) Venn diagram of loci under selection detected by three approaches used in this study, and (b) genes associated with putative main categories based on the descriptions related to each of the GO terms (see electronic supplementary material, table S3 for more detail). (Online version in colour.)

The overrepresentation analysis results showed that GO categories linked to protein binding, energetic metabolism, ion binding, glycolipid biosynthesis, epidermis development, embryogenesis/morphogenesis, renal morphology and function, cardiovascular performance, muscle activity and metabolism, nervous system, reproductive and endocrine system, behaviour, cellular processes and intracellular transport were significantly overrepresented (figure 4; electronic supplementary material, table S3). No KEGG or Reactome categories were significantly overrepresented. The REViGO results are shown in figure 4b, and the categories generated by this analysis are included in electronic supplementary material, table S3.

Figure 4.

Figure 4. Enrichment analysis using outlier SNPs recovered among Chelonia mydas Pacific shape-based morphotypes (electronic supplementary material, table S2). (a) g:GOSt multiquery Manhattan plot that illustrates the enriched terms. The x-axis shows the functional terms, and the corresponding enrichment p-values in negative log10 scale are illustrated on the y-axis. Each circle on the plot represents a single functional term. The circles are colour-coded by data source and size-scaled according to the number of annotated genes in that term (electronic supplementary material, table S3). (b) Interactive graphic view of REVIGO. Circle size indicates the frequency of the GO term in the underlying database. Similar GO terms are closer and highly similar GO terms are linked by edges in the graph. Circle colour is according to putative general categories (electronic supplementary material, table S3), and (c) bar charts showing the proportion of the enriched terms grouped according to the putative general categories included in the electronic supplementary material, table S3 (see Material and methods for more detail). GO:MF, molecular function; GO:BP, biological process, GO:CC, cellular component. (Online version in colour.)

4. Discussion

While the existence of reproductive isolation is logistically difficult to study in migratory, long-lived and endangered species, our results using about 9500 SNPs support the existence of genetic divergence between C. mydas Pacific shape-based morphotypes, demonstrating their independent evolutionary trajectory. The marked genetic structure between morphotypes found here is probably due to the natal homing behaviour of this species, the most widespread hypothesis of population connectivity in sea turtles [10]. Considerable literature using molecular markers with variable mutation rates (mtDNA, microsatellites and SNPs [7,9,4749]) in accordance with our results supports a strong natal homing behaviour by C. mydas females and low male-mediated gene flow in the Pacific Ocean. The lack of genetic structure among FGs (electronic supplementary material, figure S2b) demonstrates the multiple origins of turtles inhabiting these key habitats, which is corroborated for the first time by analysing thousands of SNPs distributed across the genome.

Genomic divergent regions between shape-based morphotypes linked to numerous biological functions were also found (figure 3; electronic supplementary material, table S3). However, more importantly, we report genes associated with pigmentation, morphogenesis, osmoregulation, digestive metabolism, reproductive system and the colder environment of the EPO. Although most of these regions are situated in introns, there is considerable evidence supporting that introns play a fundamental and evolutionarily conserved role in regulating eukaryotic gene expression. In fact, the effect of an intron on gene expression can be larger than that of the promoter of the same gene, and in many cases genes with an intact promoter are essentially not expressed at all without an intron [50]. Similarly, synonymous mutations (nucleotide changes that do not alter the encoded amino acid) have often been assumed to be neutral or silent [51]. However, recent evidence has demonstrated this type of substitution brings about changes in the phenotype by affecting splicing accuracy, translation fidelity, mRNA structure and protein folding [51,52].

Air-breathing marine vertebrates have diving responses characterized by circulatory adjustments (bradycardia, peripheral vasoconstriction and reduction in cardiac output) that redistribute the blood flow to oxygen-sensitive organs, while other tissues become anaerobic [53,54]. In prolonged dives (e.g. foraging diving), these changes may lead to peripheral hypoxia and metabolic acidosis. Such physiological responses can be accentuated in low water temperatures that occur in the EPO, and ectotherms such as sea turtles can be strongly affected [54]. Here we found a synonymous mutation in a CDS of the CREBBP gene which is associated with the response to hypoxia (GO:0001666). We also found genes related to thermoregulation (heat shock: FKBP5, GO:0031072; heat production: DNAH8, GO:0010286 and regulation of cold-induced thermogenesis: PLCL1, GO: GO:0120161). Several biological functions related to cardiovascular performance (e.g. heart development and blood pressure), muscle activity and metabolism and energetic metabolism were also enriched (figure 4; electronic supplementary material, table S3). These findings are probably associated with the dive response of these ectotherms to the cold EPO waters [53,54]. Such results could even be linked to overwintering torpor, which has been reported at water temperatures below 15°C in C. mydas populations in this region [16,20].

Melanin-based coloration in reptiles is highly variable and plays crucial roles in thermoregulation, camouflage, UV protection and sexual selection [55,56]. We found differences in allele frequency in EDA, MEF2C, CCDC88A, KRT5 and CREBBP genes between C. mydas shape-based morphotypes. These genes are involved in pigmentation (GO:0043473), ectoderm differentiation (GO:0010668), epidermis development (GO:0008544) and cellular response to ultraviolet (UV) radiation (GO:0009411) (electronic supplementary material, table S2). MEF2C is associated with melanocyte development [57] and EDA with skin patterning and early morphogenesis of all ectoderm-derived amniotic structures, including bird and crocodile scales, turtle scutes, mammalian hairs, mammary glands and avian feathers [58]. Experimental studies have demonstrated that mutations in EDA correlate with coat colour phenotype in mice [59,60]; this gene has been proposed as a pigmentation candidate gene showing positive signatures of selection in humans [61]. Variations in skin colour in human populations are adaptive, related to the regulation of ultraviolet radiation penetration in the integument [62] and regulated by a large number of genes [63]. Our findings, especially the substitution in the CREBBP gene (linked to the response to UV radiation), open new questions on the underlying basis of melanism in C. mydas populations deserving further research.

The marine habitat imposes an osmoregulatory challenge for vertebrates that maintain body fluids hyposmotic to the environment [64]. Marine reptiles including turtles have evolved specialized kidneys and salt glands that allow them to be efficient in the utilization of water and salt elimination [64]. Genes (MEF2C, DCHS1, WNK4 and ESR2; electronic supplementary material, table S2) and enrichment functions associated with renal morphology and function (e.g. nephron and renal tubule development, kidney epithelium development; electronic supplementary material, table S3) are reported and probably linked to osmoregulation, as expected for populations inhabiting in the high salinity conditions of the EPO. Several SNPs were also linked to genes related to the development and morphogenesis of anatomical structures; two were situated in CDS (CREBBP and ZNF384; electronic supplementary material, table S2). The number of divergent regions (introns and CDS) found, together with the enrichment of biological functions related to morphogenesis (figure 4; electronic supplementary material, table S3), point to an adaptive divergence previously hypothesized in this species [9,12].

Unlike the yellow morphotype, black turtles seem to have a predominantly omnivorous diet [14,16]. Although biological functions linked to digestive metabolism were not enriched, we found genes associated with lipid metabolism (HSPG2, PLCL1, MEF2C and CREBBP: GO:0006629, GO:0019915, GO:0019216), carbohydrates (MEF2C and CRYBG3: GO:0052576, GO:0030246), proteins (ERMP1 and EDA: GO:0006508), glucose (WNK4, CACNA1C, HDAC4 and MEF2C: GO:0042593, GO:0045820, GO:0071333), folic acid (EDA: GO:0046655), fatty acid omega (EDA: GO:0010430) and fat-soluble vitamins (MEF2C: GO:0033197) (electronic supplementary material, table S2), suggesting diet adaptation of the EPO morphotype. Variations in feeding ecology have been reported between black and yellow turtles in Pacific FGs [14,65], aligning with our findings.

Estimation of the age at first reproduction varies widely among C. mydas populations. Black turtles mature at smaller sizes than conspecifics in the northwestern Atlantic, Hawaii and Australia [8]. Females from this region also lay fewer eggs and show different reproductive periodicity, which has been attributed to the widely variable nutrient availability in the EPO [16,19]. Enriched categories linked to the reproductive system, endocrine system and behaviour found here (e.g. female gonad development, female sex differentiation and response to parathyroid hormone; figure 4; electronic supplementary material, table S3) constitute a baseline for future research on these morphotypes.

One limitation of a reduced-representation genome approach is the possibility that an anonymous locus showing unusual ‘outlier’ behaviour is not under selection itself, but is instead in LD with the selected site (or sites) [66]. Likewise, how complete and accurate a genome assembly is will determine its suitability for posterior analyses [67]. Although our approach is a genome-reduced representation and the C. mydas draft genome currently available is very incomplete, neutral and under putative selection SNPs reported here provide novel results on genetic structure and adaptation signatures of this species in the EPO.

Candidate genes related to cold adaptation, melanism, morphogenesis, osmoregulation, diet and reproduction are proposed. Nevertheless, our results showed similar past demography history in both morphotypes, with a population decline over 80 kyr ago (electronic supplementary material, figure S3). This bottleneck could have intensified the effect of drift in both groups, increasing their differences; thus, adaptive signal observed here must be interpreted cautiously. Future analyses based on a high-quality reference genome, together with associations or correlations between markers and traits or environmental variables, will be better evidence for adaptive significance [68]. Comparative genomic analyses will also be useful to examine in further detail genes controlling physiological or morphological sea turtle traits related to EPO adaptations, as has been studied in other Chelonians and seabirds associated with desert and climate adaptations, respectively (e.g. genes related to desert adaptations in the Gopherus species complex [69] and genes related to climate adaptation in penguins [70]). Our study challenges researchers to increase adaptive investigation in marine vertebrate species which inhabit this unique ecosystem in the world.

Ecological, genetic and morphological data have led to suggesting the C. mydas EPO populations as a regional management unit (RMU) [6] and a distinct population segment (DPS) [8]. The black turtle has even been proposed as a full species (Chelonia agassizii; Bocourt 1868) [13] or subspecies (C. mydas agassizii) [18]. Although the taxonomic evaluation of C. mydas requires extensive sampling including the entire distribution range of the species, our results based on neutral and adaptive diversity suggest that the black turtle should be managed as an independent evolutionary unit, and specific conservation strategies for each shape-based morphotype should be established (see more details in [9]).

Finally, our research contributes to disentangling the mystery of black turtle adaptation and also challenges us to think about how cosmopolitan marine species adapt to changing conditions and how can we protect their genetic reserves and adaptive potential.

Data accessibility

The authors declare that all data supporting the findings of this study are available within the paper and its electronic supplementary material [71]. The datasets generated and analysed during the current study are available at https://doi.org/10.34691/FK2/RJZUFC [72].

Authors' contributions

R.Á.-V.: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing-original draft, Writing-review and editing; N.R.-H: Data curation, Formal analysis, Methodology, Software, Validation, Writing-review and editing; M.H.: Data curation, Funding acquisition, Investigation, Resources, Writing-review and editing; C.R.: Conceptualization, Investigation, Methodology, Software, Supervision, Validation, Writing-review and editing; H.A.B.: Conceptualization, Investigation, Methodology, Resources, Supervision, Writing-review and editing; R.A.-D.: Formal analysis, Methodology, Software, Validation, Visualization, Writing-review and editing; E.R., M.L., D.A.G., J.P.M.-P., D.E.A.-R., J.A.-S., C.O.-A. and J.C.M.: Funding acquisition, Investigation, Resources, Writing-review and editing; J.A.V.: Conceptualization, Investigation, Supervision, Validation, Writing-review and editing; D.V.: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Resources, Software, Supervision, Validation, Visualization, Writing-original draft, Writing-review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Competing interests

We declare we have no competing interests.

Funding

This work was supported by CONICYT PhD Scholarship no. 21160168, Núcleo Milenio de Ecología y Manejo Sustentable de Islas Oceánicas (ESMOI) and funding sources from each associated researcher.

Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5477735.

Published by the Royal Society. All rights reserved.

References