Identifying genetic lineages through shape: An example in a cosmopolitan marine turtle species using geometric morphometrics

The green turtle (Chelonia mydas) is a globally distributed marine species whose evolutionary history has been molded by geological events and oceanographic and climate changes. Divergence between Atlantic and Pacific clades has been associated with the uplift of the Panama Isthmus, and inside the Pacific region, a biogeographic barrier located west of Hawaii has restricted the gene flow between Central/Eastern and Western Pacific populations. We investigated the carapace shape of C. mydas from individuals of Atlantic, Eastern Pacific, and Western Pacific genetic lineages using geometric morphometrics to evaluate congruence between external morphology and species’ phylogeography. Furthermore, we assessed the variation of carapace shape according to foraging grounds. Three morphologically distinctive groups were observed which aligned with predictions based on the species’ lineages, suggesting a substantial genetic influence on carapace shape. Based on the relationship between this trait and genetic lineages, we propose the existence of at least three distinct morphotypes of C. mydas. Well-defined groups in some foraging grounds (Galapagos, Costa Rica and New Zealand) may suggest that ecological or environmental conditions in these sites could also be influencing carapace shape in C. mydas. Geometric morphometrics is a suitable tool to differentiate genetic lineages in this cosmopolitan marine species. Consequently, this study opens new possibilities to explore and test ecological and evolutionary hypotheses in species with wide morphological variation and broad geographic distribution range.

Introduction Previous evidence indicates the carapace shape in chelonians possesses a significant heritable genetic component [30,31,32], and some studies have reported an association between phylogeographic differentiation and shell shape variation both tortoises and freshwater turtles [32,33,34].
Given the evolutionary history of C. mydas in the Atlantic and Pacific basins, and the heritability of this trait, in our study we expect to find congruence between carapace shape variation and genetic lineages of this species. We hypothesized the presence of three morphological distinctive groups: Atlantic Group, Eastern Pacific Group and Western Pacific Group. Here, we examined the morphological variation of C. mydas according to genetic lineages or natal origins (Atlantic, Eastern Pacific and Western Pacific) and also according to foraging grounds (Uruguay, Costa Rica, Galapagos-Ecuador, Chile and New Zealand).

Ethic statements
All procedures involving animals were carried out in accordance with approved guidelines and protocols under permits issued by national agencies from all countries involved in this study. Details on permits are described below.

Study area and data collection
Our study included C. mydas foraging grounds located in the South Western Atlantic (Uruguay), Eastern Pacific (from north to south: Costa Rica, Galapagos-Ecuador and Chile) and South Western Pacific regions (New Zealand). Specific locations are described below: Uruguay. The Coastal-Marine Protected Area of Cerro Verde (CMPA; 33.93˚S, 53.50W ) is located in north-eastern Uruguay. Juvenile turtles (n = 197) were captured alive using nets during the warmer months (December to April) between 2012 and 2016. The capture method is described in [35]. Captures were conducted by Karumbé NGO technicians and were authorized by the Fauna Department-Ministry of Cattle, Agriculture and Fishing of Uruguay (license no. 073/08, 323/11, 12/14), and the Fauna Division-Ministry of Housing, Territorial Planning and Environment of Uruguay (DF 141/16).
Costa rica. Matapalito Bay (10.93˚N, 85.78˚W) is a small inlet of the Santa Elena Peninsula located in north-western Costa Rica. Juveniles and adult turtles (juveniles, n = 22; adults, n = 20) were captured using nets between 2012 and 2017 as described in [36]. Field work and sampling were authorized by the Guanacaste Conservation Area (ACG) of the Ministry of Environment, Energy and Telecommunications (MINAE).
Galapagos. The Galapagos Archipelago (~00.66˚S, 90.55˚W) is a group of islands of volcanic origin located about 1,000 km west of mainland Ecuador. Juvenile and adult turtles n = 74; and n = 5, respectively) were captured with nets and by hand in the Galapagos foraging grounds between 2004 and 2005. Capture methods are described in [8]. Research permits were provided by the Galapagos National Park.
Chile. Bahía Salado (27.68˚S, 70.98˚W) is a bay located in the Atacama Region in northern Chile (mainland Chile). Juvenile turtles (n = 9) were captured using nets during spring 2014 and summer 2018. The capture method is described in [37]. Captures were authorized by the Chilean Sub-Secretariat of Fishing (SUBPESCA, by its Spanish abbreviation), through a Research Capture Permit granted in April 2013 (Exempt Resolution no. 917) and renewed in July 2014.
New Zealand. The study area encompassed the inshore waters of North Island, extending from Cape Reinga (34.41˚S, 174.66˚E) south to Opotiki, in the Bay of Plenty (37.98˚S, 177.28˚E). Live and dead stranded juvenile turtles (n = 25) were collected for examination between 2012 and 2017. All individuals included in this study exhibited good body condition and the cause of stranding or death was not associated with an underlying pathology (most died by interaction with fishing gear). In addition, all dead individuals were fresh (with presence of eyes, all head scales and with carapace scutes intact) [38]. Collection specifications are described in [9]. Collections were authorized by the Department of Conservation of New Zealand (authorization no. AK 30931-FAU and 52128-FAU).

Genetic lineage determination
In order to differentiate genetic lineages in Pacific C. mydas foraging grounds, mtDNA control region haplotypes were identified for each individual using primers LCM15382 (5'-GCTTA ACCCTAAAGCATTGGO3') and H950g (5'GTCTCGGATTTAGGGGTTTGO3') designed by [39]. These data were obtained according to the project's specific frameworks undertaken in each country in collaboration with external researchers (Costa Rica, Heidemeyer et al. unpublished data; Galapagos, Dutton & Zarate unpublished data, Chile, [37]; New Zealand, [40]). Thus, for the purposes of this paper, only the lineage (natal origin) corresponding to the Eastern or Western Pacific was indicated (not the specific haplotype; S1 Table). Two divergent evolutionary lineages for C. mydas have been described in the Atlantic/Mediterranean region, the "northern lineage" and the "southern lineage" [41][42][43]. The southern lineage encompasses the eastern Caribbean, South Atlantic and West African rookeries [41]. Given that about 90% of green turtles foraging in Uruguayan waters have their natal origin in the South Atlantic region [44], turtles from Uruguay of this study were classified into "Atlantic southern lineage".

Shape analysis
Geometric morphometric analyses included 352 C. mydas individuals from five foraging grounds: 197 with Atlantic natal origin (Atlantic genetic lineage-AGL, Uruguay); 105 with Eastern Pacific natal origin (Easter Pacific genetic lineage-EPGL, Costa Rica, Galapagos, Chile and New Zealand) and 50 with Western Pacific natal origin (Western Pacific genetic lineage-WPGL, Costa Rica, Galapagos and New Zealand) (S1 Table). GM analysis focused on variation in carapace shape, and was performed using dorsal photographs of individuals of distinctive genetic lineages and foraging grounds. All photographs were obtained using a reference scale. Thirty-six landmarks (Fig 1) were digitized with TPS Dig 2.30 software [45]. Landmarks were obtained between specific carapace scutes and at the borders of the marginal scutes, projected from the closest lateral scute [29]. A Procrustes superimposition was applied to the landmark data in order to remove any non-shape elements.
A multivariate regression was carried out to determine the influence of size on shape (allometry) in the dataset using centroid size (size variable) as an independent variable and shape (Procrustes coordinates) as a dependent variable [46]. Furthermore, a permutation test using 10,000 iterations was performed to assess the significance of the influence of the size on shape.
Principal component analyses (PCA) were performed using the covariance matrices of shape variation and the average shape variation in genetic lineages and foraging grounds. To visualize shape average changes and their distribution in the averaged shape space a PCA scatterplot was performed. Additionally, in order to visualize the variation in carapace shape, the average carapace shape was rendered for each genetic lineage and foraging ground.
A canonical variate analysis (CVA) was performed to have a better graphical representation of the data and to discriminate groups based on carapace shape variation in different genetic lineages and foraging grounds. The CVA is a multivariate statistical method used to find the shape characters that best distinguish among groups of specimens. The results were reported as Procrustes distances and the respective p-values for these distances, after permutation tests (10,000 iterations).
A Procrustes ANOVA was carried out to assess the significance of the differences in carapace shape between genetic lineages and between foraging grounds. All analyses were performed using MorphoJ software [47]. For these analyses, data were pooled according to genetic lineage (Atlantic, n = 197; Eastern Pacific, n = 105 and Western Pacific, n = 50), and foraging ground (Uruguay, n = 197; Costa Rica, n = 42; Galapagos, n = 79; Chile, n = 9 and New Zealand, n = 25) (S1 Table).

Carapace shape variation according to C. mydas genetic lineages
Multivariate regression showed a 21.6% of allometry with a significant permutation value (p-value = <0.0001) (Fig 2). Thus, a correction for allometry was performed and all the shape analyses (PCA and CVA) were carried out using the covariance matrix of the data corrected by size (data used from the residual of the multivariate regression). Given this allometric correction, a distinction between juveniles and adults was not performed in this study.
The first three principal components (PC) accounted for 69.7% of shape variation (PC1 = 43.6%; PC2 = 13.4% and PC3 = 12.64%). The scatterplot of the PCA (S1 Fig) showed a central cloud of points from the AGL, in contrast to a higher dispersion of points between the two Pacific lineages (EPGL and WPGL). Average carapace shape based on genetic lineage clearly varied between groups ( Fig 3A). The AGL exhibited a wider carapace, while, in contrast, an oval and triangular carapace shape were observed in the WPGL and EPGL, respectively. Differences in the second lateral scute were identified (landmarks 10-13, 15, 16 and 18- then has played an important role in the divergence of Pacific and Atlantic clades of different marine turtle species including C. mydas [4,48,49]. In the Pacific Ocean, the East Pacific Barrier (EPB), a 5000 to 8000 km deep-water extension located east of Hawaii, has been described as one of the most important barriers for dispersal separating Eastern Pacific biota from the Central Pacific and Indo West-Pacific regions [50,51]. Nevertheless, a recent study of C. mydas populations suggested that the Pacific region west of Hawaii has been a more significant barrier to gene flow than the EPB, and that the split between Central/Eastern and Western lineages in this species occurred about 340,000 years ago [6]. The uplift of the Panama Isthmus precedes the divergence in C. mydas described by Dutton et al. (2014) [6] in the Pacific Ocean (Eastern and Western Pacific lineages), and this is congruent with the degree of morphological differentiation observed in this study, where differences were more evident between basins (Atlantic-Pacific) than within the Pacific basin.
Natal homing behaviour has been demonstrated in most marine turtle species using mtDNA sequencing [2]. These maternally inherited markers show strong population structure among nesting colonies while nuclear loci reveal a contrasting pattern of male-mediated gene flow [2]. Particularly in C. mydas, studies confirm this reproductive behaviour for females and males; however, the geographic specificity of homing is uncertain, and it may vary for hundreds of kilometres among different populations [5,6,52]. In the Pacific Ocean, our results showed two distinctive morphological groups that were consistent with genetic lineages (EPGL and WPGL). This aligned with the natal homing theory that states turtles returning to their region of origin for mating and nesting, which influences the population's genetic structure. In summary, our study showed a parallel between carapace shape variation and the evolutionary history of C. mydas, as initially predicted, associated with a geographic barrier limiting gene flow between ocean basins (Panama Isthmus). In addition, it is possible this association may also be influenced by life history traits (i.e. natal homing) and oceanographic conditions (e.g. barrier west of Hawaii) that differentiate populations within the Pacific Ocean. Likewise, the congruence between phylogeography and morphology here, suggests a significant genetic influence on the carapace shape in C. mydas as reported in other chelonians [30,31,32].

Association between grouping based on carapace shape and morphotypes of C. mydas
The East Pacific form of C. mydas tends to be distinguished by conical carapace shape and dark coloration [7,15,17]. Although genetic data do not support the evolutionary distinctiveness of the black morphotype, a population level-differentiation exists between this form and the lighter form (light morphotype) [4,53]. Kamezaki & Matsui (1995) [13] examined geographic variation of skull morphology in C. mydas and they observed an exclusive distinction of black turtles (Galapagos nesting population), suggesting that due to its isolation, this group contains unique morphological characteristics. Later, Okamoto & Kamezaki (2014) [54] studied C. mydas foraging grounds in Japan (Pacific Ocean) and reported differences in carapace shape between black and light morphs, with a narrowing at the level of the eleventh marginal scute in black turtles, which remained consistent throughout growth. In this work, using geometric morphometrics, we evaluate the carapace shape of C. mydas and we observed that turtles with a Western Pacific origin (WPGL, putative "light morph"), exhibited an oval or more elongated carapace, while turtles with an Eastern Pacific origin (EPGL, "black morph") exhibited a triangular (conical) carapace. These results are consistent with the visual descriptions of both morphotypes [7,[15][16][17]. Moreover, turtles from the EPGL exhibited a narrowing of the second lateral scute and an elongation of the last central scute (Fig 3A). On the other hand, grouping based on genetic lineages was found here, differentiating the Atlantic group, and two groups of the Pacific (Fig 3B). Such grouping did not correspond with the worldwide recognized assignation of morphotypes [7,15,17], due to our data showed the presence of a distinctive Atlantic morphotype that differs from the light morphotype that occurs in the Pacific Ocean (Fig 3).
Until now, differences in carapace shape between C. mydas from the Atlantic and Pacific lineages and intra-Pacific lineages had not been tested. Colour alone should not be considered as a diagnostic tool to distinguish between populations because this character is highly variable throughout the range of C. mydas [9,55]. In fact, for this reason we did not include this character in our analysis.   Our results show that carapace shape could enable us to differentiate intraspecific genetic lineages in this cosmopolitan species. Based on these results, we propose the existence of at least three distinct morphotypes: Atlantic, Eastern Pacific and Western Pacific. Nevertheless, further research incorporating other evolutionary lineages (e.g. "northern lineage" from the Atlantic/Mediterranean) may provide more insight into carapace shape variation and the designation of other morphotypes, globally.

Carapace shape variation through foraging grounds in the Pacific Ocean: Conservation implications
Green turtles have a circumglobal distribution with hundreds of nesting beaches and foraging grounds making up a complex network of migratory routes [1]. As described previously, a marked genetic differentiation has been observed between green turtle populations at a larger scale in the Pacific Ocean (corresponding to EPGL and WPGL [6]), which is most likely associated to oceanographic conditions and natal homing behavior. At a finer scale, using mitochondrial markers, genetic structuring has also been observed, which has been used to identify distinctive management units (MUs) [6,10]. Management units or stocks correspond to populations that exchange so few migrants that are genetically distinct and demographically independent [56]. Specifically, in the Central and Eastern Pacific region, five MUs have been designated: Northwest Hawaii, Revillagigedo, Michoacan, Costa Rica, and Galapagos-Machalilla [6,10]. Most of these MUs have been proposed based on genetic data exclusively from nesting sites [6]. Just for the case of Galapagos and Northwest Hawaii genetic data from foraging grounds have also been considered [10,57].
Our results based on foraging grounds showed well differentiated groups for Costa Rica, Galapagos and New Zealand. Thus, although the natal origin (EPGL and WPGL) had great relevance in the morphological differentiation of populations, as previously discussed, foraging sites also could have an effect on carapace shape variation in this species.
Although empirical fitness data would be required to properly asses the adaptive value of the carapace shape in C. mydas' lineages, further research examining the relationship between carapace morphology, specific environmental conditions in foraging grounds, and non-neutral genetic variation, may provide more insight about selection pressures on the carapace shape of this species. In this context, and given the longevity and conservation status of C. mydas that restrict experimental studies, genomic tools may be useful to address these questions.
Regarding Chile, the lack of population differentiation could be due to Chilean waters just constituting foraging habitats for C. mydas (no nesting exists) and the natal origin of individuals is mainly Galapagos-Ecuador [37]. The latter, is supported by the low differentiation between both populations (Chile-Galapagos, Fig 4 and Table 3). In any case, research increasing the sample size in Chile could allow us to observe a clearer pattern of morphological variation.
Our results based on foraging grounds reveal the importance of studying the role of selection on the morphology of C. mydas, and the relevance to incorporate environmental and genetic information from these habitats when a MU is defined. In this way, by integrating data from the key habitats of the green turtle's life cycle, the evolutionary potential of their threatened populations can be preserved.

Conclusions
Our study shows that the carapace shape in C. mydas is markedly associated with the species' lineages suggesting a substantial genetic influence on this trait. Based on the relationship between carapace shape and genetic lineages found here, we propose the existence of at least three distinct morphotypes of C. mydas: Atlantic, Eastern Pacific and Western Pacific. Welldifferentiated groups in some foraging grounds may suggest an effect of ecological or environmental operating conditions on morphological variations of C. mydas. Likewise, these results highlight the importance to integrate data from rookeries and foraging grounds to define MUs in order to conserve the evolutionary potential of distinctive populations. This is the first study using geometric morphometrics to evaluate the congruence between phylogeography and morphological variation in marine turtles. Our results, based on this emergent tool, open new possibilities to test ecological and evolutionary hypotheses in morphologically variable and widely distributed species.
Supporting information S1 Table. List of individuals used in this study including foraging ground (country), specific location, year of collection and genetic lineage or natal origin (haplotype origin).